diff --git a/examples/workflows/plot_wildfire_smoke.py b/examples/workflows/plot_wildfire_smoke.py new file mode 100644 index 0000000000..85d3328c0c --- /dev/null +++ b/examples/workflows/plot_wildfire_smoke.py @@ -0,0 +1,176 @@ +""" +Tracking wildfire smoke with AirNow +----------------------------------- + +Canadian wildfire smoke routinely degrades air quality across the US +Midwest and Great Lakes. This example uses ACT's AirNow discovery function +``act.discovery.get_airnow_bounded_obs`` as the single data source to build +a picture of one such event, and plots it entirely with ACT's own displays +combined into one four-panel figure: + +* a regional map of the latest hourly surface PM2.5 at every AirNow station + (``GeographicPlotDisplay``) on top, and +* 24-hour PM2.5 time series for one site each in Duluth MN, Toronto ON, and + Chicago IL (``TimeSeriesDisplay``) on a shared y-axis below, with EPA AQI + category bands shaded behind them. + +Because the US EPA AirNow feed also re-serves Canadian provincial monitors +(for example the Ontario Ministry of the Environment network), a single API +path covers stations on both sides of the border. + +The date window below is fixed to a July 2026 smoke episode so the example is +reproducible; AirNow serves this range from its rolling archive. Set your free +AirNow token (https://docs.airnowapi.org/) in the ``AIRNOW_API`` environment +variable before running. + +Author: Scott Collis +Author: Claude (Anthropic) +""" + +import os + +import cartopy.crs as ccrs +import matplotlib.pyplot as plt +import numpy as np +import xarray as xr +from matplotlib.patches import Patch + +import act + +# --------------------------------------------------------------------------- +# Fixed date window (UTC). AirNow expects 'YYYY-MM-DDTHH'. Pinning the dates +# (rather than using the current time) keeps the example reproducible. +# --------------------------------------------------------------------------- +START_DATE = '2026-07-14T15' +END_DATE = '2026-07-15T16' + +token = os.getenv('AIRNOW_API') + +# Map colorbar ceiling (ug/m^3). +PM25_MAX = 250.0 + +# EPA PM2.5 AQI category bands: (low, high, color, label), ug/m^3. Shaded +# behind the time series so the severity is legible at a glance. +EPA_BANDS = [ + (0.0, 12.0, '#2c7fb8', 'Good'), + (12.0, 35.4, '#fdae61', 'Moderate'), + (35.4, 55.4, '#d7301f', 'Unhealthy (SG)'), + (55.4, 150.4, '#7a0177', 'Unhealthy'), + (150.4, 600.0, '#4d004b', 'Hazardous'), +] + +if token is not None and len(token) > 0: + # ----------------------------------------------------------------------- + # 1. Pull every station in the map region through ACT. data_type='B' + # returns both AQI and concentrations; mon_type=2 includes permanent + # and mobile monitors. The result is a (time, sites) xarray.Dataset. + # ----------------------------------------------------------------------- + map_bounds = '-104,40,-74,50' + ds_map = act.discovery.get_airnow_bounded_obs( + token, START_DATE, END_DATE, map_bounds, 'PM25', mon_type=2, data_type='B' + ) + + # Reduce to the latest valid PM2.5 per site. AirNow flags missing as -999. + pm = ds_map['PM2.5'].values.copy() + pm[pm < 0] = np.nan + latest = np.array( + [col[np.where(~np.isnan(col))[0][-1]] if np.any(~np.isnan(col)) else np.nan for col in pm.T] + ) + + ds_latest = xr.Dataset( + {'PM2.5': ('sites', latest)}, + coords={ + 'latitude': ('sites', ds_map['latitude'].values), + 'longitude': ('sites', ds_map['longitude'].values), + }, + ) + ds_latest['PM2.5'].attrs = {'long_name': 'PM2.5', 'units': 'ug/m^3'} + + # ----------------------------------------------------------------------- + # 2. Pull a tight bounding box around one site in each of three cities. + # ----------------------------------------------------------------------- + cities = { + 'Duluth, MN': ('-92.3,46.6,-91.9,47.0', 'West Duluth'), + 'Toronto, ON': ('-79.65,43.58,-79.25,43.85', 'Toronto Downtown'), + 'Chicago, IL': ('-87.9,41.7,-87.5,42.05', 'CHI_COM'), + } + + ds_cities = {} + for city, (bounds, site) in cities.items(): + ds = act.discovery.get_airnow_bounded_obs( + token, START_DATE, END_DATE, bounds, 'PM25', mon_type=2, data_type='B' + ) + ds_cities[city] = ds.sel(sites=site) + + # ----------------------------------------------------------------------- + # 3. Build one four-panel figure. ACT's GeographicPlotDisplay and + # TimeSeriesDisplay each normally own a whole figure, so we place both + # on a single matplotlib figure by hand: geoplot draws onto the current + # figure and returns its GeoAxes (plus a colorbar), which we reposition + # into the top band; the three shared-y series go in a GridSpec row + # below. geoplot calls plt.axes() on the *current* figure, so we build + # the display first, then create our figure so it is the current one. + # ----------------------------------------------------------------------- + geo = act.plotting.GeographicPlotDisplay(ds_latest) + plt.close('all') # drop the stray figure the ctor made + fig = plt.figure(figsize=(13, 10)) # our figure is now the current one + + map_ax = geo.geoplot( + 'PM2.5', + lat_field='latitude', + lon_field='longitude', + projection=ccrs.PlateCarree(), + cartopy_feature=['LAND', 'OCEAN', 'LAKES', 'STATES', 'BORDERS'], + title='US + Canada surface PM2.5 via AirNow (data pulled with ACT)', + cmap='turbo', + vmin=0, + vmax=PM25_MAX, + marker='o', + s=40, + edgecolor='white', + linewidth=0.3, + ) + map_ax.set_extent([-104, -74, 39.5, 50.5], crs=ccrs.PlateCarree()) + + # geoplot added [GeoAxes, colorbar]; move both into the top band. + cbar_ax = [a for a in fig.axes if a is not map_ax][0] + map_ax.set_position([0.06, 0.44, 0.82, 0.52]) + cbar_ax.set_position([0.90, 0.47, 0.018, 0.44]) + + # Bottom row: three shared-y series with EPA category bands behind them. + gs = fig.add_gridspec(1, 3, left=0.06, right=0.97, bottom=0.14, top=0.34, wspace=0.12) + bottom_axes = [fig.add_subplot(gs[0, i]) for i in range(3)] + + tsd = act.plotting.TimeSeriesDisplay(ds_cities, subplot_shape=(1, 3)) + plt.close(tsd.fig) # discard its throwaway figure + tsd.fig = fig # point it at our figure and axes + tsd.axes = np.array([bottom_axes]) + for idx, (city, (_, site)) in enumerate(cities.items()): + ax = bottom_axes[idx] + for lo, hi, color, _ in EPA_BANDS: + ax.axhspan(lo, hi, color=color, alpha=0.13, zorder=0, linewidth=0) + tsd.plot( + 'PM2.5', + dsname=city, + subplot_index=(0, idx), + y_rng=(0, 600), + force_line_plot=True, + marker='o', + set_title=f'{city} \u2014 {site}', + ) + if idx > 0: # label the y-axis once + ax.set_ylabel('') + ax.tick_params(labelleft=False) + + handles = [Patch(facecolor=c, alpha=0.5, label=lbl) for _, _, c, lbl in EPA_BANDS] + fig.legend( + handles=handles, + loc='lower center', + ncol=5, + frameon=False, + fontsize=9, + bbox_to_anchor=(0.5, 0.02), + title='EPA PM2.5 AQI category (bands)', + ) + + plt.show()