This page was generated from notebooks/aridity.ipynb. Interactive online version: Binder badge

Aridity

PyDaymet computes Potential Evaportranspiration based on FAO Penman-Monteith equation using tmin (deg c), tmax (deg c), vp (Pa), srad (W/m^2), and dayl (s) variables from Daymet:

image.png

Basin

Let’s start by getting the basin geometry of USGS-01031500 station (Piscataquis River near Dover-Foxcroft, Maine)

[1]:
import os
import warnings
from pathlib import Path

import geopandas as gpd
import pynhd as nhd

warnings.filterwarnings("ignore", message=".*initial implementation of Parquet.*")
root = Path("input_data")
os.makedirs(root, exist_ok=True)
BASE_PLOT = {"facecolor": "k", "edgecolor": "b", "alpha": 0.2, "figsize": (18, 9)}
CRS = "esri:102008"
station_id = "01031500"

nldi = nhd.NLDI()
cfile = Path(root, f"basin_{station_id}.feather")
if cfile.exists():
    basin = gpd.read_feather(cfile)
else:
    basin = nldi.get_basins(station_id)
    basin.to_feather(cfile)
[2]:
ax = basin.to_crs(CRS).plot(**BASE_PLOT)
ax.axis("off")
ax.margins(0)
../_images/notebooks_aridity_4_0.png

Precipitation and Potential Evapotranspiration

[3]:
import pydaymet as daymet
from tqdm.notebook import tqdm

years = list(range(2006, 2016))
geometry = basin.iloc[0].geometry

for yr in tqdm(years, desc="Getting Climate"):
    cfile = Path(root, f"clm_{yr}.nc")
    if cfile.exists():
        continue
    daymet.get_bygeom(geometry, yr, variables="prcp", pet="hargreaves_samani").to_netcdf(cfile)
[4]:
import xarray as xr

clm = xr.open_mfdataset(Path(root).glob("clm_*.nc"), coords="minimal")
clm["elevation"] = clm["elevation"].isel(time=0, drop=True)
clm
[4]:
<xarray.Dataset>
Dimensions:    (time: 3650, y: 37, x: 40)
Coordinates:
  * time       (time) datetime64[ns] 2006-01-01T12:00:00 ... 2015-12-31T12:00:00
  * y          (y) float32 719.0 718.0 717.0 716.0 ... 686.0 685.0 684.0 683.0
  * x          (x) float32 2.21e+03 2.211e+03 2.212e+03 ... 2.248e+03 2.249e+03
Data variables:
    dayl       (time, y, x) float32 dask.array<chunksize=(365, 37, 40), meta=np.ndarray>
    lat        (time, y, x) float64 dask.array<chunksize=(365, 37, 40), meta=np.ndarray>
    lon        (time, y, x) float64 dask.array<chunksize=(365, 37, 40), meta=np.ndarray>
    prcp       (time, y, x) float32 dask.array<chunksize=(365, 37, 40), meta=np.ndarray>
    srad       (time, y, x) float32 dask.array<chunksize=(365, 37, 40), meta=np.ndarray>
    tmax       (time, y, x) float32 dask.array<chunksize=(365, 37, 40), meta=np.ndarray>
    tmin       (time, y, x) float32 dask.array<chunksize=(365, 37, 40), meta=np.ndarray>
    vp         (time, y, x) float32 dask.array<chunksize=(365, 37, 40), meta=np.ndarray>
    elevation  (y, x) float64 dask.array<chunksize=(37, 40), meta=np.ndarray>
    pet        (time, y, x) float64 dask.array<chunksize=(365, 37, 40), meta=np.ndarray>
Attributes: (12/17)
    start_year:          2006
    source:              Daymet Software Version 4.0
    Version_software:    Daymet Software Version 4.0
    Version_data:        Daymet Data Version 4.0
    Conventions:         CF-1.6
    citation:            Please see http://daymet.ornl.gov/ for current Dayme...
    ...                  ...
    geospatial_lon_max:  -69.13309404831536
    crs:                 +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-10...
    nodatavals:          0.0
    transform:           [ 1.00000e+00  0.00000e+00  2.20625e+03  0.00000e+00...
    res:                 [1. 1.]
    bounds:              [2209.24651856  682.07851575 2250.2123966   719.9719...
[5]:
clm.elevation.plot()
[5]:
<matplotlib.collections.QuadMesh at 0x1cd0fae30>
../_images/notebooks_aridity_8_1.png
[6]:
import hvplot.xarray  # noqa

clm.hvplot.violin(y="pet", by="time.month")