Here is an approach for looping.
Environment
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
You can also use
ESRI Satellitefor satellite images. I used this: