Documentation Index
Fetch the complete documentation index at: https://mintlify.com/NASAARSET/PM2.5_AQ_Online_2026/llms.txt
Use this file to discover all available pages before exploring further.
The ARSET PM2.5 training workflow connects three independently sourced datasets — satellite-derived PM2.5 (SatPM), MERRA-2 CNN reanalysis, and ground-level monitor observations — into a unified analysis pipeline. Each stage produces structured outputs that feed directly into the next, and every intermediate result can be saved and reloaded from Google Drive to support multi-session work.
1. Define the Study Domain
All downstream downloads and aggregations are scoped to a single bounding box and time window. Set these variables at the top of your notebook before running any other cells.
# Spatial bounding box
n_lat_min = 18.0
n_lat_max = 32.0
n_lon_min = -92.0
n_lon_max = -75.0
# Time window
t_start = '2023-01-01'
t_end = '2023-12-31'
# Local timezone for diel-cycle and daily-average computations
s_timezome = 'America/Chicago'
# Grid aggregation quality threshold
n_minimum_samples_per_grid = 3
# Threshold for above/below analysis
n_threshold = 35.0 # µg/m³ (US 24-hr NAAQS)
The variable s_timezome (single-z spelling) matches the notebook source. Keep this spelling consistent when referencing it in helper functions so that downstream local_hour_of_day and local_day columns are computed correctly.
2. Download Data
Download SatPM Monthly NetCDF Files
SatPM files are stored as monthly NetCDF archives on S3. The download loop iterates over every calendar month in [t_start, t_end], fetches the corresponding file, and subsets it to the bounding box before writing to disk.
import s3fs, xarray as xr, pandas as pd
s3 = s3fs.S3FileSystem(anon=True)
months = pd.date_range(t_start, t_end, freq='MS')
sat_files = []
for m in months:
s3_key = f"arset-pm25/SatPM/{m.strftime('%Y%m')}_SatPM25.nc"
local_path = f"data/SatPM_{m.strftime('%Y%m')}.nc"
s3.get(s3_key, local_path)
ds = xr.open_dataset(local_path)
ds_sub = ds.sel(lat=slice(n_lat_min, n_lat_max),
lon=slice(n_lon_min, n_lon_max))
ds_sub.to_netcdf(local_path.replace('.nc', '_subset.nc'))
sat_files.append(local_path.replace('.nc', '_subset.nc'))
Output: one subsetted NetCDF per month, e.g., data/SatPM_202301_subset.nc.
Download MERRA-2 CNN Hourly NetCDF via earthaccess
MERRA-2 CNN data are retrieved through the earthaccess library, which handles NASA Earthdata authentication and direct-to-memory streaming. Files are hourly at 0.5° × 0.625° resolution.
import earthaccess
earthaccess.login(strategy='netrc') # or strategy='interactive'
results = earthaccess.search_data(
short_name='M2T1NXAER',
temporal=(t_start, t_end),
bounding_box=(n_lon_min, n_lat_min, n_lon_max, n_lat_max)
)
files = earthaccess.download(results, local_path='data/MERRA2/')
If your Earthdata credentials are already stored in ~/.netrc, strategy='netrc' avoids interactive prompts when running in a scheduled or headless environment.
Output: hourly NetCDF files under data/MERRA2/, one per day.
Download Ground Monitor Data
Ground monitor data can come from two sources depending on the study region.
OpenAQ API (international / US monitors)
import requests, pandas as pd
def fetch_openaq(lat_min, lat_max, lon_min, lon_max, date_from, date_to):
params = dict(
parameter='pm25',
date_from=date_from,
date_to=date_to,
limit=10000,
coordinates=f"{(lat_min+lat_max)/2},{(lon_min+lon_max)/2}",
radius=500000
)
r = requests.get('https://api.openaq.org/v2/measurements', params=params)
r.raise_for_status()
return pd.json_normalize(r.json()['results'])
f_data_monitors_raw = fetch_openaq(n_lat_min, n_lat_max, n_lon_min, n_lon_max,
t_start, t_end)
SINCA Scraper (Chilean monitor network)
# SINCA returns CSV pages; the scraper iterates station IDs
import requests, pandas as pd, io
SINCA_BASE = 'https://sinca.mma.gob.cl/cgi-bin/broker.exe'
def fetch_sinca_station(station_id, year):
params = dict(datatype='D', station=station_id,
parameter='MP25', year=year)
r = requests.get(SINCA_BASE, params=params)
return pd.read_csv(io.StringIO(r.text), sep=';', decimal=',')
frames = [fetch_sinca_station(sid, yr)
for sid in s_station_ids
for yr in pd.date_range(t_start, t_end, freq='YS').year]
f_data_monitors_raw = pd.concat(frames)
Output: f_data_monitors_raw DataFrame with columns including lat, lon, time, and raw PM2.5 values.
3. Resume from Saved Subsets
Long download sessions can be interrupted. Save and reload intermediate DataFrames from Google Drive without re-downloading.
from google.colab import drive
import pickle, os
drive.mount('/content/drive')
SAVE_DIR = '/content/drive/MyDrive/ARSET_PM25/'
os.makedirs(SAVE_DIR, exist_ok=True)
# Save
with open(SAVE_DIR + 'f_data_monitors.pkl', 'wb') as fh:
pickle.dump(f_data_monitors, fh)
# Reload in a new session
with open(SAVE_DIR + 'f_data_monitors.pkl', 'rb') as fh:
f_data_monitors = pickle.load(fh)
Reload cells must be run before any analysis cells. If you skip the reload and f_data_monitors is undefined, downstream grid-snapping and aggregation steps will raise a NameError.
4. Analysis Pipeline
Each analysis step builds on the cleaned monitor DataFrame and the gridded satellite/reanalysis arrays. Run them in order or independently after reloading saved subsets.
4a. Map Comparisons
Overlay gridded satellite PM2.5 fields with point-level monitor observations to visually assess spatial agreement.
F_plot_map(a_data_plot=a_data_SatPM['PM25_filtered'],
a_data_plot_point=f_data_monitors['PM25_filtered'],
n_lat_min=n_lat_min,
n_lat_max=n_lat_max,
n_lon_min=n_lon_min,
n_lon_max=n_lon_max,
Colorbar_Label='PM$_{2.5}$ [µg/m$^{3}$]',
Plot_Title=f'Comparing SatPM to Monitors\n(Average {t_start} to {t_end})',
Colorbar_Limits=[0, 40])
F_plot_map signature and parameter reference.
4b. Monthly Averages
Aggregate hourly monitor observations to monthly means, then look up the corresponding SatPM and MERRA-2 CNN values for each (lat, lon, time_month_start) triplet.
import numpy as np
# Specify the minimum number of valid hourly samples for a monthly average:
# 24*20 = 480 hours (at least 20 days of data per month)
n_minimum_samples_per_month = 24*20
# Create a working copy and stamp each row with the month start:
f_data_monitors_monthly = f_data_monitors.copy().reset_index()
f_data_monitors_monthly['time_month_start'] = [
np.datetime64(t_time).astype('datetime64[M]')
for t_time in f_data_monitors_monthly['time']
]
# Compute monthly count, mean, and standard deviation:
f_data_combined_monthly = (
f_data_monitors_monthly
.groupby(['lat', 'lon', 'time_month_start'])[['PM25_filtered']]
.count()
.rename(columns={'PM25_filtered': 'Count'})
)
f_data_combined_monthly['Average'] = (
f_data_monitors_monthly
.groupby(['lat', 'lon', 'time_month_start'])[['PM25_filtered']]
.mean()
)
f_data_combined_monthly['StDev'] = (
f_data_monitors_monthly
.groupby(['lat', 'lon', 'time_month_start'])[['PM25_filtered']]
.std()
)
# Drop cells that don't meet the minimum sample threshold:
f_data_combined_monthly = f_data_combined_monthly.where(
f_data_combined_monthly['Count'] >= n_minimum_samples_per_month
)
a_data_SatPM and a_data_MERRA2_CNN.
4c. Hourly Grid-Snapped Comparisons
Snap each monitor to the nearest MERRA-2 grid cell (0.5° lat × 0.625° lon) and require at least n_minimum_samples_per_grid measurements per cell-hour before including that cell.
n_minimum_samples_per_grid = 3
f_data_monitors_gridded = f_data_monitors.copy().reset_index()
# Snap to MERRA-2 grid resolution
f_data_monitors_gridded['lat_grid'] = (
round(f_data_monitors_gridded['lat'] / 0.5) * 0.5
)
f_data_monitors_gridded['lon_grid'] = (
round(f_data_monitors_gridded['lon'] / 0.625) * 0.625
)
# Aggregate within each grid cell
f_data_combined_gridded = (
f_data_monitors_gridded
.groupby(['lat_grid', 'lon_grid', 'time'])[['PM25_filtered']]
.count()
.rename(columns={'PM25_filtered': 'Count'})
)
f_data_combined_gridded['Average'] = (
f_data_monitors_gridded
.groupby(['lat_grid', 'lon_grid', 'time'])[['PM25_filtered']]
.mean()
)
f_data_combined_gridded['StDev'] = (
f_data_monitors_gridded
.groupby(['lat_grid', 'lon_grid', 'time'])[['PM25_filtered']]
.std()
)
# Apply quality filter
f_data_combined_gridded = (
f_data_combined_gridded
.where(f_data_combined_gridded['Count'] >= n_minimum_samples_per_grid)
.dropna()
)
The Count threshold prevents isolated single-monitor grid cells from dominating the statistics. Increase n_minimum_samples_per_grid in dense urban networks to require stronger consensus.
4d. Daily Averages
Collapse the hourly gridded DataFrame to daily means for threshold and time-series analysis.
import pytz
tz = pytz.timezone(s_timezome)
f_data_combined_gridded = f_data_combined_gridded.reset_index()
f_data_combined_gridded['local_time'] = (
f_data_combined_gridded['time']
.dt.tz_localize('UTC')
.dt.tz_convert(tz)
)
f_data_combined_gridded['local_day'] = (
f_data_combined_gridded['local_time'].dt.date
)
f_data_combined_gridded_daily = (
f_data_combined_gridded
.groupby(['lat_grid', 'lon_grid', 'local_day'])[['Average', 'MERRA2_CNN']]
.mean()
.reset_index()
)
4e. Threshold Analysis
Split daily averages above and below a regulatory threshold to evaluate model performance separately in polluted vs. clean conditions.
f_above = f_data_combined_gridded_daily[
f_data_combined_gridded_daily['Average'] >= n_threshold
]
f_below = f_data_combined_gridded_daily[
f_data_combined_gridded_daily['Average'] < n_threshold
]
print(f'Days above {n_threshold} µg/m³: {len(f_above)}')
print(f'Days below {n_threshold} µg/m³: {len(f_below)}')
4f. Diel Cycle
Group by local hour of day to reveal the diurnal PM2.5 pattern. See the Visualization guide for the accompanying plot.
f_data_combined_gridded['local_hour_of_day'] = (
f_data_combined_gridded['local_time'].dt.hour
)
f_data_combined_gridded_by_timeofday = f_data_combined_gridded.groupby('local_hour_of_day')[['Average']].median()
f_data_combined_gridded_by_timeofday['Count'] = f_data_combined_gridded.groupby('local_hour_of_day')[['Average']].count()
f_data_combined_gridded_by_timeofday['75th_percentile'] = f_data_combined_gridded.groupby('local_hour_of_day')[['Average']].quantile(0.75)
f_data_combined_gridded_by_timeofday['25th_percentile'] = f_data_combined_gridded.groupby('local_hour_of_day')[['Average']].quantile(0.25)
4g. Time Series — Single Grid Cell
Extract one grid cell and plot its daily PM2.5 trajectory over the full study period. The notebook uses plt.figure() / add_subplot() and labels the variable f_data_combined_gridded_daily_at_Location:
# Select only the data for the grid cell corresponding to the location of interest:
f_data_combined_gridded_daily_at_Location = f_data_combined_gridded_daily.loc[
(lat_cell, lon_cell)
]
# Plot the time series:
o_fig = plt.figure(figsize=(15, 5))
o_ax = o_fig.add_subplot()
o_ax.plot(f_data_combined_gridded_daily_at_Location.index,
f_data_combined_gridded_daily_at_Location['Average'],
color='black', label='Monitor Average')
o_ax.plot(f_data_combined_gridded_daily_at_Location.index,
f_data_combined_gridded_daily_at_Location['MERRA2_CNN'],
color='red', label='MERRA-2 CNN Average')
o_ax.legend()
o_ax.set_ylabel('Daily Average PM$_{2.5}$ [$\mu$g/m$^{3}$]')
o_ax.set_title(f'Location: {s_location}')
Data Flow Summary
| Stage | Input | Output |
|---|
| Domain definition | User parameters | Bounding box + time range variables |
| SatPM download | S3 monthly NetCDF | Subsetted monthly NetCDF files |
| MERRA-2 download | earthaccess / NASA EarthData | Hourly NetCDF under data/MERRA2/ |
| Monitor download | OpenAQ API or SINCA | f_data_monitors_raw DataFrame |
| Grid snapping | f_data_monitors | f_data_combined_gridded (hourly) |
| Daily aggregation | f_data_combined_gridded | f_data_combined_gridded_daily |
| Threshold split | f_data_combined_gridded_daily | f_above, f_below DataFrames |
| Diel cycle | f_data_combined_gridded | f_data_combined_gridded_by_timeofday |
| Time series | f_data_combined_gridded_daily | Single-cell line plot |
After completing the grid-snapping step, save f_data_combined_gridded to Google Drive. It is the most expensive object to reconstruct and serves as the input for steps 4d–4g.
