3DEP Elevation Grid-to-line intersection

Introduction

In this example, a set of target NHD stream-segments is developed using the HyRiver package. A USGS 3DEP Elevation dataset is pulled using the ClimateR-Catalog. Finally gdptools InterpGen class is used to interpolate both the stream-segments and then interpolate elevatoin from the 3DEP dataset onto those points.

Create conda environment

• gdptools-examples can be downloaded from https://code.usgs.gov/wma/nhgf/toolsteam/gdptools/-/raw/develop/environment-examples.yml?inline=false

Imports

import warnings
warnings.filterwarnings("ignore")
import geopandas as gpd
import pandas as pd
import hvplot.pandas
import hvplot.xarray

import xarray as xr
import pynhd
from pynhd import NLDI, NHDPlusHR, WaterData
from gdptools import UserTiffData
from gdptools import InterpGen
import geoviews as gv

gv.extension('bokeh', 'matplotlib')
# gv.extension()
# gv.output(dpi=120, fig='png')
# gv.output(backend='bokeh')

Create target lines

nldi = NLDI()
station_id = "06730200"

basin = nldi.get_basins(station_id)

flw_main = nldi.navigate_byid(
    fsource="nwissite",
    fid=f"USGS-{station_id}",
    navigation="upstreamMain",
    source="flowlines",
    distance=1000,
)

nhdp_mr = WaterData("nhdflowline_network")
comids = [int(c) for c in flw_main.nhdplus_comid.to_list()]
nhd_mr_fl = nhdp_mr.byid("comid", comids)
nhd_mr_fl.columns.to_list()
flw = pynhd.prepare_nhdplus(nhd_mr_fl, 0, 0, purge_non_dendritic=True)
flw2 = flw.rename(columns={"comid": "ID", "tocomid": "toID"})
strm_segs = pynhd.topoogical_sort(flw2)

# copy the sorted comids into a list and remove the last element
# List will be used to sort ouput into topological order later
sort_ids = strm_segs[0]

# create a new dataframe with the sorted comids
mr_fl = nhd_mr_fl.set_index("comid").loc[sort_ids[:-1]]
mr_fl = mr_fl.reset_index()

mr_fl.head()

	comid	geometry	fdate	resolution	gnis_id	gnis_name	lengthkm	reachcode	flowdir	wbareacomi	...	qc_12	vc_12	qe_12	ve_12	lakefract	surfarea	rareahload	rpuid	vpuid	enabled
