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

foliumgdalpython

I am a newbie in programming and wanted to ask if someone could help me with example of loading/reading multiple rasters using for loop (or maybe another alternatives to deal with multiple files?). Basically, it loads several rasters with gdal, then plot and export the map as html file using folium. Here is the code :

import gdal
import folium

#poland dataset
pol = gdal.Open('ISG/POLAND_KTH-PL-GEOID2015_15_G_20160729.isg') 

#Get coordinates, cols and rows
polgeotransform = pol.GetGeoTransform()
polcols = pol.RasterXSize 
polrows = pol.RasterYSize 

#Get extent
polxmin=polgeotransform[0]
polymax=polgeotransform[3]
polxmax=polxmin+polcols*polgeotransform[1]
polymin=polymax+polrows*polgeotransform[5]

#Get Central point
polcenterx=(polxmin+polxmax)/2
polcentery=(polymin+polymax)/2

#Raster convert to array
polbands = pol.RasterCount
polband = pol.GetRasterBand(1)
poldataset= polband.ReadAsArray(0,0,polcols,polrows)
poldataimage=poldataset
poldataimage[poldataimage[:,:]==-340282346638528859811704183484516925440.000]=0
#-------------------------------------------------------------------------
#netherlands dataset
nl = gdal.Open('ISG/NLGEO218_20191011.isg') 

#Get coordinates, cols and rows
nlgeotransform = nl.GetGeoTransform()
nlcols = nl.RasterXSize 
nlrows = nl.RasterYSize 

#Get extent
nlxmin=nlgeotransform[0]
nlymax=nlgeotransform[3]
nlxmax=nlxmin+nlcols*nlgeotransform[1]
nlymin=nlymax+nlrows*nlgeotransform[5]

#Get Central point
nlcenterx=(nlxmin+nlxmax)/2
nlcentery=(nlymin+nlymax)/2

#Raster convert to array
nlbands = nl.RasterCount
nlband = nl.GetRasterBand(1)
nldataset= nlband.ReadAsArray(0,0,nlcols,nlrows)
nldataimage=nldataset
nldataimage[nldataimage[:,:]==-340282346638528859811704183484516925440.000]=0
#-------------------------------------------------------------------------
#Visualization in folium
map= folium.Map(location=[centery, centerx], zoom_start=3)

folium.TileLayer('Stamen Terrain').add_to(map)
folium.TileLayer('cartodbpositron').add_to(map)

#--------------------------------------------------------------------------
#geoid layers visualization
#poland
pol = folium.raster_layers.ImageOverlay(
    image=poldataimage,
    bounds=[[polymin, polxmin], [polymax, polxmax]],
    opacity=1,
    show=False,
    
).add_to(map)

#netherlands
nl = folium.raster_layers.ImageOverlay(
    image=nldataimage,
    bounds=[[nlymin, nlxmin], [nlymax, nlxmax]],
    opacity=1,
    show=False,
    
).add_to(map)

# add name to the layers
pol.layer_name = 'Poland'
nl.layer_name = 'Netherlands'

folium.LayerControl().add_to(map)

#Save to html
map.save('html/plotISG.html')
m = folium.Map()
html_string = m.get_root().render()

Best Answer

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
Related Question