This notebook-based case study applies the full satellite PM₂.₅ analysis workflow to one of South America’s most dramatic recent air-quality events: the February 2023 Central Chile megafires. Working with the SatPM V6.GL.03 monthly product and the hourly MERRA-2 CNN dataset, you will download data, compare them to ground monitors sourced directly from Chile’s SINCA network, and perform a suite of map, time-series, and statistical analyses across a domain spanning Santiago to Concepción.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.
This notebook runs entirely inside Google Colab and requires (1) a free NASA Earthdata account and (2) a free OpenAQ API key. Both are obtained online in a few minutes.
Open in Google Colab
🇺🇸 Exercise (EN)
English exercise notebook — fill-in-the-blank cells omitted for training
🇺🇸 Complete (EN)
English complete notebook — all cells pre-filled
🇪🇸 Ejercicio (ES)
Cuaderno de ejercicios en español
🇪🇸 Completo (ES)
Cuaderno completo en español
Notebook authors
- Carl Malings — Morgan State University & NASA Goddard Space Flight Center
- Sebastián Diez — Centro de Investigación en Tecnologías para la Sociedad, Universidad del Desarrollo (Chile)
Notebook walkthrough
The notebook opens with context on why this event is scientifically interesting. The February 2023 megafires in Central Chile produced extreme PM₂.₅ episodes visible from both satellites and ground monitors. Because the fires generated data that many monitoring networks initially flagged as not validated, ground-truth retrieval requires care — motivating the SINCA-vs-OpenAQ comparison covered in a later step.
The setup cell installs and imports every dependency, mounts your Google Drive so downloaded data can persist, and prompts for credentials.
# Install & Import the required packages:
!pip install --quiet cartopy earthaccess
import os, getpass
from google.colab import drive
import warnings; warnings.filterwarnings('ignore')
import numpy as np
import pandas as pd
import xarray as xr
import earthaccess, requests, urllib, re, time
from urllib.parse import urljoin
from bs4 import BeautifulSoup
import matplotlib.pyplot as plt
import cartopy
import seaborn as sns
# Connect to Google Drive:
drive.mount('/content/drive')
# Store login information:
earthaccess.login()
s_key_openAQ = getpass.getpass('Enter OpenAQ API key:')
Save a copy of this notebook to your own Google Drive (
File → Save a Copy in Drive) before running it. Create two sub-folders — SatPM_Data and MERRA2_CNN_Data — inside your project folder so downloaded files stay organised.Four helper functions are defined here and reused throughout the notebook. All are documented inline with full docstrings.
F_subset_and_combinexarray objectF_get_OpenAQ_from_API(lat, lon, time)-indexed pandas DataFrameF_plot_mapxarray) or point (pandas) PM₂.₅ data on a Cartopy PlateCarree map with coastlines, borders, and a colour barF_compute_metrics/api/f-compute-metricsdef F_compute_metrics(f_data,
s_time = 'time_month_start',
s_value = 'Average',
s_estimate = 'SatPM',
b_print = True):
d_metrics = {}
f_data['error'] = f_data[s_estimate] - f_data[s_value]
d_metrics['Estimated Mean'] = f_data[s_estimate].mean()
d_metrics['Measured Mean'] = f_data[s_value].mean()
d_metrics['Bias'] = f_data['error'].mean()
d_metrics['Mean Normalized Bias'] = d_metrics['Bias'] / d_metrics['Measured Mean']
d_metrics['RMSE'] = np.sqrt((f_data['error']**2).mean())
d_metrics['Mean Normalized RMSE'] = d_metrics['RMSE'] / d_metrics['Measured Mean']
d_metrics['R^2'] = (
1 - (f_data['error']**2).sum()
/ ((f_data[s_value] - d_metrics['Measured Mean'])**2).sum()
)
# ... spatial and temporal correlation loops omitted for brevity
return d_metrics
Chile’s national air-quality monitoring network SINCA (Sistema de Información Nacional de Calidad del Aire) is the preferred ground-truth source for this case study. A complete adapter —
F_get_SINCA_from_web — is defined here and returns data in exactly the same (lat, lon, time) indexed format as F_get_OpenAQ_from_API, so the rest of the notebook is data-source agnostic.Why not just use OpenAQ for Chile?OpenAQ ingests mostly validated records. Chilean SINCA stations report fire-episode data as preliminary or not validated, so those peak values are largely absent from OpenAQ — precisely the extreme concentrations this case study targets. Additionally, OpenAQ ingests SINCA data as 24-hour rolling averages rather than hourly averages, further masking fire-episode peaks.SINCA itself exposes three parallel series per station —
validated, preliminary, and not_validated — and the parser below consumes all three (preferring validated → preliminary → not_validated). In general, always check the original data source; local researchers and experts who collected the data will know best how to use it appropriately.The scraper navigates SINCA’s public website using
requests and BeautifulSoup, identifies PM₂.₅ stations in the requested regions, downloads hourly CSV data through the SINCA CGI endpoint, converts Chilean local time (America/Santiago) to UTC, and clips results to the bounding box.def F_get_SINCA_from_web(n_lon_min, n_lat_min, n_lon_max, n_lat_max,
t_start, t_end, regions, l_variables=["PM2.5"]):
"""Drop-in replacement for F_get_OpenAQ_from_API, reading from SINCA.
Returns a pandas DataFrame with a 3-level (lat, lon, time) index and a
'PM2.5' column, with time in UTC — exactly the format the rest of the
notebook expects. SINCA serves whole regions, so pass the region IDs that
cover the bounding box; results are then clipped to the box.
"""
frames = []
for region_id in as_list(regions):
df_st = get_sinca_region_stations(region_id)
for _, row in df_st.iterrows():
d = _download_one_station_pm25(row, pd.to_datetime(t_start),
pd.to_datetime(t_end))
if d is not None and len(d):
frames.append(d)
df = pd.concat(frames, ignore_index=True)
# Convert Chilean local time → UTC:
df["time"] = (df["datetime_local"]
.dt.tz_localize("America/Santiago",
ambiguous="NaT", nonexistent="NaT")
.dt.tz_convert("UTC"))
df = df[df["lat"].between(n_lat_min, n_lat_max)
& df["lon"].between(n_lon_min, n_lon_max)]
f_data = (df[["lat", "lon", "time", "PM2.5"]]
.dropna(subset=["PM2.5"])
.groupby(["lat", "lon", "time"]).mean())
return f_data
The SINCA region IDs used in this case study are Roman-numeral region codes:
V (Valparaíso), M (Metropolitana), VI (O’Higgins), VII (Maule), XVI (Ñuble), VIII (Biobío), IX (Araucanía).The spatial domain covers the area from Santiago southward to Concepción, and the time window spans the full austral summer surrounding the February 2023 fires.
# Define the bounding box:
n_lat_min = -38
n_lat_max = -32
n_lon_min = -73.75
n_lon_max = -69.75
# Define the start and end time:
t_start = np.datetime64('2022-12-01 00:00:00')
t_end = np.datetime64('2023-03-31 23:59:59')
# Specify the local time zone:
s_timezome = 'America/Santiago'
# Preview the selected region:
F_plot_map(n_lat_min=n_lat_min, n_lat_max=n_lat_max,
n_lon_min=n_lon_min, n_lon_max=n_lon_max)
All timestamps in the notebook are stored as UTC. The
s_timezome variable is used only when converting UTC timestamps to local time for diel-cycle and daily-average visualisations.Monthly SatPM V6.GL.03 files are fetched directly from the SatPM S3 bucket. The loop constructs the URL for each month in the analysis window, skips files already on disk, and calls
F_subset_and_combine to crop every file to the bounding box before concatenating.s_folder_SatPM = '/content/drive/MyDrive/Colab Notebooks/ARSET/PM25_Training_2026/SatPM_Data'
l_files_SatPM = []
for t_month in np.arange(t_start.astype('datetime64[M]'),
t_end.astype('datetime64[M]') + np.timedelta64(1,'M'),
np.timedelta64(1,'M')):
s_year = str(t_month)[0:4]
s_month = str(t_month)[5:7]
s_url = (f'https://s3.amazonaws.com/satpmdata/V6GL03/FineResolution/GL/'
f'Monthly/{s_year}/V6GL03.CNNPM25.GL.{s_year}{s_month}-{s_year}{s_month}.nc')
s_download_path = os.path.join(s_folder_SatPM, s_url.split('/')[-1])
if not os.path.exists(s_download_path):
urllib.request.urlretrieve(s_url, s_download_path)
l_files_SatPM.append(s_download_path)
def F_time_from_SatPM_url(s_url):
s_YYYYMM = s_url.split('.')[-2].split('-')[0]
return np.datetime64(f'{s_YYYYMM[0:4]}-{s_YYYYMM[4:6]}-01 00:00:00')
a_data_SatPM = F_subset_and_combine(
l_files_SatPM, n_lon_min, n_lat_min, n_lon_max, n_lat_max,
t_start, t_end, F_time_from_file=F_time_from_SatPM_url)
# Filter out negative values:
a_data_SatPM['PM25_filtered'] = a_data_SatPM['PM25'].where(a_data_SatPM['PM25'] > 0)
a_data_SatPM.to_netcdf(os.path.join(s_folder_SatPM, 'SatPM_data_subset.nc'))
After downloading, the cell renders a quick map plot so you can visually verify coverage before proceeding.
Hourly MERRA-2 CNN files are searched and downloaded with the
earthaccess Python library, which handles NASA Earthdata authentication automatically.s_folder_MERRA2_CNN = '/content/drive/MyDrive/Colab Notebooks/ARSET/PM25_Training_2026/MERRA2_CNN_Data'
l_results = earthaccess.search_data(
short_name='MERRA2_CNN_HAQAST_PM25',
bounding_box=(n_lon_min, n_lat_min, n_lon_max, n_lat_max),
temporal=(str(t_start), str(t_end))
)
l_files_MERRA2_CNN = earthaccess.download(l_results, s_folder_MERRA2_CNN)
a_data_MERRA2_CNN = F_subset_and_combine(
l_files_MERRA2_CNN, n_lon_min, n_lat_min, n_lon_max, n_lat_max, t_start, t_end)
# Keep only quality-flagged values (QFLAG >= 3):
a_data_MERRA2_CNN['PM25_filtered'] = (
a_data_MERRA2_CNN['MERRA2_CNN_Surface_PM25']
.where(a_data_MERRA2_CNN['QFLAG'] >= 3))
a_data_MERRA2_CNN.to_netcdf(os.path.join(s_folder_MERRA2_CNN, 'MERRA2_CNN_data_subset.nc'))
MERRA-2 CNN data are available at 0.5° latitude × 0.625° longitude resolution with hourly time steps. The
QFLAG >= 3 filter retains only retrievals the algorithm considers reliable.Both data sources are queried for the same domain and time window, then plotted side-by-side at a single Santiago station to illustrate the validation-status issue.
s_folder_Main = '/content/drive/MyDrive/Colab Notebooks/ARSET/PM25_Training_2026'
# Query OpenAQ (reference monitors only):
f_data_monitors_OpenAQ = F_get_OpenAQ_from_API(
n_lon_min, n_lat_min, n_lon_max, n_lat_max, t_start, t_end)
# Query SINCA directly:
f_data_monitors_SINCA = F_get_SINCA_from_web(
n_lon_min, n_lat_min, n_lon_max, n_lat_max,
t_start, t_end,
regions=['V', 'M', 'VI', 'VII', 'XVI', 'VIII', 'IX'])
The comparison time-series plot at Parque O’Higgins (lat −33.464142, lon −70.660797) shows that OpenAQ data are much smoother (24-hour rolling values) and missing the sharpest peaks from the February fires. SINCA hourly data capture those peaks clearly. After this comparison the notebook sets
f_data_monitors = f_data_monitors_SINCA and filters out physically implausible values (PM2.5 < 0 or > 2000 µg/m³).# Side-by-side time-series comparison at Parque O'Higgins:
o_fig = plt.figure(figsize=(15, 5))
o_ax = o_fig.add_subplot()
o_ax.plot(f_data_monitors_OpenAQ.loc[(-33.464142, -70.660797)].index,
f_data_monitors_OpenAQ.loc[(-33.464142, -70.660797)]['PM2.5'],
color='blue', label="Data from OpenAQ")
o_ax.plot(f_data_monitors_SINCA.loc[(-33.464142, -70.660797)].index,
f_data_monitors_SINCA.loc[(-33.464142, -70.660797)]['PM2.5'],
color='orange', label="Data from SINCA")
o_ax.legend()
o_ax.set_ylabel('PM$_{2.5}$ [µg/m$^{3}$]')
o_ax.set_title("Data comparison at Parque O'Higgins")
The analysis section opens by reloading any subset files cached to Drive (handy when returning to the notebook later), then generates two overlay maps — SatPM and MERRA-2 CNN with SINCA monitor points superimposed as coloured scatter points.
# Overlay monitors on the SatPM map:
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, 50])
Monitor data are aggregated to monthly averages (requiring at least
24 × 20 = 480 valid hourly samples per station-month) and matched to the co-located SatPM grid point. A scatter plot colour-coded by month and F_compute_metrics output quantify agreement.The 480-sample threshold corresponds to 20 full days of valid hourly data in a month. Adjust
n_minimum_samples_per_month if your region has sparser coverage.Direct hourly comparison is possible because MERRA-2 CNN provides hourly output. To enable a spatially fair comparison, monitor stations are first snapped to the nearest MERRA-2 grid cell (0.5° lat × 0.625° lon) and averaged together — a MERRA-2 pixel often covers multiple stations. A grid cell must contain at least 3 valid monitor samples at that hour to be included.
n_minimum_samples_per_grid = 3
# Create a new data frame for helping to compute the grid averages:
f_data_monitors_gridded = f_data_monitors.copy().reset_index()
# Assign each monitor to the nearest MERRA-2 grid centre:
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
# Group by grid cell and hour — compute count, mean, and standard deviation:
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()
)
# Enforce the minimum-samples filter:
f_data_combined_gridded = (
f_data_combined_gridded
.where(f_data_combined_gridded['Count'] >= n_minimum_samples_per_grid)
.dropna()
)
Local time-of-day and local date columns are appended at this stage for the diel-cycle and daily-average analyses that follow.
The 24-hour median and interquartile range are computed for both monitor grids and MERRA-2 CNN to visualise whether the model captures the diurnal PM₂.₅ pattern (e.g., overnight wood-burning peaks in Chilean winter).
Gridded hourly data are further aggregated to local calendar days and a second scatter-plot / metrics block assesses performance at the daily scale — typically showing improved R² relative to hourly comparison because averaging reduces noise.
A configurable daily-average threshold (
n_threshold = 100 µg/m³ for this fire event) splits the dataset into “above threshold” (fire-affected) and “below threshold” (background) periods. Performance metrics are computed separately for each regime to examine whether the satellite products track background conditions better than extreme fire episodes.n_threshold = 100 # µg/m³ daily average
f_data_combined_gridded['Daily_Average_Above_Threshold'] = False
for i_row in range(len(f_data_combined_gridded)):
n_lat = f_data_combined_gridded.index[i_row][0]
n_lon = f_data_combined_gridded.index[i_row][1]
t_time = f_data_combined_gridded.index[i_row][2]
t_day = f_data_combined_gridded.loc[(n_lat, n_lon, t_time), 'local_day']
if f_data_combined_gridded_daily.loc[(n_lat, n_lon, t_day), 'Average'] > n_threshold:
f_data_combined_gridded.loc[(n_lat, n_lon, t_time),
'Daily_Average_Above_Threshold'] = True
# Metrics for non-fire days only:
F_compute_metrics(
f_data_combined_gridded[
f_data_combined_gridded['Daily_Average_Above_Threshold'] == False],
s_time='time', s_value='Average', s_estimate='MERRA2_CNN')
The final cell plots daily-average time series for three pre-defined MERRA-2 grid locations, contrasting monitor observations (black) against MERRA-2 CNN estimates (red). The February 2023 fire signal is clearest at Talcahuano y Concepción, which sits at the epicentre of the fires.
d_location = {
'Santiago': (-33.5, -70.625),
'Viña del Mar': (-33.0, -71.25),
'Talcahuano y Concepción': (-37.0, -73.125),
}
s_location = 'Talcahuano y Concepción'
f_ts = f_data_combined_gridded_daily.loc[d_location[s_location]]
o_fig = plt.figure(figsize=(15, 5))
o_ax = o_fig.add_subplot()
o_ax.plot(f_ts.index, f_ts['Average'], color='black', label='Monitor Average')
o_ax.plot(f_ts.index, f_ts['MERRA2_CNN'], color='red', label='MERRA-2 CNN Average')
o_ax.legend()
o_ax.set_ylabel('Daily Average PM$_{2.5}$ [µg/m$^{3}$]')
o_ax.set_title(f'Location: {s_location}')
Key concepts and cross-references
| Concept | Details |
|---|---|
| Ground-truth network | SINCA — Chile’s national air-quality monitoring network; exposes hourly validated, preliminary, and not-validated PM₂.₅ data |
| Aggregator comparison | OpenAQ — ingests mostly validated records; for fire episodes this means missing not-validated peaks |
| Monthly satellite product | SatPM V6.GL.03 — 0.01° resolution, global monthly PM₂.₅ estimates |
| Hourly satellite product | MERRA-2 CNN (DOI: 10.5067/OCKK5HCFW5N3) — 0.5° × 0.625°, hourly |
| Performance metrics | F_compute_metrics — computes Bias, RMSE, MNB, MNRMSE, R², spatial and temporal correlations |
| MERRA-2 grid snapping | Monitors snapped to nearest 0.5° lat × 0.625° lon grid centre; minimum 3 valid samples required per grid cell per hour |