Generación de elevaciones topobatimétricas para modelamiento de inundaciones con Python geoespacial

Los mapas de elevación que representan la superficie y la batimetría de los ríos son una entrada esencial para el modelamiento de inundaciones en software como HEC-RAS. Incluso con las últimas versiones de software GIS de código abierto de alto perfil, como QGIS, la combinación de un mapa de elevación de la superficie y un mapa de elevación del fondo del río es un desafío que requiere muchos ajustes, conversiones y trabajo manual. Hemos desarrollado un script útil que trabaja con las elevaciones de la superficie y el fondo del río de manera "inteligente" y crea un ráster geoespacial con la elevación topobatimétrica utilizando bibliotecas de Python geoespacial como Shapely y Rasterio. El script también incluye algunos pasos clave para identificar el cuerpo del río y tratar los valores faltantes en el mapa de batimetría.

Tutorial

Código

import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import rasterio
from rasterio.plot import show
from shapely.geometry import Polygon, Point
from rasterio.features import shapes
from rtree import index

#open the raster utm datasets
riverBat = rasterio.open('../Rst/riverBathymetry.tif')
surfElev = rasterio.open('../Rst/surfaceElevation_15N.tif')
show(riverBat)

#get no data value
rivNodata = riverBat.nodata
rivNodata

-3.4028230607370965e+38

#create an array of the active river
bath = riverBat.read(1)
activeRiver = np.zeros(bath.shape)
activeRiver[bath > riverBat.nodata] = 1
plt.imshow(activeRiver)

#creating a contour over the active river
polys = shapes(activeRiver, connectivity=8, transform=riverBat.transform)

#looping over the polygon generator to create a geospatial polygon
hugePolys = []

for poly in polys:
    #print(poly)
    shpPol = Polygon(poly[0]['coordinates'][0])
    if shpPol.area > 100000 and poly[1]==1:
        #print(poly[0])
        print(shpPol.area)
        hugePolys.append(shpPol)

1187999.2424999955

# show resulting polygon
hugePolys[0]

svg

#open the aoi
aoi = gpd.read_file('../Shp/aoi_15N.shp')
aoi.iloc[0].geometry

svg

#explore the aoi bounds
bounds = aoi.total_bounds
bounds

array([ 430657.62111644, 4338698.67453457,  433523.37037566,
       4340610.00927973])

xMin, xMax = bounds[0], bounds[2]
yMin, yMax = bounds[1], bounds[3]

res = 25 #output raster resolution of 1m

#define the cell centers arrays
xArray = np.arange(xMin,xMax,res)
yArray = np.arange(yMin,yMax,res)
yArray[:5]

array([4338698.67453457, 4338723.67453457, 4338748.67453457,
       4338773.67453457, 4338798.67453457])

#get nrows and ncols
nrow = yArray.shape[0]
ncol = xArray.shape[0]
print(nrow,ncol)

77 115

#get no data value
riverNodata = riverBat.nodata

#create sample array
topoBath = np.zeros([nrow,ncol])

i = 0

for row in range(nrow):
    for col in range(ncol):
        
        colX = xArray[col]
        rowY = yArray[row]

        # create a point object of the xy
        pixel = Point(colX,rowY)
        #apply value of bathymetry if available
        if hugePolys[0].contains(pixel):
            #extract raster values
            bathValue = list(riverBat.sample([(colX,rowY)]))[0]
            bathValue = bathValue.tolist()[0]
            #check if sampled values are empty
            if bathValue == riverNodata:
                topoBath[row,col] = topoBath[row,col-1]
            else:
                topoBath[row,col] = bathValue
        else:
            surfValue = list(surfElev.sample([(colX,rowY)]))[0]
            surfValue = surfValue.tolist()[0]
            topoBath[row,col] = surfValue

        #counter
        if i % 100 ==0:
            print(f'\r %d'%i, end="", flush=True)

        i+=1

#preview of the topo array (inverted)
plt.imshow(topoBath)

#definition of the raster transform array
from rasterio.transform import Affine
transform = Affine.translation(xArray[0]-res/2, yArray[0]-res/2)*Affine.scale(res,res)
transform

Affine(25.0, 0.0, 430645.12111643935,
       0.0, 25.0, 4338686.1745345695)

#get crs as wkt
from rasterio.crs import CRS
rasterCrs = CRS.from_epsg(26915)
rasterCrs.data

{'proj': 'utm', 'zone': 15, 'datum': 'NAD83', 'units': 'm', 'no_defs': True}

#definition, register and close of interpolated raster
interpRaster = rasterio.open('../Rst/topoBath_'+str(res)+'.tif',
                                'w',
                                driver='GTiff',
                                height=nrow,
                                width=ncol,
                                count=1,
                                dtype=topoBath.dtype,
                                crs=rasterCrs,
                                transform=transform,
                                )
interpRaster.write(topoBath,1)
interpRaster.close()