Computing PET
This page was generated from pet_validation.ipynb.
Interactive online version:
Computing PET#
[1]:
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
[2]:
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.
[3]:
station = "01013500"
coords = (-68.58264, 47.23739)
dates = ("2000-01-01", "2009-01-12")
root = Path("input_data")
[4]:
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)
[5]:
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.
[6]:
camels = pd.read_parquet(Path(root, f"{station}_05_model_output.parquet"))
camels = camels.loc[clm_c.index, "PET"]
style["CAMELS"] = "-"
[7]:
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.
[8]:
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")
[8]:
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.
[9]:
pet.plot(figsize=(12, 6), ylabel="PET (mm/day)", style=style)
[9]:
<AxesSubplot:xlabel='time', ylabel='PET (mm/day)'>

[10]:
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)
[11]:
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"
)
