Elevation Differencing

Contents

Elevation Differencing¶

Canopy Height and Snow Depth¶

Canopy Height Model is a popular raster product from lidar data. It can be derived from differencing DSM and DEM. Snow depth is another raster product of interest. This is derived by differencing snow-on DEM and snow-off DEM.

We will difference the appropriate elevation model to derive vegetation height and snow depth. To begin, let’s load the necessary packages for this tutorial.

#import packages

import os

import numpy as np 

#plotting packages
import matplotlib.pyplot as plt
import seaborn as sns
import hvplot.xarray

#geospatial packages
import geopandas as gpd #for vector data
import xarray as xr
import rioxarray #for raster data
import rasterstats as rs

More Analysis¶

Lidar provides highly accurate DEM, vegetation height and snow depth. Our goal here is to quickly demonstrate some of the analysis we can do with these products. We will sample snow depth, elevation and vegetation height at 100 locations to quantify snow distribution over the the subset region

#import point location data
cam_100_points = gpd.read_file('/tmp/input/cameron_100_points.gpkg')

#create a dir to store the output
if not os.path.exists("./output"):
    os.makedirs("./output")

#create a buffer around the point feature and save to file
cam_100_points.buffer(20).to_file('./output/cam_100_points_buffer.gpkg')

Let’s see the sampled points

cam_100_poly = gpd.read_file('./output/cam_100_points_buffer.gpkg')
fig, ax = plt.subplots(figsize=(10,12))
cam_sd_roi.plot(ax = ax)
cam_100_poly.plot(ax = ax, color = 'black')
plt.show()

../../_images/2_elevation_differencing_17_0.png

#Extract raster values for the vegetation height
cam_VH_stat = rs.zonal_stats('./output/cam_100_points_buffer.gpkg', 
            cam_vh_roi.squeeze().values,
            nodata = -9999,
            geojson_out=True,
            affine = cam_vh_roi.rio.transform(), 
            stats = 'count mean')
cam_VH_stat_df =  gpd.GeoDataFrame.from_features(cam_VH_stat) #turn the geojson to a dataframe
cam_VH_stat_df.rename(columns = {'mean': 'VH_mean', 'count': 'VH_count'}, inplace = True) #rename the columns

#Extract raster values for the snow depth
cam_SD_stat = rs.zonal_stats(cam_VH_stat_df,
            cam_sd_roi.squeeze().values,
            nodata = -9999,
            geojson_out=True,
            affine = cam_sd_roi.rio.transform(),
            stats = 'count mean')
cam_SD_stat_df =  gpd.GeoDataFrame.from_features(cam_SD_stat) #turn the geojson to a dataframe
cam_SD_stat_df.rename(columns = {'mean': 'SD_mean', 'count': 'SD_count'}, inplace = True) #rename the columns

#Extract raster values for elevation
cam_DEM_stat = rs.zonal_stats(cam_SD_stat_df,
            cam_dem_roi.squeeze().values,
            nodata = -9999,
            geojson_out=True,
            affine = cam_dem_roi.rio.transform(),
            stats = 'count mean')
cam_DEM_stat_df =  gpd.GeoDataFrame.from_features(cam_DEM_stat) #turn the geojson to a dataframe
cam_DEM_stat_df.rename(columns = {'mean': 'ELEV_mean', 'count': 'ELEV_count'}, inplace = True) #rename the columns

#rename the dataframe
vh_sd_dem = cam_DEM_stat_df

# Set font size and font family
plt.rcParams.update({'font.size': 18})
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = ['Times New Roman'] + plt.rcParams['font.serif']

#create figure and gridspec object
fig = plt.figure(figsize=(20,12), constrained_layout=True)
gspec = fig.add_gridspec(ncols=5, nrows=5)

ax0 = fig.add_subplot(gspec[0, 0])
ax1 = fig.add_subplot(gspec[0, 1])
ax2 = fig.add_subplot(gspec[0, 2])
ax3 = fig.add_subplot(gspec[1:5, 0])
ax4 = fig.add_subplot(gspec[1:5, 1])
ax5 = fig.add_subplot(gspec[1:5, 2])
ax6 = fig.add_subplot(gspec[0:5, 3])
ax7 = fig.add_subplot(gspec[0:5, 4])

# plot the boxplot of elevation, vegetation and snow thickness
sns.boxplot(x= vh_sd_dem['ELEV_mean'], ax=ax0, palette= 'colorblind')
ax0.set_xlabel('')
sns.boxplot(x= vh_sd_dem['VH_mean'], ax=ax1, palette= 'colorblind')
ax1.set_xlabel('')
sns.boxplot(x= vh_sd_dem['SD_mean'], ax=ax2, palette= 'colorblind')
ax2.set_xlabel('')

# plot the histogram and boxplot of elevation, vegetation and snow thickness
sns.histplot(data = vh_sd_dem, x ='ELEV_mean', ax=ax3, bins=20)
ax3.set_xlabel('Elevation (m)')
sns.histplot(data = vh_sd_dem, x ='VH_mean', ax=ax4, bins=20, binrange=(1,10))
ax4.set_xlabel('Vegetation Height (m)')
sns.histplot(data = vh_sd_dem, x ='SD_mean', ax=ax5, bins=30, binrange=(0,1.5))
ax5.set_xlabel('Snow Depth (m)')

#make the scatter plot of elevation vs snow thickness and vegetation vs snow thickness
sns.kdeplot(data = vh_sd_dem, x = 'SD_mean', y = 'ELEV_mean', shade=True, ax=ax6, cmap='Blues')
ax6.set_xlabel('Snow Depth (m)')
ax6.set_ylabel('Elevation (m)')
#set xlimit and y limit
#ax6.set_xlim(0, 2.5)
ax6.set_ylim(2900, 3300)

sns.kdeplot(data = vh_sd_dem, x = 'SD_mean', y = 'VH_mean', shade=True, ax=ax7, cmap='Blues')
ax7.set_xlabel('Snow Depth (m)')
ax7.set_ylabel('Vegetation Height (m)')
#set xlimit and y limit
#ax7.set_xlim(0, 2.5)
ax7.set_ylim(-3, 15)

plt.show()

../../_images/2_elevation_differencing_19_0.png

previous

Introduction to Lidar Remote Sensing

next

Common Pitfall in Elevation differencing