0	2889852	MULTILINESTRING Z ((-105.63211 40.05203 0.0000...	2000-12-22T05:00:00Z	Medium	178467	North Boulder Creek	0.373	10190005000122	With Digitized	2888392	...	0.376	-9998.00000	0.376	-9998.00000	1.0	50667.042715	10.004544	10c	10L	1
1	2889186	MULTILINESTRING Z ((-105.62824 40.05363 0.0000...	1999-05-28T04:00:00Z	Medium	178467	North Boulder Creek	0.552	10190005000121	With Digitized	0	...	0.498	0.57306	0.498	0.57306	NaN	NaN	NaN	10c	10L	1
2	2889770	MULTILINESTRING Z ((-105.62207 40.05389 0.0000...	1999-05-28T04:00:00Z	Medium	178467	North Boulder Creek	0.478	10190005000120	With Digitized	2888388	...	0.592	-9998.00000	0.592	-9998.00000	1.0	52090.136727	6.691098	10c	10L	1
3	2889188	MULTILINESTRING Z ((-105.61719 40.05555 0.0000...	1999-05-28T04:00:00Z	Medium	178467	North Boulder Creek	0.396	10190005000119	With Digitized	0	...	0.659	0.62836	0.659	0.62836	NaN	NaN	NaN	10c	10L	1
4	2889772	MULTILINESTRING Z ((-105.61573 40.05258 0.0000...	1999-05-28T04:00:00Z	Medium	178467	North Boulder Creek	0.257	10190005000118	With Digitized	2888396	...	0.785	-9998.00000	0.785	-9998.00000	1.0	48345.911611	4.708286	10c	10L	1

5 rows × 138 columns

Plot the target lines

from bokeh.models import PrintfTickFormatter
formatter = PrintfTickFormatter(format="%i")
gv.tile_sources.EsriTerrain * gv.Path(mr_fl.to_crs(4326), vdims="comid").opts(color='comid',cmap='category20', width=600, height=300, aspect="equal", line_width=5, show_legend=True, colorbar=True, colorbar_opts={'formatter': formatter}, clabel="segment ID", title="NHDPlus MR Flowlines upstream of USGS Gage 06730200")

Define source data

import pprint
climater_cat = "https://mikejohnson51.github.io/climateR-catalogs/catalog.parquet"
cat = pd.read_parquet(climater_cat)

_id = "USGS 3DEP"
_asset = "10m CONUS DEM"
_vars = ["elevation"]
# create list of catalog parameters for each variable
cat_params = [
    cat.query("id == @_id & asset == @_asset & variable == @_vars").to_dict(orient="records")[0]
    for _var in _vars
]
# create dictionary of variable names and catalog parameters
cat_dict = dict(zip(_vars, cat_params))
pprint.pprint(cat_dict)

{'elevation': {'T_name': None,
               'URL': '/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/USGS_Seamless_DEM_13.vrt',
               'X1': -180.0005555560936,
               'X_name': None,
               'Xn': 180.00000026239937,
               'Y1': -180.0005555560936,
               'Y_name': None,
               'Yn': 72.0005555562949,
               'asset': '10m CONUS DEM',
               'crs': '+proj=longlat +datum=NAD83 +no_defs',
               'description': '10m CONUS DEM',
               'duration': None,
               'ensemble': None,
               'id': 'USGS 3DEP',
               'interval': None,
               'model': None,
               'nT': nan,
               'ncols': 3888006.0,
               'nrows': 939612.0,
               'resX': 9.259259266022042e-05,
               'resY': 9.25925926599094e-05,
               'scenario': None,
               'tiled': '',
               'toptobottom': None,
               'type': 'vrt',
               'units': 'meters',
               'variable': 'elevation',
               'varname': 'elevation'}}

Load source data

import rioxarray as rxr
ds_3dep = rxr.open_rasterio(cat_dict["elevation"]["URL"])
ds_3dep

<xarray.DataArray (band: 1, y: 939612, x: 3888006)>
[3653217093672 values with dtype=float32]
Coordinates:
  * band         (band) int64 1
  * x            (x) float64 -180.0 -180.0 -180.0 -180.0 ... 180.0 180.0 180.0
  * y            (y) float64 72.0 72.0 72.0 72.0 ... -15.0 -15.0 -15.0 -15.0
    spatial_ref  int64 0
Attributes:
    _FillValue:    -999999.0
    scale_factor:  1.0
    add_offset:    0.0

xarray.DataArray

• band: 1
• y: 939612
• x: 3888006

• ...

  [3653217093672 values with dtype=float32]

• Coordinates: (4)
  • band
    (band)
    int64
    1

    array([1])

  • x
    (x)
    float64
    -180.0 -180.0 ... 180.0 180.0

    array([-180.000509, -180.000417, -180.000324, ...,  179.999769,  179.999861,
            179.999954])

  • y
    (y)
    float64
    72.0 72.0 72.0 ... -15.0 -15.0

    array([ 72.000509,  72.000417,  72.000324, ..., -15.000324, -15.000417,
           -15.000509])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101004,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4269"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314140356

    inverse_flattening :

    298.257222101004

    reference_ellipsoid_name :

    GRS 1980

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    NAD83

    horizontal_datum_name :

    North American Datum 1983

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101004,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4269"]]

    GeoTransform :

    -180.0005555560936 9.259259266022042e-05 0.0 72.0005555562949 0.0 -9.25925926599094e-05

    array(0)

• Indexes: (3)
  • band
    PandasIndex

    PandasIndex(Index([1], dtype='int64', name='band'))

  • x
    PandasIndex

    PandasIndex(Index([-180.00050925979727,  -180.0004166672046, -180.00032407461194,
            -180.0002314820193, -180.00013888942664, -180.00004629683397,
            -179.9999537042413, -179.99986111164864, -179.99976851905598,
           -179.99967592646334,
           ...
             179.9991206327691,  179.99921322536173,  179.99930581795437,
            179.99939841054706,   179.9994910031397,   179.9995835957324,
            179.99967618832503,  179.99976878091772,  179.99986137351036,
              179.999953966103],
          dtype='float64', name='x', length=3888006))

  • y
    PandasIndex

    PandasIndex(Index([  72.00050925999857,    72.0004166674059,   72.00032407481325,
             72.00023148222058,   72.00013888962792,   72.00004629703527,
              71.9999537044426,   71.99986111184995,   71.99976851925729,
             71.99967592666462,
           ...
            -14.99967598843763, -14.999768581030281, -14.999861173622946,
           -14.999953766215611, -15.000046358808262, -15.000138951400928,
           -15.000231543993578, -15.000324136586244, -15.000416729178909,
            -15.00050932177156],
          dtype='float64', name='y', length=939612))

• Attributes: (3)

  _FillValue :

  -999999.0

  scale_factor :

  1.0

  add_offset :

  0.0

Define UserTiffData object

varname = "elevation"
tx_name = "x"
ty_name = "y"
band = 1
bandname = "band"
crs = 4269
f_crs = 4326
id_feature = "comid"
user_data = UserTiffData(
    var=varname,
    ds=ds_3dep,
    proj_ds=crs,
    x_coord=tx_name,
    y_coord=ty_name,
    band=band,
    bname=bandname,
    f_feature=mr_fl,
    id_feature=id_feature,
    proj_feature=crs
)

Define InterpGen object and generate interpolation

interpObject1 = InterpGen(user_data, pt_spacing=100, stat='all', interp_method="linear")
# stats1, pts1 = interpObject1.calc_interp()

stats1, pts1 = interpObject1.calc_interp()

Plot the sub-setted 3DEP 10m DEM

ds_ss = user_data.get_source_subset(key="").to_dataset(name="elevation")
ds_ss.elevation.hvplot.quadmesh(x="x", y="y", geo=True, project=True, rasterize=True,  cmap="viridis", width=600, height=300, aspect="equal", clabel="3DEP 10m DEM", xlabel="Longitude", ylabel="Latitude", title="3DEP 10m DEM")

Show the interpolated points

pts1.head(10)

	band	spatial_ref	values	dist	comid	geometry	varname
pt
0	1	0	3620.148240	0.0	2889186	POINT (-105.62824 40.05363)	elevation
1	1	0	3620.536627	100.0	2889186	POINT (-105.62721 40.05347)	elevation
2	1	0	3605.426680	200.0	2889186	POINT (-105.62617 40.05330)	elevation
3	1	0	3589.268640	300.0	2889186	POINT (-105.62513 40.05317)	elevation
4	1	0	3578.737563	400.0	2889186	POINT (-105.62408 40.05319)	elevation
5	1	0	3571.558674	500.0	2889186	POINT (-105.62307 40.05350)	elevation
6	1	0	3561.002174	600.0	2889186	POINT (-105.62208 40.05386)	elevation
0	1	0	3562.425944	0.0	2889188	POINT (-105.61719 40.05555)	elevation
1	1	0	3560.132308	100.0	2889188	POINT (-105.61623 40.05523)	elevation
2	1	0	3524.454447	200.0	2889188	POINT (-105.61593 40.05429)	elevation

Show The statistics

stats1.head(10)

	comid	varname	mean	median	std	min	max
0	2889186	elevation	3592.382657	3589.268640	21.869675	3561.002174	3620.536627
0	2889188	elevation	3530.434489	3542.293378	35.513002	3474.725256	3562.425944
0	2889198	elevation	3448.163403	3448.163403	5.088806	3443.074596	3453.252209
0	2889214	elevation	1558.948776	1558.909865	1.177185	1556.995965	1560.632596
0	2889222	elevation	3383.079127	3384.111898	20.186191	3357.855991	3407.269491
0	2889228	elevation	1562.467624	1562.329515	1.200751	1560.514595	1564.282919
0	2889258	elevation	1566.994821	1567.276414	1.405481	1564.551665	1569.605381
0	2889260	elevation	1570.061471	1569.926945	1.177150	1568.664712	1572.449691
0	2889264	elevation	1572.780278	1572.780278	0.221672	1572.558607	1573.001950
0	2889276	elevation	1574.920845	1574.883089	1.282822	1572.987986	1576.906272

Plot the mean elevation by comid

# sort stats by topological stream order
stats = stats1.set_index("comid").loc[sort_ids[:-1]]
mr_fl["tmean"] = stats["mean"].values
pmin = mr_fl["tmean"].min()
pmax = mr_fl["tmean"].max()
pmin, pmax
import geoviews as gv

gv.extension("bokeh")
gv.tile_sources.EsriTerrain * gv.Path(mr_fl.to_crs(4326), vdims=["tmean"]).opts(width=600, height=300, aspect="equal",color="tmean", cmap="viridis", line_width=5, colorbar=True, clim=(pmin, pmax), clabel="Mean Elevation (m)", title="Mean Elevation by Stream Segment")

Plot the interpolated data as a scatter plot by comid

import plotly.graph_objects as go
fig = go.Figure()
pts = pts1.groupby("comid")
for index in range(len(pts)):
    group = pts.get_group(sort_ids[index])
    fig.add_trace(go.Scatter(x=group["dist"], y=group["values"], mode='markers', marker=dict(size=5, color=group["values"], colorscale='Viridis', cmin=pmin, cmax=pmax), name=str(sort_ids[index])))
fig.update_layout(title="Stream-segment Interpolated Elevation 10m DEM", xaxis_title="Distance (m)", yaxis_title="Elevation (m)")
fig.show()

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[13], line 1
----> 1 import plotly.graph_objects as go
      2 fig = go.Figure()
      3 pts = pts1.groupby("comid")

ModuleNotFoundError: No module named 'plotly'