Compare AirNow and NAQFC

airnowcompare.py

__version__ = "1.0.0"
__all__ = ['open_remote', 'open_naqfc', 'open_airnow', 'open_airfuse']
__doc__ = """
Utility for Comparing AirNow Surfaces
=====================================

Using AirNow grib2 AQI files as reference. Can compare to NAQFC from NOAA or
the AirFuse from AirNow. For simplicity, NAQFC and AirFuse are converted to
a Nowcast by averaging the last 3 time steps. The Nowcast is then put into AQI
units for comparison.

Requires:
  python>=3.7
  matplotlib, pandas, xarray, cfgrib, pyproj

Example
-------
  #wget -N https://gist.github.com/barronh/945df1b465905b0f5454258b317cee1b/raw/airnowcompare.py
  #python -m pip install netcdf4 cfgrib eccodes pycno
  from airnowcompare import open_airnow, open_airfuse, compare
  import pandas as pd
  date = pd.to_datetime('2025-04-23T16Z')
  spc = 'o3'
  anf = open_airnow(date, spc)
  othf = open_airfuse(date, spc, aqi=True)
  fig = compare(anf[spc], othf[spc])
  fig.savefig('test.png')
"""
import pyproj
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import pycno

aqicolors = ['green', 'yellow', 'orange', 'red', 'purple', 'maroon']
aqiedges = [0, 50.5, 100.5, 150.5, 200.5, 300.5]
edges = {
    'pm25': [0, 9.1, 35.5, 54.5, 125.5, 225.5],
    'o3': [0, 55, 71, 86, 106, 201],
}

naqfccrs = pyproj.CRS.from_cf({
  'grid_mapping_name': 'lambert_conformal_conic',
  'latitude_of_projection_origin': 25.0,
  'longitude_of_central_meridian': 265.0 - 360,
  'standard_parallel': 25.0,
  'earth_radius': 6371229.0
})
naqfcproj4 = '+proj=lcc +lat_1=25 +lat_0=25 +lon_0=-95 +k_0=1 +x_0=0 +y_0=0 +R=6371229 +to_meter=1000 +no_defs'
naqfcproj = pyproj.Proj(naqfcproj4)


def naqfc2airnow(nq, an):
  import numpy as np
  import xarray as xr

  nq = nq.copy()
  out = an * np.nan
  if 'x' not in nq.coords:
    cx, cy = naqfcproj(nq.longitude, nq.latitude)
    nq.coords['x'] = cx.mean(0)
    nq.coords['y'] = cy.mean(1)
  if 'longitude' in nq.coords:
    del nq.coords['longitude']
    del nq.coords['latitude']
  LAT, LON = xr.broadcast(an.latitude, an.longitude)
  x, y = naqfcproj(LON, LAT)
  X = xr.full_like(an, fill_value=np.nan)
  X[:] = x
  Y = xr.full_like(an, fill_value=np.nan)
  Y[:] = y
  out = nq.interp(x=X, y=Y, method='linear')
  return out


def toaqi(Z, spc=None):
  """
  Arguments
  ---------
  Z : DataArray

  Returns
  -------
  out : same as Z
    Interpolated to aqi number space
  """
  import numpy as np
  out = Z * np.nan
  spc = spc or getattr(Z, 'name')
  out[:] = np.interp(Z, edges[spc], aqiedges)
  return out


def open_remote(url, dest=None):
  """
  Download a url to dest and open as an xarray dataset.

  Arguments
  ---------
  url : str
    url to data (e.g., https://...)
  dest : str
    Defaults to last element of url delimited by /
  
  Returns
  -------
  f : xarray.Dataset
  """
  import os
  import requests
  import xarray as xr  
  if dest is None:
    dest = url.split('/')[-1]
  if not os.path.exists(dest):
    os.makedirs(os.path.dirname(dest), exist_ok=True)
    print('downloading...', end='',  flush=True)
    with open(dest, 'wb') as outf, requests.get(url) as r:
      outf.write(r.content)
  print('opening...', end='', flush=True)
  f = xr.open_dataset(dest, decode_times=False)
  if 'longitude' in f.coords:
    lon = f.coords['longitude'] - 360
    f.coords['longitude'] = lon
  else:
    Y, X = xr.broadcast(f.y, f.x)
    lon, lat = naqfcproj(X, Y, inverse=True)
    f.coords['longitude'] = ('y', 'x'), lon
    f.coords['latitude'] = ('y', 'x'), lat
  return f


