Ctrl+K

This page was generated from us_coasts.ipynb. Interactive online version: Binder badge

CONUS Coasts#

US States and Coastlines#

Census TIGER (Topologically Integrated Geographic Encoding and Referencing database)

[1]:

import warnings
from pathlib import Path

import geopandas as gpd

warnings.filterwarnings("ignore", message=".*initial implementation of Parquet.*")
root = Path("input_data")
root.mkdir(exist_ok=True)
BASE_PLOT = {"facecolor": "k", "edgecolor": "b", "alpha": 0.2, "figsize": (12, 6)}
CRS = "epsg:5070"


def tiger_shp_2022(data):
    url = f"https://www2.census.gov/geo/tiger/TIGER2022/{data.upper()}/tl_2022_us_{data}.zip"
    return gpd.read_file(url)


cfile = Path(root, "us_state.feather")
if cfile.exists():
    state = gpd.read_feather(cfile)
else:
    state = tiger_shp_2022("state")
    state.to_feather(cfile)

[2]:

ax = state.to_crs(CRS).plot(**BASE_PLOT)
ax.axis("off")
ax.margins(0)
../../_images/examples_notebooks_us_coasts_3_0.png 
[3]:

cfile = Path(root, "us_coastline.feather")
if cfile.exists():
    coastline = gpd.read_feather(cfile)
else:
    coastline = tiger_shp_2022("coastline")
    coastline.to_feather(cfile)

[4]:

ax = coastline.to_crs(CRS).plot(figsize=(12, 6))
ax.axis("off")
ax.margins(0)
../../_images/examples_notebooks_us_coasts_5_0.png

Clip to CONUS#

[5]:

from shapely.geometry import box

conus_bounds = box(-125, 24, -65, 50)

cfile = Path(root, "conus_states.feather")
if cfile.exists():
    conus_states = gpd.read_feather(cfile)
else:
    conus_states = state[state.within(conus_bounds)]
    conus_states.to_feather(cfile)

[6]:

ax = conus_states.to_crs(CRS).plot(**BASE_PLOT)
ax.axis("off")
ax.margins(0)
../../_images/examples_notebooks_us_coasts_8_0.png 
[7]:

conus_coastline = coastline[coastline.within(conus_bounds)]

cfile = Path(root, "us_coast_states.feather")
if cfile.exists():
    coast_states = gpd.read_feather(cfile)
else:
    coast_states = state[state.intersects(conus_coastline.unary_union)]
    coast_states.to_feather(cfile)

[8]:

ax = coast_states.to_crs(CRS).plot(**BASE_PLOT)
conus_coastline.to_crs(CRS).plot(ax=ax, edgecolor="b", lw=2, zorder=1)
ax.axis("off")
ax.margins(0)
../../_images/examples_notebooks_us_coasts_10_0.png

Tidal and Estuary USGS stations#

58575595ff9641b980af424fbff5010a

We need to look at the Water Services API.

[9]:

import pandas as pd

from pygeohydro import NWIS, ZeroMatchedError

cfile = Path(root, "coast_stations.feather")
if cfile.exists():
    coast_stations = gpd.read_feather(cfile)
else:
    queries = [
        {
            "stateCd": s.lower(),
            "siteType": "ST-TS,ES",
            "hasDataTypeCd": "dv",
            "outputDataTypeCd": "dv",
        }
        for s in coast_states.STUSPS
    ]
    nwis = NWIS()

    def get_info(q):
        try:
            return nwis.get_info(q, True)
        except ZeroMatchedError:
            return None

    sites = pd.concat(get_info(q) for q in queries if q is not None)

    coast_stations = gpd.GeoDataFrame(
        sites,
        geometry=gpd.points_from_xy(sites.dec_long_va, sites.dec_lat_va),
        crs="epsg:4269",
    )
    coast_stations.to_feather(cfile)

st = coast_stations[["site_no", "site_tp_cd", "geometry"]].to_crs(CRS)
ts = st[st.site_tp_cd == "ST-TS"].drop_duplicates()
es = st[st.site_tp_cd == "ES"].drop_duplicates()

[10]:

ax = conus_states.to_crs(CRS).plot(**BASE_PLOT)
coastline[coastline.within(conus_bounds)].to_crs(CRS).plot(ax=ax, edgecolor="g", lw=2, zorder=1)
ts.plot(ax=ax, lw=3, c="r")
es.plot(ax=ax, lw=3, c="b")
ax.legend(["Coastline", f"ST-TS ({ts.shape[0]})", f"ES ({es.shape[0]})"], loc="best")
ax.axis("off")
ax.margins(0)
ax.figure.savefig(Path("_static", "us_coasts.png"), dpi=300, bbox_inches="tight", facecolor="w")
../../_images/examples_notebooks_us_coasts_14_0.png

Mean daily discharge for all stations#

[11]:

import numpy as np
import pandas as pd

cfile = Path(root, "discharge.parquet")
dates = ("2000-01-01", "2015-12-31")

if cfile.exists():
    discharge = pd.read_parquet(cfile)
else:
    nwis = NWIS()
    discharge = nwis.get_streamflow(
        coast_stations.site_no,
        dates,
    )
    discharge[discharge < 0] = np.nan
    discharge.to_parquet(cfile)

[12]:

ax = discharge.plot(legend=False, lw=0.8, figsize=(7, 3.5))
ax.set_ylabel("Q (cms)")
ax.set_xlabel("")
ax.margins(x=0)
../../_images/examples_notebooks_us_coasts_17_0.png

Let’s find the station that has the largest peak value and see it on a map.

[13]:

station_id = discharge.max().idxmax().split("-")[1]

coast_stations[coast_stations.site_no == station_id].iloc[0]["station_nm"]

[13]:

'COLUMBIA RIVER AT PORT WESTWARD, NEAR QUINCY, OR'

[14]:

import pygeohydro as gh

lat, lon = coast_stations[coast_stations.site_no == station_id].iloc[0][
    ["dec_lat_va", "dec_long_va"]
]
bbox = (lon - 0.2, lat - 0.2, lon + 0.2, lat + 0.2)
nwis_kwds = {"hasDataTypeCd": "dv", "outputDataTypeCd": "dv", "parameterCd": "00060"}

station_map = gh.interactive_map(bbox, nwis_kwds=nwis_kwds)

[15]:

station_map

[15]:
Make this Notebook Trusted to load map: File -> Trust Notebook

River network data#

75d7cf0bdd8247aea37e452d93d251b0

Basin#

[16]:

import pynhd as nhd

nldi = nhd.NLDI()
cfile = Path(root, "basin.feather")
if cfile.exists():
    basin = gpd.read_feather(cfile)
else:
    basin = nldi.get_basins(station_id)
    basin.to_feather(cfile)

[17]:

ax = basin.to_crs(CRS).plot(**BASE_PLOT)
ax.axis("off")
ax.margins(0)
../../_images/examples_notebooks_us_coasts_25_0.png 
[18]:

import folium

folium.GeoJson(basin.geometry).add_to(station_map)
station_map

[18]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Main river network#

[19]:

cfile = Path(root, "flowline_main.feather")
if cfile.exists():
    flw_main = gpd.read_feather(cfile)
else:
    flw_main = nldi.navigate_byid(
        fsource="nwissite",
        fid=f"USGS-{station_id}",
        navigation="upstreamMain",
        source="flowlines",
        distance=4000,
    )
    flw_main.to_feather(cfile)

Tributaries#

[20]:

cfile = Path(root, "flowline_trib.feather")
if cfile.exists():
    flw_trib = gpd.read_feather(cfile)
else:
    flw_trib = nldi.navigate_byid(
        fsource="nwissite",
        fid=f"USGS-{station_id}",
        navigation="upstreamTributaries",
        source="flowlines",
        distance=4000,
    )
    flw_trib.to_feather(cfile)
flw_trib["nhdplus_comid"] = flw_trib["nhdplus_comid"].astype("float").astype("Int64")

[21]:

BASE_PLOT["figsize"] = (8, 5)
ax = basin.plot(**BASE_PLOT)
flw_trib.plot(ax=ax)
flw_main.plot(ax=ax, lw=3, color="r")
ax.legend(["Tributaries", "Main"])
ax.axis("off")
ax.margins(0)
../../_images/examples_notebooks_us_coasts_31_0.png

NHDPlus Value Added Attributes (VAA)#

[22]:

cfile = Path(root, "nhdplus_vaa.parquet")
nhdplus_vaa = nhd.nhdplus_vaa(cfile)

flw_trib = flw_trib.merge(nhdplus_vaa, left_on="nhdplus_comid", right_on="comid", how="left")
flw_trib.loc[flw_trib.slope < 0, "slope"] = np.nan

[23]:

ax = basin.plot(**BASE_PLOT)
flw_trib.plot(
    ax=ax,
    column="slope",
    scheme="natural_breaks",
    k=6,
    cmap="plasma",
    legend=True,
    legend_kwds={"title": "Slope [-]"},
)
ax.axis("off")
ax.margins(0)
../../_images/examples_notebooks_us_coasts_34_0.png
On this page
Edit on GitHub
Show Source