Spatiotemporal Distribution of Dams¶

[1]:

import os
import shutil

import geopandas as gpd
import matplotlib
import matplotlib.pyplot as plt
import pygeohydro as gh
from matplotlib.font_manager import fontManager
from tqdm.auto import tqdm


Get the data from the National Inventory of Dams using PyGeoHydro package within contiguous US.

[2]:

CRS = "esri:102008"
nid = gh.NID()
conus_geom = world[world.iso_a3 == "USA"].geometry.iloc[0].geoms[0]

dam_list = nid.get_byfilter([{"nation": ["USA"]}])
conus_dams = dam_list[0][dam_list[0].is_valid]
conus_dams = conus_dams[conus_dams.within(conus_geom)]
conus_dams = nid.inventory_byid(conus_dams.id)


Since we’re going to use dam coordinates and completion year to filter the data, let’s check number of missing data to find the total available dams that include these two data.

[3]:

missing = ~(conus_dams.is_valid).sum() + conus_dams.yearCompleted.isna().sum()
conus_dams.shape[0] - missing

[3]:

165439


Let’s plot the number of missing data by state.

[4]:

ax = (
conus_dams[~conus_dams.is_valid | conus_dams.yearCompleted.isna()]
.groupby("state")
.size()
.plot.bar(figsize=(20, 5))
)
ax.set_title("Number of dams with missing coordinate or completion year data")
for p in ax.patches:
ax.annotate(
p.get_height(),
(p.get_x() + p.get_width() / 2, p.get_height() + 15),
ha="center",
va="center",
)
ax.figure.set_dpi(300)


Plot the frames.

[5]:

column = "damHeight"
cmap = "viridis"
min_q, max_q = 0.1, 0.9

label = nid.fields_meta[nid.fields_meta["name"] == column].label.iloc[0]
label = "\n".join([label, f"{min_q} - {max_q} Quantile"])
norm = plt.Normalize(
vmin=conus_dams[column].quantile(min_q), vmax=conus_dams[column].quantile(max_q)
)

dpi = 300.0
figsize = (1920.0 / dpi, 1080.0 / dpi)
font = "Lato"
indent = "\n        "
if font in {f.name for f in fontManager.ttflist}:
matplotlib.rcParams["font.sans-serif"] = font
matplotlib.rcParams["font.family"] = "sans-serif"
plt.ioff()

os.makedirs("tmp", exist_ok=True)

conus = gpd.GeoSeries([conus_geom], crs=world.crs).to_crs(CRS)
conus_dams = conus_dams.to_crs(CRS)

def get_ax():
ax = conus.plot(figsize=figsize, facecolor="none", edgecolor="k")
ax.axis(False)
fig = ax.figure
fig.set_dpi(dpi)
[
ax.get_position().x1 + 0.01,
ax.get_position().y0,
0.02,
ax.get_position().height,
]
)
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
fig.colorbar(sm, cax=cax, label=label)
return ax

yr_min = conus_dams.yearCompleted.astype("Int32").min()
yr_max = conus_dams.yearCompleted.astype("Int32").max()
years = range(yr_min + 1, yr_max + 1)
with tqdm(total=len(years), desc="Plotting") as pbar:
for year in years:
pbar.set_postfix_str(f"Year: {year}")
dams = conus_dams[conus_dams.yearCompleted <= year]
ax = get_ax()
dams.plot(
ax=ax,
column=column,
cmap=cmap,
norm=norm,
alpha=0.3,
markersize=3,
)
ax.set_title(f"Dams Completed Up to {year}\nTotal = {len(dams):,}")
h_max = dams[column].max()
name_max = dams.iloc[dams[column].argmax()]["name"].title()
ax.annotate(
f"Largest Dam:{indent}Height: {h_max:.1f} ft{indent}Name: {name_max}",
xy=(0, 0),
xycoords=ax.transAxes,
)
ax.figure.savefig(f"tmp/{year}.png", bbox_inches="tight", dpi=dpi, facecolor="w")
plt.close("all")
pbar.update(1)

Plotting: 100%|██████████| 2021/2021 [10:38<00:00,  3.17it/s, Year: 2021]


Repeat the last frame 100 times.

[6]:

for i in range(1, 100):
shutil.copy(f"tmp/{years[-1]}.png", f"tmp/{years[-1] + i}.png")


Convert the frames to a video file.

[7]:

!ffmpeg -hide_banner -loglevel panic -start_number 1641 -i tmp/%04d.png -pix_fmt yuv420p -vf scale=1920:-2 -y input_data/NID_2019.mp4


You can check out the video here:

[8]:

YouTubeVideo("AM2f9MaBjiQ")

[8]: