Python GDAL – How to Loop for Creating a Multi-Page TIFF (Layer Stack)

gdalpython

So far I use python to create a layer stack for each pair of TIF files now l am looking for the solution to loop this process.

This is my code for creating a layer stack.

import os
from osgeo import gdal
imgB01Path = "E:/eyerStack/Raw2Tif_2016/tif16/2016.01.01.Smb01.tif"
imgB02Path = "E:/eyerStack/Raw2Tif_2016/tif16/2016.01.01.Smb02.tif"
outImgPath = "E:/LeyerStack2016/2016.01.01.tif"

oB01Band = gdal.Open(imgB01Path)
oB02Band = gdal.Open(imgB02Path)

gDriver = oB01Band.GetDriver()
outImg = gDriver.Create(outImgPath, oB01Band.RasterXSize,oB01Band.RasterYSize, 2, gdal.GDT_UInt16, ["PHOTOMETRIC=RGB"])

outImg.SetGCPs(tuple(oB01Band.GetGCPs()),oB01Band.GetGCPProjection()) #Set GCPs
outImg.SetGeoTransform(tuple(oB01Band.GetGeoTransform())) #Set GeoTransform
outImg.SetProjection(oB01Band.GetProjection()) #Set Projection

B01 = oB01Band.GetRasterBand(1).ReadAsArray()
B02 = oB02Band.GetRasterBand(1).ReadAsArray()

outImg.GetRasterBand(1).WriteArray(B01)
outImg.GetRasterBand(2).WriteArray(B02)

oB01Band = None #Clear Memory
oB02Band = None #Clear Memory
outImg = None #Clear Memory

The example file name of band1

2015.01.01.Smb01.tif
2015.01.09.Smb01.tif
2015.01.17.Smb01.tif
2015.01.25.Smb01.tif
2015.02.02.Smb01.tif

The example file name of band2

2015.01.01.Smb02.tif
2015.01.09.Smb02.tif
2015.01.17.Smb02.tif
2015.01.25.Smb02.tif
2015.02.02.Smb02.tif

This is the complete adapted code.

import os
from osgeo import gdal

tifFolder = r'E:/ProcessNdvi2015_2021/05LayerStack/tif/'

tifs = [os.path.join(root, file) for root, folder, files in os.walk(tifFolder) for file in files if file.endswith('tif') and 'smb' in file.lower()]
tifs.sort()
count = 0
for t1, t2 in zip(tifs[::2], tifs[1::2]):
    ##print(t1,t2)
    oB01Band = gdal.Open(tifs[count])
    oB02Band = gdal.Open(tifs[count+1])
    filename = tifs[count].split("/")[-1]
    
    outImgPath = "E:/ProcessNdvi2015_2021/05LayerStack/LayerStack2/" + filename 

    gDriver = oB01Band.GetDriver()
    outImg = gDriver.Create(outImgPath, oB01Band.RasterXSize,oB01Band.RasterYSize, 2, gdal.GDT_Int16, ["PHOTOMETRIC=RGB"])

    outImg.SetGCPs(tuple(oB01Band.GetGCPs()),oB01Band.GetGCPProjection()) #Set GCPs
    outImg.SetGeoTransform(tuple(oB01Band.GetGeoTransform())) #Set GeoTransform
    outImg.SetProjection(oB01Band.GetProjection()) #Set Projection

    B01 = oB01Band.GetRasterBand(1).ReadAsArray()
    B02 = oB02Band.GetRasterBand(1).ReadAsArray()

    outImg.GetRasterBand(1).WriteArray(B01)
    outImg.GetRasterBand(2).WriteArray(B02)

    oB01Band = None #Clear Memory
    oB02Band = None #Clear Memory
    outImg = None #Clear Memory
    count = count+2

Best Answer

I like os walk in combination with collections defaultdict list for problems like this. You can store each date as a key, and the two files as values:

import os
from collections import defaultdict

tifFolder = r'/home/bera/Desktop/GIStest/input_files/'

tifs = defaultdict(list)

for root, folders, files in os.walk(tifFolder): #Find all tif files with smb in filename in the tifFolder
    for file in files:
        if 'smb' in file.lower() and file.endswith('tif'):
            tifs[file.split('.S')[0]].append(os.path.join(root, file)) #Add the full path to the file as value and date as key to the defaultdict
           
#defaultdict(list,
            #{'2015.01.01': ['/home/bera/Desktop/GIStest/input_files/2015.01.01.Smb01.tif', '/home/bera/Desktop/GIStest/input_files/2015.01.01.Smb02.tif'],
             #'2015.01.09': ['/home/bera/Desktop/GIStest/input_files/2015.01.09.Smb01.tif', '/home/bera/Desktop/GIStest/input_files/2015.01.09.Smb02.tif'],
             #'2015.01.17': ['/home/bera/Desktop/GIStest/input_files/2015.01.17.Smb02.tif', '/home/bera/Desktop/GIStest/input_files/2015.01.17.Smb01.tif']})

for date, tiffiles in tifs.items():
    smb01, smb02 = sorted(tiffiles) #You will get ValueError here if there's missing files
    print(date)
    print(smb01)
    print(smb02)
    #oB01Band = gdal.Open(smb01)
    #oB02Band = gdal.Open(smb02)
    # ...
    

