Plots from EMBER

ember.py

__version__ = '0.1.0'
__doc__ = """
Demonstrating how to work with EMBER files from doi:10.5281/zenodo.13737754

Requires requests, xarray, pandas, netcdf4, pyproj, pycno, and pyrsig

For example, install with pip:
python -m pip install -qq xarray pandas netcdf4 pyproj pycno pyrsig

Functions
---------
open_url : function
    Open files using a url or keyword (base, zfire, zcafire)


Example
-------

from ember import open_url
import matplotlib.colors as mc
import matplotlib.pyplot as plt
import pyproj
import pycno
import pyrsig

# Open the Basecase (all emissions)
bcf = open_url('base')

# Open the zero-fire case (no fire emissions)
zff = open_url('zfire')
# zcf = open_url('zcafire')

# Create a fires contribution
ff = bcf - zff
for vk, v in ff.data_vars.items():
  v.attrs['long_name'] = 'Fire ' + v.attrs['long_name']

# Make some maps
proj = pyproj.Proj(bcf.crs_proj4)
clist = [
    "#FFFFFF", "#A020F0", "#00B2EE", "#00FF00", "#FFFF00", "#FFA500",
    "#FF0000", "#A52A2A"
]
cmap = mc.LinearSegmentedColormap.from_list('kpc', clist, 50)
dates = ['2023-05-25', '2023-05-26']
for date in dates:
  bcZ = bcf['ATOTIJ_AVG'].sel(TSTEP=date)
  ffZ = ff['ATOTIJ_AVG'].sel(TSTEP=date)
  aqdf = pyrsig.to_dataframe('aqs.pm25_daily_average', bdate=date)
  aqdf['x'], aqdf['y'] = proj(aqdf['LONGITUDE(deg)'], aqdf['LATITUDE(deg)'])
  fig, axx = plt.subplots(2, 1, figsize=(6, 8))
  qm = bcZ.plot(vmin=0, vmax=30, cmap=cmap, extend='neither', ax=axx[0])
  axx[0].scatter(aqdf.x, aqdf.y, c=aqdf['pm25_daily_average(ug/m3)'],norm=qm.norm, cmap=qm.cmap)
  qm = ffZ.plot(vmin=0, vmax=20, cmap=cmap, extend='neither', ax=axx[1])
  pycno.cno(proj).drawstates(ax=axx)
  plt.setp(
      axx, facecolor='gainsboro', aspect=1, xlim=(15, 150), ylim=(30, 120),
      xticklabels=[], yticklabels=[], xlabel='', ylabel=''
  )
"""
import pyproj
import pycno
import pyrsig
import xarray as xr
xr.set_options(keep_attrs=True)


def open_url(url):
  """
  Arguments
  ---------
  url : str
    URL for file to download or base, zfire, zcafire

  Returns
  -------
  outf : xarray.Dataset
  """
  import requests
  import os
  # Links for Zenodo files
  root = 'https://zenodo.org/records/13737754/files/'
  root = f'{root}/HR2DAY_CMAQv54_cb6r5_ae7_aq.36US3.35.basecase2023_'
  bcurl =   f'{root}2023_O3_PM25_11Apr2023to29Sep2023.nc?download=1'
  zfurl =   f'{root}nofires_2023_O3_PM25_11Arp2023to29Sep2023.nc?download=1'
  zcafurl = f'{root}noCAfires_2023_O3_PM25_11Apr2023to29Sep2023.nc?download=1'
  if url == 'base':
    url = bcurl
  elif url == 'zfire':
    url = zfurl
  elif url == 'zcafire':
    url = zcafurl
  localpath = url.split('/')[-1].split('?')[0]
  if not os.path.exists(localpath):
    r = requests.get(url, stream=True)
    r.raise_for_status()
    with open(localpath, 'wb') as lf:
      lf.write(r.content)
  outf = pyrsig.open_ioapi(localpath)
  # removes spaces fromin attributes
  for vk, v in outf.data_vars.items():
    for pk, pv in v.attrs.items():
      if isinstance(pv, str):
        v.attrs[pk] = pv.strip()

  return outf

To run, you have to install prerequisites, download ember.py, then use ember.open_url

Install Prerequisites

python -m pip install -qq xarray pandas netcdf4 pyproj pycno pyrsig
wget -N https://gist.githubusercontent.com/barronh/c9daef9c04f5bfcfaecf3fd83a08725c/raw/ember.py

Custom Code

This example makes plots for PM and Fire PM. It is cut and paste from ember.__doc__

from ember import open_url
import matplotlib.colors as mc
import matplotlib.pyplot as plt
import pyproj
import pycno
import pyrsig
# Open the Basecase (all emissions)
bcf = open_url('base')
# Open the zero-fire case (no fire emissions)
zff = open_url('zfire')
# zcf = open_url('zcafire')
# Create a fires contribution
ff = bcf - zff
for vk, v in ff.data_vars.items():
  v.attrs['long_name'] = 'Fire ' + v.attrs['long_name']
# Make some maps
proj = pyproj.Proj(bcf.crs_proj4)
clist = [
    "#FFFFFF", "#A020F0", "#00B2EE", "#00FF00", "#FFFF00", "#FFA500",
    "#FF0000", "#A52A2A"
]
cmap = mc.LinearSegmentedColormap.from_list('kpc', clist, 50)
dates = ['2023-05-25', '2023-05-26']
for date in dates:
  bcZ = bcf['ATOTIJ_AVG'].sel(TSTEP=date)
  ffZ = ff['ATOTIJ_AVG'].sel(TSTEP=date)
  aqdf = pyrsig.to_dataframe('aqs.pm25_daily_average', bdate=date)
  aqdf['x'], aqdf['y'] = proj(aqdf['LONGITUDE(deg)'], aqdf['LATITUDE(deg)'])
  fig, axx = plt.subplots(2, 1, figsize=(6, 8))
  qm = bcZ.plot(vmin=0, vmax=30, cmap=cmap, extend='neither', ax=axx[0])
  axx[0].scatter(aqdf.x, aqdf.y, c=aqdf['pm25_daily_average(ug/m3)'],norm=qm.norm, cmap=qm.cmap)
  qm = ffZ.plot(vmin=0, vmax=20, cmap=cmap, extend='neither', ax=axx[1])
  pycno.cno(proj).drawstates(ax=axx)
  plt.setp(
      axx, facecolor='gainsboro', aspect=1, xlim=(15, 150), ylim=(30, 120),
      xticklabels=[], yticklabels=[], xlabel='', ylabel=''
  )

Output from example: