EQUATES Evaluation

README.md

EQUATES Evaluation Quick Hits

equates_trend_eval.py

__version__ = "0.1.0"
__doc__ = """
Script to Compare EQUATES Hourly Predictions to AQS Observations
----------------------------------------------------------------

---
  author: Barron H. Henderson
  last updated: 2025-05-22
---

This RSIG script is designed specifically for comparing EQUATES to AQS, but
can be adapted for any hourly data sources. The default figure produced is
a yearly average for comparison.

Requires Python 3.9+ with pyrsig, pyproj, pandas, xarray, and netcdf4

```bash
python -m pip install pyrsig pyproj pandas xarray netcdf4
```

"""
# User Input
# ----------
# wdir : str
#   name for region defined by bbox
# bbox : tuple
#   swlon, swlat, nelon, nelat in degrees east and north of 0,0
# ovkey : str
#   name of species (ozone or no2)
# years : list
#    years to retrieve
# months : list
#   months within each year to retrieve
# lsthours : list
#   lst hours to include in plot
wdir = 'BakersfieldCA'
bbox = (-119.21, 35.22, -118.80, 35.44)
ovkey = 'ozone'  # ozone or no2
years = list(range(2005, 2020))
months = [6, 7, 8]
lsthours = [10, 11, 12, 13, 14, 15, 16, 17, 18, 19]

# Import Libraries
# ----------------

import pyrsig
import pandas as pd
import pyproj
import xarray as xr
xr.set_options(keep_attrs=True)

# convert user input
# ------------------
# get RSIG names and keys

okey = f'aqs.{ovkey}'
mvkey = ovkey.upper().replace('OZONE', 'O3')
mkey = f'cmaq.equates.conus.aconc.{mvkey}'
ekey = f'EQUATES_{mvkey}'

# Retrieve Data
# -------------
# Get a month of data at a time
# If using a big bbox, consider getting smaller time chunks

adfs = []
edss = []
for YYYY in years:
  for MM in months:
    bdate = pd.to_datetime(f'{YYYY}-{MM:02d}-01T00:00:00')
    edate = bdate + pd.to_timedelta(bdate.daysinmonth * 24 * 3600 - 1, unit='s')
    print(f'\r{bdate:%FT%H%M%S} to {edate:%FTT%H%M%S}', end='.', flush=True)
    r = pyrsig.RsigApi(bbox=bbox, bdate=bdate, edate=edate, gridfit=True, workdir=wdir)
    eds = r.to_ioapi(mkey)
    edss.append(eds)
    adf = r.to_dataframe(okey, unit_keys=False, parse_dates=True)
    adfs.append(adf)

adf = pd.concat(adfs)
eds = xr.concat(edss, dim='TSTEP')

# add equates sampled at measurement sites
# ----------------------------------------
# Convert lon/lat to COL/ROW, and time to TSTEP
# Then sample the model and add prediction to dataframe

proj = pyproj.Proj(eds.crs_proj4)
adf['COL'], adf['ROW'] = proj(adf['LONGITUDE'], adf['LATITUDE'])
adf['TSTEP'] = adf['time'].dt.tz_localize(None)
idx = adf[['ROW', 'COL', 'TSTEP']].to_xarray()
idx['LAY'] = ('LAY',), [1]
adf[ekey] = eds.isel(LAY=0)[mvkey].sel(ROW=idx.ROW, COL=idx.COL, TSTEP=idx.TSTEP, method='nearest').to_dataframe()[mvkey]

# Make an annual plot
# -------------------
# Filter out the local hours of interest (00:00 to 00:59 is LST 0)
# Convert time to approximate local solar

adf['time_lst'] = adf['time'] + pd.to_timedelta(adf['LONGITUDE'] / 15, unit='h')
# Query for hours of interest
plotdf = adf.query(f'time_lst.dt.hour.isin({lsthours})').groupby(pd.Grouper(key='time', freq='YS'))[[ovkey, ekey]].mean()
ax = plotdf.plot(marker='x')
ax.figure.savefig(f'{wdir}/{wdir}_{okey}_vs_{ekey}.png')