Background
Author: Diane Menuz
Date: 7/23/2026
Goal: Download forcing variables from NLDAS and use to calculate hourly ETo with refet package
Downloads are based on this nasa tutorial
NLDAS data are returned as UTC and refet uses UTC for calculations.
Initial Set Up
Register for an Earthdata account
Log in to account and generate a token
Copy token as parameter into script below; note that tokens are good for a limited period of time
If call_nldas_time_series throws an error, you may need to turn on some apps under the Authorized Apps in your Earthdata account. I experimented with turning on the following but in the end did not need any of them:
Giovanni
NASA GESDISC DATA ARCHIVE
Hyrax in the cloud
New Harmony Giovanni Adapter
Import Libraries
[1]:
import os
import sys
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import refet
sys.path.insert(0, os.path.abspath("../../src")) # adjust path as needed
from fluxfootprints import fetch_nldas_forcing_dataset
Parameters
Copy your token from Earthdata into the parameters below
[2]:
import os
token = os.environ.get('EARTHDATA_TOKEN')
[33]:
# station parameters
# note that the time start and end should be in UTC
lat = 37.7353
lon = -111.5708
elevation_m = 1729.3
time_start = '2021-06-01'
time_end = '2026-10-01'
utc_offset = -7
STATIONID = 'US-UTE'
Download NLDAS Data and Format
[34]:
# format datetime and d
time_start = pd.to_datetime(time_start)
time_start = time_start.strftime("%Y-%m-%dT%H:%M:%S")
time_end = pd.to_datetime(time_end)
time_end = time_end.strftime("%Y-%m-%dT%H:%M:%S")
[35]:
# obtain nldas data for site and date range
nldas_raw = fetch_nldas_forcing_dataset(lat, lon, time_start, time_end, token)
nldas_raw['timestamp_local'] = nldas_raw.index + pd.Timedelta(hours=utc_offset)
Fetching temp_K...
Status Code: 200
Fetching spec_hum...
Status Code: 200
Fetching pressure_pa...
Status Code: 200
Fetching wind_u10...
Status Code: 200
Fetching wind_v10...
Status Code: 200
Fetching solar_rad...
Status Code: 200
[36]:
# convert variables as needed and calculate hourly ref et
df = nldas_raw.copy()
# Compute 10m wind speed magnitude from NLDAS zonal/meridional vectors
df["wind_10m"] = np.sqrt(df["wind_u10"] ** 2 + df["wind_v10"] ** 2)
# Compute actual vapor pressure (ea) in kPa from specific humidity (q) & surface pressure (P)
p_kpa = df["pressure_pa"] / 1000.0
df["ea_kpa"] = (df["spec_hum"] * p_kpa) / (
0.622 + (1.0 - 0.622) * df["spec_hum"]
)
# run hourly model
refet_model = refet.Hourly(
tmean=df["temp_K"], # Temperature in Kelvin
ea=df["ea_kpa"], # Actual vapor pressure in kPa
rs=df["solar_rad"], # Solar radiation in W/m^2
uz=df["wind_10m"], # Wind speed in m/s at 10m
zw=10, # Tell refet wind is measured at 10m
elev=elevation_m, # Elevation in meters
lat=lat, # Latitude in degrees
lon=lon, # Longitude in degrees
doy=df.index.dayofyear, # Day of year (1-366)
time=df.index.hour, # UTC hour (0-23)
method="asce",
input_units={
"tmean": "K", # Kelvin -> refet converts to C internally
"rs": "w/m2", # W/m2 -> refet converts to MJ/m2/hr internally
"lat": "deg", # Degrees -> refet converts to rad internally
"uz": "m/s",
},
)
# xtract short grass (ETo) in mm/hr
df["eto_grass"] = refet_model.eto()
# Compute each hour's fraction relative to the daily total
daily_total_eto = df.groupby(df.timestamp_local.dt.date)["eto_grass"].transform("sum")
df["eto_hourly_fraction"] = df["eto_grass"] / daily_total_eto
#df["timestamp_local"] = df.index#.to_series().dt.strftime("%Y-%m-%d %H:%M:%S")
df["timestamp_end"] = df["timestamp_local"] + pd.Timedelta(hours=1)
NLDAS_ET_PATH= f"M:/Shared drives/UGS_Flux/Data_Processing/final_processing_workflow/footprints/nldas/nldas_results/{STATIONID}_nldas.csv"
df.to_csv(NLDAS_ET_PATH + f"{STATIONID}_nldas.csv", index=True)
[22]:
# plot a single day- make sure time alignment is correct
second_date = df["timestamp_local"].dt.date.unique()[1]
one_day = df[df["timestamp_local"].dt.date == second_date]
one_day.plot(x="timestamp_local", y="eto_grass", grid=True, marker="o")
plt.title("Hourly Reference ET (mm/hr)")
plt.ylabel("ETo [mm/hr]")
plt.xlabel("Local Time")
plt.xticks(rotation=45) # Rotates timestamps so they don't overlap
plt.tight_layout()
plt.show()
[37]:
STATIONID = 'US-UTE'
NLDAS_ET_PATH= f"M:/Shared drives/UGS_Flux/Data_Processing/final_processing_workflow/footprints/nldas/nldas_results/{STATIONID}_nldas.csv"
pd.read_csv(NLDAS_ET_PATH).head()
[37]:
| temp_K | spec_hum | pressure_pa | wind_u10 | wind_v10 | solar_rad | wind_10m | ea_kpa | eto_grass | timestamp_local | eto_hourly_fraction | timestamp_end | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 295.18 | 0.003043 | 81664.4 | 0.77 | 1.65 | 677.926 | 1.820824 | 0.398788 | 0.466893 | 2024-06-18 16:00:00 | 0.401319 | 2024-06-18 17:00:00 |
| 1 | 295.43 | 0.003050 | 81644.8 | 1.11 | 1.95 | 475.875 | 2.243791 | 0.399608 | 0.341527 | 2024-06-18 17:00:00 | 0.293560 | 2024-06-18 18:00:00 |
| 2 | 294.06 | 0.003275 | 81629.1 | 1.27 | 1.80 | 297.647 | 2.202930 | 0.428998 | 0.210657 | 2024-06-18 18:00:00 | 0.181071 | 2024-06-18 19:00:00 |
| 3 | 292.69 | 0.003501 | 81613.6 | 1.42 | 1.66 | 126.529 | 2.184491 | 0.458462 | 0.071179 | 2024-06-18 19:00:00 | 0.061182 | 2024-06-18 20:00:00 |
| 4 | 291.32 | 0.003728 | 81598.2 | 1.58 | 1.51 | 0.000 | 2.185521 | 0.487998 | 0.031942 | 2024-06-18 20:00:00 | 0.027456 | 2024-06-18 21:00:00 |