This page was generated from aridity.ipynb. Interactive online version:

Aridity

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

from pathlib import Path

import geopandas as gpd

import pynhd as nhd

root = Path("input_data")
root.mkdir(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)

ax = basin.to_crs(CRS).plot(**BASE_PLOT)
ax.axis("off")
ax.margins(0)

Precipitation and Potential Evapotranspiration

from tqdm.notebook import tqdm

import pydaymet as daymet

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)

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

<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...

xarray.Dataset

• Dimensions:
  • time: 3650
  • y: 37
  • x: 40
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2006-01-01T12:00:00 ... 2015-12-...

    standard_name :

    time

    bounds :

    time_bnds

    long_name :

    24-hour day based on local time

    _ChunkSizes :

    1

    _CoordinateAxisType :

    Time

    array(['2006-01-01T12:00:00.000000000', '2006-01-02T12:00:00.000000000',
           '2006-01-03T12:00:00.000000000', ..., '2015-12-29T12:00:00.000000000',
           '2015-12-30T12:00:00.000000000', '2015-12-31T12:00:00.000000000'],
          dtype='datetime64[ns]')

  • y
    (y)
    float32
    719.0 718.0 717.0 ... 684.0 683.0

    units :

    km

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    array([719., 718., 717., 716., 715., 714., 713., 712., 711., 710., 709., 708.,
           707., 706., 705., 704., 703., 702., 701., 700., 699., 698., 697., 696.,
           695., 694., 693., 692., 691., 690., 689., 688., 687., 686., 685., 684.,
           683.], dtype=float32)

  • x
    (x)
    float32
    2.21e+03 2.211e+03 ... 2.249e+03

    units :

    km

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    array([2209.75, 2210.75, 2211.75, 2212.75, 2213.75, 2214.75, 2215.75, 2216.75,
           2217.75, 2218.75, 2219.75, 2220.75, 2221.75, 2222.75, 2223.75, 2224.75,
           2225.75, 2226.75, 2227.75, 2228.75, 2229.75, 2230.75, 2231.75, 2232.75,
           2233.75, 2234.75, 2235.75, 2236.75, 2237.75, 2238.75, 2239.75, 2240.75,
           2241.75, 2242.75, 2243.75, 2244.75, 2245.75, 2246.75, 2247.75, 2248.75],
          dtype=float32)

