How to extract max elevation values for each road links in a shapefile using Python libraries?

Question

I want to extract values from a raster DEM file for each feature (polyline) in a shapefile. For each feature in a road shapefile, I want to extract the max and min elevation within a buffer (say 10 ft) of that line.

I am looking for a function similar to AddSurfaceInformation_3d from ESRI ArcGIS.

I attempted using GeoPandas and rasterio with no success. I am receiving only null values.

import geopandas as gpd
gdf = gpd.read_file(pth_to_road_shapefile)

from rasterstats import zonal_stats
with rasterio.open("path_to_test.tif") as src:
    affine = src.transform
    array = src.read(1)
    df_zonal_stats = pd.DataFrame(zonal_stats(gdf, array, affine=affine))

gdf2 = pd.concat([gdf, df_zonal_stats], axis=1) 

Example DEM and vector shapefile here

Answer 1

Your raster DEM and vector road shapefile have different coordinate systems. They need a common coordinate reference system for the rasterstats package to compute the zonal statistics.

import rasterio
import geopandas as gpd
import rasterstats as rs
src = rasterio.open('example_dem.tif')
gdf = gpd.read_file('example_road.shp')
gdf.crs

{'init': 'epsg:4326'}

src.meta['crs']

CRS.from_wkt('PROJCS["NAD83 / Texas South Central (ftUS)",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",27.8333333333333],PARAMETER["central_meridian",-99],PARAMETER["standard_parallel_1",30.2833333333333],PARAMETER["standard_parallel_2",28.3833333333333],PARAMETER["false_easting",1968500],PARAMETER["false_northing",13123333.333],UNIT["US survey foot",0.304800609601219,AUTHORITY["9003","EPSG"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]')

As you can see these two CRS are not the same. To rectify, we can reproject the road shapefile to the DEM CRS.

# Project it
proj_gdf = gdf.to_crs(src.meta['crs'])
# Check the new CRS
proj_gdf.crs

RS.from_wkt('PROJCS["NAD83 / Texas South Central (ftUS)",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",27.8333333333333],PARAMETER["central_meridian",-99],PARAMETER["standard_parallel_1",30.2833333333333],PARAMETER["standard_parallel_2",28.3833333333333],PARAMETER["false_easting",1968500],PARAMETER["false_northing",13123333.333],UNIT["US survey foot",0.304800609601219,AUTHORITY["9003","EPSG"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]')

Now you can see it is the same! We are now ready to perform the zonal statistics:

stats=['min', 'max']
arr = src.read(1)
affine = src.transform
zonal_stats = (rs.zonal_stats(proj_gdf, arr, stats=stats, affine=affine))

Now inspect zonal_stats:

[{'min': 0.1312335878610611, 'max': 40.10498809814453}, {'min': 5.505249500274658, 'max': 30.498685836791992}, {'min': 16.5944881439209, 'max': 38.684383392333984}]