def open_airnow(date, spc):
  """
  Open an AirNow interpolated surface.

  Arguments
  ---------
  date : date-like
    date including utc hour
  spc : str
    pm25, pm10 or ozone
  
  Returns
  -------
  f : xarray.Dataset
    nowcast-like dataset with spc as a variable in AQI units
  """
  if spc in ('o3', 'ozone'):
    sfx = ''
  else:
    sfx = '_' + spc.lower()
  dest = f'{date:%Y/%Y%m%d/US-%y%m%d%H}{sfx}.grib2'
  root = 'https://files.airnowtech.org/airnow'
  f = open_remote(f'{root}/{dest}', dest)
  key = list(f.data_vars)[0]
  f = f.rename(**{key: spc})
  f[spc].attrs['long_name'] = f'AirNow {spc}'
  return f


def open_naqfc(date, spc, bc=True, aqi=True):
  """
  Open an NAQFC interpolated surface using the average of the last
  three times as an approximation of the Nowcast. If aqi is True, then
  return the spc in aqi units.

  Arguments
  ---------
  date : date-like
    date including utc hour
  spc : str
    pm25 or ozone
  bc : bool
    Use KFAN bias correction
  
  Returns
  -------
  f : xarray.Dataset
    nowcast-like dataset with spc as a variable
  """
  fdate = date - pd.to_timedelta(3, unit='h')
  while fdate.hour not in (6, 12):
    fdate = fdate - pd.to_timedelta(1, unit='h')
  sfx = {True: '_bc', False: ''}[bc]
  root = 'https://noaa-nws-naqfc-pds.s3.amazonaws.com'
  dest = f'AQMv7/CS/{fdate:%Y%m%d/%H/aqm.t%Hz}.ave_1hr_{spc}{sfx}.{fdate:%Y%m%d}.227.grib2'
  f = open_remote(f'{root}/{dest}', dest).rename(time='reftime', step='time')
  key = list(f.data_vars)[0]
  f = f.rename(**{key: spc})
  reftime = pd.to_datetime('1970-01-01T00+0000')
  f.coords['time'] = reftime + pd.to_timedelta(f.valid_time.data, unit='s')
  f = f.sel(time=slice(date - pd.to_timedelta('3h'), date)).mean('time')
  f[spc] = toaqi(f[spc], spc=spc)
  f[spc].attrs['long_name'] = f'NAQFC {spc}'
  return f


def open_airfuse(date, spc, aqi=True):
  root = 'https://airnow-navigator-layers.s3.us-east-2.amazonaws.com'
  if spc in ('o3', 'ozone'):
    dest = f'{date:fusion/O3/%Y/%m/%d/%H/Fusion_O3_NAQFC_airnow_%Y-%m-%dT%H}Z.nc'
    key = 'aVNA'
  else:
    dest = f'{date:fusion/PM25/%Y/%m/%d/%H/Fusion_PM25_NAQFC_%Y-%m-%dT%H}Z.nc'
    key = 'FUSED_aVNA'
  f = open_remote(f'{root}/{dest}', dest)
  f = f.rename(**{key: spc})
  if aqi:
    h = pd.to_timedelta('1h')
    fm1 = open_airfuse(date - h, spc, aqi=False)
    fm2 = open_airfuse(date - h - h, spc, aqi=False)
    f = xr.concat([f, fm1, fm2], dim='time').mean('time')
    f[spc] = toaqi(f[spc], spc=spc)
  f[spc].attrs['long_name'] = f'AirFuse {spc}'
  return f