#2015.01.01
#/home/bera/Desktop/GIStest/input_files/2015.01.01.Smb01.tif
#/home/bera/Desktop/GIStest/input_files/2015.01.01.Smb02.tif

#2015.01.09
#/home/bera/Desktop/GIStest/input_files/2015.01.09.Smb01.tif
#/home/bera/Desktop/GIStest/input_files/2015.01.09.Smb02.tif

#2015.01.17
#/home/bera/Desktop/GIStest/input_files/2015.01.17.Smb01.tif
#/home/bera/Desktop/GIStest/input_files/2015.01.17.Smb02.tif

You can also try something like this, which is easier:

import os

tifFolder = r'C:\Users\jksdhfd\Desktop\Temp'

tifs = [os.path.join(root, file) for root, folder, files in os.walk(tifFolder) for file in files if file.endswith('tif') and 'smb' in file.lower()]
tifs.sort()

for t1, t2 in zip(tifs[::2], tifs[1::2]):
    print(t1, t2)

Related Solutions

[GIS] How to create tiles from tiff file using python

The individual image tiles which gdal2tiles.py creates and stores into directory structure are not georeferenced. Therefore for example gdalinfo sees that every tile has its origin at (0 0). However, the auxiliary tilemapresource.xml contains all the information that is needed for georeferencing.

There is a TileLayer plugin for QGIS which can interpret the tilemapresource.xml file and show your tiles correctly. Find details about how to use it with gdal2tiles from gdal2tiles and qgis.

Python – Making Loop for Reading Multiple Rasters Using Python and Folium

Here is an approach for looping.

Environment

Windows 10
ipynb

Methods

Since gdal can only read ISG files I cannot serve testdata. So I made testdata with Geotif (.tif). Change that filetype to ISG, and also the path of the files for your purposes.

If you dont have ipynb, then the folder needs to be a link, where are available files and there have to be a writeable folder in /var/www/ by usergroup wwwdata for the converted files (set crs to openstreetmap crs).

According to this, map location cannot be changed after creation. But maps cannot be added if map is notbcreated. So if you want get for example the mid point of all layers for location, you have to iterate the rasters first and calculate the mid point.

I tested your code for bounding box in another crs. Here it wasnt working correctly. Maybe I missed something. So I decided to append another gdal method from here, which gave me expected results.

Related to this I think folium crs is 3857 for all I/O. So I added another gdal method to change crs from here.

The layer names in the sidebar cannot be changed within the function the adds the layer to the map I think. So I added some generic declaration.

Code

import folium
from osgeo import gdal, osr
import os

def GetExtent(ds):
    """ Return list of corner coordinates from a gdal Dataset """
    xmin, xpixel, _, ymax, _, ypixel = ds.GetGeoTransform()
    width, height = ds.RasterXSize, ds.RasterYSize
    xmax = xmin + width * xpixel
    ymin = ymax + height * ypixel

    return (xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin)

def ReprojectCoords(coords,src_srs,tgt_srs):
    """ Reproject a list of x,y coordinates. """
    trans_coords=[]
    transform = osr.CoordinateTransformation( src_srs, tgt_srs)
    for x,y in coords:
        x,y,z = transform.TransformPoint(x,y)
        trans_coords.append([x,y])
    return trans_coords

map = folium.Map(location=[49.51, 11.23],default_zoom_start=14,crs='EPSG3857')

def run_exec(name,fnum,path,temp_files_path,layername):
    input_raster = gdal.Open(path+'/'+name)
    output_raster = temp_files_path+'/' + name

    #change epsg to openstreetmap epsg
    gdal.Warp(output_raster,input_raster,dstSRS='EPSG:3857')
    nl = gdal.Open(temp_files_path+'/' + name)
    
    #get layer extent
    ext=GetExtent(nl)
    src_srs=osr.SpatialReference()
    src_srs.ImportFromWkt(nl.GetProjection())
    tgt_srs = src_srs.CloneGeogCS()
    geo_ext=ReprojectCoords(ext, src_srs, tgt_srs)

    #Get data
    nlband = nl.GetRasterBand(1)
    nldataset= nlband.ReadAsArray(0,0,nl.RasterXSize,nl.RasterYSize)

    exec("""
#variables must be generic as layer_name cannot be set in following function
nl"""+ str(fnum) +""" = folium.raster_layers.ImageOverlay(
    image=nldataset,
    #bounds=[[nlymin, nlxmin], [nlymax, nlxmax]],
    bounds=geo_ext,
    opacity=1,
    show=False,
).add_to(map)

nl"""+ str(fnum) +""".layer_name = layername
    """)

#paths in ipynb, layernames for map
#order of layernames is lexographically to filenames in path
path="./0409"
temp_files_path="./1028"
layernames=["layer1","layer2"]

fnum=0
for root, dirs, files in os.walk(path, topdown=False):
    for name in files:
        if ".tif" in name:
            layername=layernames[fnum]
            fnum+=1
            run_exec(name,fnum,path,temp_files_path,layername)


folium.LayerControl().add_to(map)
map

Best Answer

Related Solutions

[GIS] How to create tiles from tiff file using python

Python – Making Loop for Reading Multiple Rasters Using Python and Folium

Related Question