Elevation maps that represent surface and river bathymetry are an essential input for flood modeling in software like HEC-RAS. Even with the latest version of high profile open source GIS software as QGIS the combination of a surface elevation map and a river bottom elevation map is a challenge that requires many turnaround, conversion and manual labour. We have developed a useful script that works with surface and river bottom elevation in a "smart" way and creates a geospatial raster with the topobathymetric elevation by the use of geospatial Python libraries as Shapely and Rasterio. The script also includes some key steps to identify the river body and treat the missing values on the bathymetry map.

Tutorial

Code

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)

<matplotlib.image.AxesImage at 0x262ac72a2a0>

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