def compare(an, nq, diff=None):
  """
  Arguments
  ---------
  an : DataArray
  nq : DataArray
  diff : DataArray
  Returns
  -------
  """
  if diff is None:
    diff = naqfc2airnow(nq, an) - an
    diff.attrs['long_name'] = f'{nq.long_name} - {an.long_name}'
  cmap, norm = plt.matplotlib.colors.from_levels_and_colors(aqiedges, aqicolors, extend='max')
  gskw = dict(left=.05, right=.98)
  fig, axx = plt.subplots(1, 3, figsize=(16, 4), gridspec_kw=gskw)
  qm = axx[0].pcolormesh(an.longitude, an.latitude, an, norm=norm, cmap=cmap)
  fig.colorbar(qm, label=f'{an.long_name} AQI')
  qm = axx[1].pcolormesh(nq.longitude, nq.latitude, nq, norm=norm, cmap=cmap)
  fig.colorbar(qm, label=f'{nq.long_name}')
  dnorm = plt.matplotlib.colors.TwoSlopeNorm(0)
  dcmap = 'seismic'
  qm = axx[2].pcolormesh(diff.longitude, diff.latitude, diff, norm=dnorm, cmap=dcmap)
  axx[2].set(facecolor='gainsboro')
  fig.colorbar(qm, label=f'{diff.long_name}')
  pycno.cno().drawstates(ax=axx, resnum=1)
  return fig


if __name__ == "__main__":
  import argparse
  prsr = argparse.ArgumentParser()
  prsr.add_argument('--figpath', default=None)
  prsr.add_argument('--other', choices={'naqfc', 'airfuse'}, default='airfuse')
  prsr.add_argument('spc', choices={'o3', 'pm25', 'pm25', 'pm10'})
  prsr.add_argument('date', help='date', type=pd.to_datetime)
  args = prsr.parse_args(['o3', '2025-04-24T16+0000'])
  anf = open_airnow(args.date, args.spc)
  if args.other == 'naqfc':
    nqf = open_naqfc(args.date, args.spc, aqi=True)
  elif args.other == 'airfuse':
    nqf = open_airfuse(args.date, args.spc, aqi=True)
  fig = compare(anf[args.spc], nqf[args.spc])
  if args.figpath is None:
    datestr = args.date.strftime('%FT%HZ')
    args.figpath = f'airnow_vs_{args.other}_{datestr}.png'
  fig.savefig(args.figpath)

Install module and prerequisites:

wget -N https://gist.github.com/barronh/945df1b465905b0f5454258b317cee1b/raw/airnowcompare.py
python -m pip install netcdf4 cfgrib eccodes pycno pyrsig

Compare AirFuse to AirNow for ozone:

import matplotlib.pyplot as plt
import pyrsig
import pandas as pd
from airnowcompare import open_airnow, open_airfuse, open_naqfc, compare
from airnowcompare import aqiedges, aqicolors, toaqi

datestr = '2025-04-24T16Z'
date = pd.to_datetime(datestr)
spc = 'pm25'
anf = open_airnow(date, spc)
aff = open_airfuse(date, spc, aqi=True)
nqf = open_naqfc(date, spc, aqi=True)
# Get observations from RSIG
df = pyrsig.to_dataframe('airnow.pm25', bdate=date - pd.to_timedelta('3.25h'), edate=date + pd.to_timedelta('3599s'), unit_keys=False)
# Average them to a location (i.e., monitor)
mondf = df.groupby(['LONGITUDE', 'LATITUDE'], as_index=False).mean(numeric_only=True)
mondf['aqi'] = toaqi(mondf['pm25'], 'pm25')

# Add them to Figure 2
cmap, norm = plt.matplotlib.colors.from_levels_and_colors(aqiedges, aqicolors, extend='max')


fig1 = compare(anf[spc], aff[spc])
wndw = dict(xlim=(-76, -72), ylim=(38, 42))
plt.setp(fig1.axes[:3], **wndw)
fig2 = compare(anf[spc], nqf[spc])
plt.setp(fig2.axes[:3], **wndw)
for fig in [fig1, fig2]:
  for ax in fig.axes[:2]:
    mondf.plot.scatter(x='LONGITUDE', y='LATITUDE', c='aqi', edgecolor='k', ax=ax, norm=norm, cmap=cmap, colorbar=False)
    ax.set(xlabel='', ylabel='')
fig1.savefig(f'airnow_vs_airfuse_{spc}_{datestr}.png')
fig2.savefig(f'airnow_vs_naqfc_{spc}_{datestr}.png')