1

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()

2
  • Can you clarify what problems have occurred using your code? Aug 14, 2021 at 12:24
  • hi, the code works fine. I just want to make it more efficient to handle multiple rasters. By looking at several readings, using 'for loop' might work, but currently unable to do so (I'm pretty beginner) and I do not find (yet) specific source of how to import multiple rasters with efficient way. That is why I tried to ask here
    – duns
    Aug 14, 2021 at 13:06

1 Answer 1

1

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
1
  • Thank you for the answer! This is what I need. I use Spyder, but I could try tinker with your solution a bit to make it work for me.
    – duns
    Aug 15, 2021 at 14:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.