USGS 3DEP Lidar elevation data with xarray, rasterio, and matplotlib

lidar3dep.py

import matplotlib.pyplot as plt
import xarray as xr
import glob

# The code below is for 4 DEMs that cover the Montford Dam that forms Lake Alan Henry.
# Download geotiffs here. 
: https://viewer.nationalmap.gov/basic/?basemap=b1&category=ned,nedsrc&q=&zoom=4&bbox=-139.74609375,10.14193169,-54.22851563,61.14323525&preview=&avail=&refpoly=

fns = glob.glob('/Users/ebruning/Downloads/*.tif')
# ds=xr.open_rasterio('/Users/ebruning/Downloads/USGS_one_meter_x31y366_TX_West_Central_B12_2018.tif')
all_ds = [xr.open_rasterio(f)[0,:,:].to_dataset(name='dep') for f in fns]

fig = plt.figure(figsize=(17,17), dpi=300)
ax = fig.add_subplot(111)

for ds in all_ds:
    xsel = (ds.x<.314e6) & (ds.x>.306e6)
    ysel = (ds.y>3.657e6) & (ds.y<3.663e6)
    img = ds.dep[ysel,xsel]
#         NGVD 29 = NAVD 88 –3.6 feet. Emergency spillway crest is 2240, operational spillway 2220, top of dam is 2263 in NGVD29
    img = ax.pcolormesh(img.x,img.y, (img.data*3.28084 - 3.6), vmin=2100, vmax=2400)

ax.set_title('Lake Alan Henry, Montford Dam, 2018, 3DEP lidar elevation (feet)')
ax.set_xlabel('UTM east (meters)')
ax.set_ylabel('UTM north (meters)')
fig.tight_layout()
plt.colorbar(img)