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

Computing PET

import calendar
from pathlib import Path

import pandas as pd
import rioxarray as rxr
from hymod import compute_kge

import pydaymet as daymet
from pynhd import WaterData

import warnings

warnings.filterwarnings("ignore", message=".*Index.*")

PyDaymet offers three methods for computing PET: penman_monteith, priestley_taylor, and hargreaves_samani. Let’s compute PET using PyDaymet with these three methods and compare them with the CAMELS.

We choose station ID of 01013500 and time period of 2000-01-01 to 2009-01-12.

station = "01013500"
coords = (-68.58264, 47.23739)
dates = ("2000-01-01", "2009-01-12")
root = Path("input_data")

clm_par = Path(root, f"{station}_clm_single_pixel.parquet")
clm_nc = Path(root, f"{station}_clm_basin.nc")

if clm_par.exists() and clm_nc.exists():
    clm_c = pd.read_parquet(clm_par)
    clm_b = rxr.open_rasterio(clm_nc)
else:
    wd = WaterData("gagesii_basins")
    basin = wd.byid("gage_id", station)
    clm_c = daymet.get_bycoords(coords, dates)
    clm_b = daymet.get_bygeom(basin.geometry[0], dates)
    clm_c.to_parquet(clm_par)
    clm_b.to_netcdf(clm_nc)

single, gridded, style = {}, {}, {}
methods = {"HS": "hargreaves_samani", "PM": "penman_monteith", "PT": "priestley_taylor"}
for n, m in methods.items():
    pet_c = daymet.potential_et(clm_c, coords, method=m)
    single[f"{n} (Single Pixel)"] = pet_c["pet (mm/day)"]
    style[f"{n} (Single Pixel)"] = "-"
    pet_b = daymet.potential_et(clm_b, method=m)
    gridded[f"{n} (Areal Average)"] = pet_b.pet.mean(dim=["x", "y"]).values.flatten()
    style[f"{n} (Areal Average)"] = "--"

Now, let’s get the CAMELS dataset and extract PET for our target station and time period.

camels = pd.read_parquet(Path(root, f"{station}_05_model_output.parquet"))
camels = camels.loc[clm_c.index, "PET"]
style["CAMELS"] = "-"

pet = pd.DataFrame({**single, **gridded, "CAMELS": camels})
pet[pet.lt(0)] = 0.0

Let’s compute Kling-Gupta Efficiency (KGE) for each method in comparison to the CAMELS’ results.

kge = {
    m: compute_kge(pet["CAMELS"].to_numpy("f8"), pet[m].to_numpy("f8"))
    for m in pet
    if m != "CAMELS"
}
pd.Series(kge, name="Kling-Gupta Efficiency")

HS (Single Pixel)     0.926985
PM (Single Pixel)     0.857869
PT (Single Pixel)     0.904256
HS (Areal Average)    0.914182
PM (Areal Average)    0.865927
PT (Areal Average)    0.926194
Name: Kling-Gupta Efficiency, dtype: float64

It seems that PT (Areal Average) performs better than the rest. We can plot the data for daily and mean-monthly PET.

pet.plot(figsize=(12, 6), ylabel="PET (mm/day)", style=style)

<AxesSubplot:xlabel='time', ylabel='PET (mm/day)'>

month_abbr = dict(enumerate(calendar.month_abbr))
mean_month = pet.groupby(pd.Grouper(freq="M")).sum()
mean_month = mean_month.groupby(mean_month.index.month).mean()
mean_month.index = mean_month.index.map(month_abbr)

ax = mean_month.plot(figsize=(12, 6), ylabel="PET (mm/month)", style=style)
ax.figure.savefig(
    Path("_static", "pet_validation.png"), dpi=300, bbox_inches="tight", facecolor="w"
)