• Data variables: (10)
  • dayl
    (time, y, x)
    float32
    dask.array<chunksize=(365, 37, 40), meta=np.ndarray>

    long_name :

    daylength

    units :

    s/day

    grid_mapping :

    lambert_conformal_conic

    cell_methods :

    area: mean

    _ChunkSizes :

    [ 16777216 -402456576 -402456576]

    crs :

    epsg:4326

    nodatavals :

    0.0

    Array Chunk Bytes 20.61 MiB 2.06 MiB Shape (3650, 37, 40) (365, 37, 40) Count 30 Tasks 10 Chunks Type float32 numpy.ndarray

  • lat
    (time, y, x)
    float64
    dask.array<chunksize=(365, 37, 40), meta=np.ndarray>

    units :

    degrees_north

    long_name :

    latitude coordinate

    standard_name :

    latitude

    _CoordinateAxisType :

    Lat

    crs :

    epsg:4326

    nodatavals :

    0.0

    Array Chunk Bytes 41.21 MiB 4.12 MiB Shape (3650, 37, 40) (365, 37, 40) Count 40 Tasks 10 Chunks Type float64 numpy.ndarray

  • lon
    (time, y, x)
    float64
    dask.array<chunksize=(365, 37, 40), meta=np.ndarray>

    units :

    degrees_east

    long_name :

    longitude coordinate

    standard_name :

    longitude

    _CoordinateAxisType :

    Lon

    crs :

    epsg:4326

    nodatavals :

    0.0

    Array Chunk Bytes 41.21 MiB 4.12 MiB Shape (3650, 37, 40) (365, 37, 40) Count 40 Tasks 10 Chunks Type float64 numpy.ndarray

  • prcp
    (time, y, x)
    float32
    dask.array<chunksize=(365, 37, 40), meta=np.ndarray>

    long_name :

    daily total precipitation

    units :

    mm/day

    grid_mapping :

    lambert_conformal_conic

    cell_methods :

    area: mean time: sum

    _ChunkSizes :

    [ 16777216 -402456576 -402456576]

    crs :

    epsg:4326

    nodatavals :

    0.0

    Array Chunk Bytes 20.61 MiB 2.06 MiB Shape (3650, 37, 40) (365, 37, 40) Count 30 Tasks 10 Chunks Type float32 numpy.ndarray

  • srad
    (time, y, x)
    float32
    dask.array<chunksize=(365, 37, 40), meta=np.ndarray>

    long_name :

    daylight average incident shortwave radiation

    units :

    W/m2

    grid_mapping :

    lambert_conformal_conic

    cell_methods :

    area: mean time: mean

    _ChunkSizes :

    [ 16777216 -402456576 -402456576]

    crs :

    epsg:4326

    nodatavals :

    0.0

    Array Chunk Bytes 20.61 MiB 2.06 MiB Shape (3650, 37, 40) (365, 37, 40) Count 30 Tasks 10 Chunks Type float32 numpy.ndarray

  • tmax
    (time, y, x)
    float32
    dask.array<chunksize=(365, 37, 40), meta=np.ndarray>

    long_name :

    daily maximum temperature

    units :

    degrees C

    grid_mapping :

    lambert_conformal_conic

    cell_methods :

    area: mean time: maximum

    _ChunkSizes :

    [ 16777216 -402456576 -402456576]

    crs :

    epsg:4326

    nodatavals :

    0.0

    Array Chunk Bytes 20.61 MiB 2.06 MiB Shape (3650, 37, 40) (365, 37, 40) Count 30 Tasks 10 Chunks Type float32 numpy.ndarray

  • tmin
    (time, y, x)
    float32
    dask.array<chunksize=(365, 37, 40), meta=np.ndarray>

    long_name :

    daily minimum temperature

    units :

    degrees C

    grid_mapping :

    lambert_conformal_conic

    cell_methods :

    area: mean time: minimum

    _ChunkSizes :

    [ 16777216 -402456576 -402456576]

    crs :

    epsg:4326

    nodatavals :

    0.0

    Array Chunk Bytes 20.61 MiB 2.06 MiB Shape (3650, 37, 40) (365, 37, 40) Count 30 Tasks 10 Chunks Type float32 numpy.ndarray

  • vp
    (time, y, x)
    float32
    dask.array<chunksize=(365, 37, 40), meta=np.ndarray>

    long_name :

    daily average vapor pressure

    units :

    Pa

    grid_mapping :

    lambert_conformal_conic

    cell_methods :

    area: mean time: mean

    _ChunkSizes :

    [ 16777216 -402456576 -402456576]

    crs :

    epsg:4326

    nodatavals :

    0.0

    Array Chunk Bytes 20.61 MiB 2.06 MiB Shape (3650, 37, 40) (365, 37, 40) Count 30 Tasks 10 Chunks Type float32 numpy.ndarray

  • elevation
    (y, x)
    float64
    dask.array<chunksize=(37, 40), meta=np.ndarray>

    units :

    meters

    crs :

    epsg:4326

    nodatavals :

    0.0

    Array Chunk Bytes 11.56 kiB 11.56 kiB Shape (37, 40) (37, 40) Count 41 Tasks 1 Chunks Type float64 numpy.ndarray

  • pet
    (time, y, x)
    float64
    dask.array<chunksize=(365, 37, 40), meta=np.ndarray>

    units :

    mm/day

    crs :

    epsg:4326

    nodatavals :

    0.0

    Array Chunk Bytes 41.21 MiB 4.12 MiB Shape (3650, 37, 40) (365, 37, 40) Count 30 Tasks 10 Chunks Type float64 numpy.ndarray

• Attributes: (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 Daymet data citation information

  references :

  Please see http://daymet.ornl.gov/ for current information on Daymet references

  History :

  Translated to CF-1.0 Conventions by Netcdf-Java CDM (CFGridWriter2) Original Dataset = /daacftp/daymet/Daymet_Daily_V4/data/daymet_v4_daily_na_dayl_2006.nc; Translation Date = 2021-05-26T07:00:14.233Z

  geospatial_lat_min :

  44.964817634118106

  geospatial_lat_max :

  45.5636129082089

  geospatial_lon_min :

  -69.95790349094842

  geospatial_lon_max :

  -69.13309404831536

  crs :

  +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=km +no_defs

  nodatavals :

  0.0

  transform :

  [ 1.00000e+00 0.00000e+00 2.20625e+03 0.00000e+00 -1.00000e+00 7.29500e+02 0.00000e+00 0.00000e+00 1.00000e+00]

  res :

  [1. 1.]

  bounds :

  [2209.24651856 682.07851575 2250.2123966 719.97190299]

clm.elevation.plot()

<matplotlib.collections.QuadMesh at 0x1b8c9a980>

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

Aridity

variables = {"pet": "Mean Annual PET", "prcp": "Mean Annual P"}
data = {}
for v, n in variables.items():
    data[v] = clm[v].groupby("time.year").sum().mean(dim="year")
    data[v] = data[v].where(data[v] > 0)
    data[v] = data[v].rename(n)
    data[v] = data[v].assign_attrs(unit="mm/year")

aridity = data["pet"] / data["prcp"]
aridity = aridity.rename("Aridity Index")

ax = aridity.plot(size=6)
ax.figure.savefig(Path("_static", "aridity.png"), dpi=300, bbox_inches="tight", facecolor="w")