Skip to content

Instantly share code, notes, and snippets.

@akhuff157

akhuff157/tutorial_goes_abi_pm25.ipynb

Last active March 25, 2025 00:52
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save akhuff157/5b95859cd4223b801f9c7fe4023cbf8b to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save akhuff157/5b95859cd4223b801f9c7fe4023cbf8b to your computer and use it in GitHub Desktop.
Download ZIP
tutorial_goes_abi_pm25.ipynb
Display the source blob
Display the rendered blob
Raw
tutorial_goes_abi_pm25.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"provenance": [],
"authorship_tag": "ABX9TyMjiljM0sx48PppNuIN0fwQ",
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
},
"language_info": {
"name": "python"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/akhuff157/5b95859cd4223b801f9c7fe4023cbf8b/tutorial_goes_abi_pm25.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "markdown",
"source": [
"# **Short Tutorial: How to Open, Read, and Visualize Satellite-Estimated Hourly Surface PM2.5 Data Files**\n",
"\n",
"This tutorial was written in August 2024 by Dr. Amy Huff, IMSG at NOAA/NESDIS/STAR (amy.huff@noaa.gov). Please **credit the NOAA/NESDIS/STAR Aerosols and Atmospheric Composition Science Team** if you use any of the tutorial code in your work.\n",
"\n",
"A description of the algorithm used to generate hourly surface PM2.5 estimates from GOES-East & -West Advanced Baseline Imager (ABI) aerosol optical depth (AOD) is given in [Zhang and Kondragunta, Earth and Space Science, 2021](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2020EA001599)."
],
"metadata": {
"id": "25ks3Xnkivop"
}
},
{
"cell_type": "markdown",
"source": [
"## **About Google Colab**\n",
"\n",
"Google Colab is a hosted Jupyter Notebook service that requires no setup to use and provides free access to computing resources.\n",
"\n",
"Jupyter Notebook is an open-source web application that supports >40 programming languages, including Python. It allows code to be broken into “blocks” that run independently, which makes it ideal for tutorials.\n",
"\n",
"While Colab is free, powerful, and easy to use, it has some limitations that users should keep in mind:\n",
"\n",
"\n",
"1. Colab sessions are temporary. A session will expire after 12 hours of continuous use or after 90 minutes of idle time. *Any output, including files you download or generate, will be lost after the session expires.* Therefore, any files users want to save must be downloaded to the user's local computer or Google Drive account.\n",
"2. It is not possible to run Colab using a virtual or conda configuration environment. Therefore, in general, users must work with the existing Colab configuration.\n",
"\n",
"**The Python code demonstrated in this tutorial is universal.** Specific lines of code or functions will run in any Python IDE (e.g., Spyder, Visual Studio Code), Jupyter Notebook, or the Python interpreter."
],
"metadata": {
"id": "5lhQyrdvkpwK"
}
},
{
"cell_type": "markdown",
"source": [
"### **Import Python modules and packages**\n",
"\n",
"- [datetime](https://docs.python.org/3/library/datetime.html): module to manipulate dates and times\n",
"- [pathlib](https://docs.python.org/3/library/pathlib.html): module to set file system paths corresponding to the user's operating system\n",
"- [io](https://docs.python.org/3/library/io.html): module for working with I/O streams\n",
"- [xarray](https://docs.xarray.dev/en/stable/index.html): library to work with labeled multi-dimensional arrays\n",
"- [requests](https://requests.readthedocs.io/en/latest/): library to send HTTP requests\n",
"- [matplotlib](https://matplotlib.org/stable/index.html): library to make plots\n",
"- [cartopy](https://scitools.org.uk/cartopy/docs/latest/index.html): library to make maps\n",
"\n",
"The Colab configuration does not include the `cartopy` package, so it needs to be installed before it is imported.\n",
"\n",
"**You may ignore any messages about package dependency conflicts; they will not impact the tutorial.**"
],
"metadata": {
"id": "I44pITOYlAsN"
}
},
{
"cell_type": "code",
"source": [
"# Install cartopy package in Google Colab quietly (no progress notifications)\n",
"\n",
"!pip install --quiet cartopy"
],
"metadata": {
"id": "fOMHxdqmkE-_"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"# Import modules and packages\n",
"\n",
"import datetime\n",
"\n",
"from pathlib import Path\n",
"\n",
"from io import BytesIO\n",
"\n",
"import requests\n",
"\n",
"import xarray as xr\n",
"\n",
"import matplotlib as mpl\n",
"from matplotlib import pyplot as plt\n",
"\n",
"from cartopy import crs as ccrs\n",
"import cartopy.feature as cfeature"
],
"metadata": {
"id": "y45jRl3-mdIm"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## **Archive of Satellite-Estimated Hourly Surface PM2.5 Data Files**\n",
"\n",
"A rolling 10-day archive of the satellite-estimated hourly surface PM2.5 data files are available on the [STAR FTP](https://www.star.nesdis.noaa.gov/pub/smcd/hzhang/GOES-16/NRT/CONUS/pm25gwr/). \n",
"\n",
"The files are in [netCDF4 format](https://docs.unidata.ucar.edu/netcdf-c/current/netcdf_data_model.html) and contain latitude, longitude, and PM2.5 variables derived from [GOES-East and GOES-West ABI CONUS sector](https://www.star.nesdis.noaa.gov/atmospheric-composition-training/satellite_data_abi_scan_sectors.php#conus) observations.\n",
"\n",
"Hourly surface PM2.5 estimates correspond to the \"starting hour\" 1-hour average surface PM2.5 concentration, e.g. the 16:00 UTC file corresponds to the average PM2.5 for 16:00-16:59 UTC."
],
"metadata": {
"id": "jdWlQcxEm71Z"
}
},
{
"cell_type": "markdown",
"source": [
"## **Example: Satellite-Estimated Surface PM2.5 for 16 August 2024 at 16:00 UTC**\n",
"\n",
"The Python code below will open the satellite-estimated surface PM2.5 file for 16:00 UTC on 16 August 2024, process the data, and plot the estimated PM2.5 concentrations on a map.\n",
"\n",
"**Ignore the `DownloadWarning` messages**, which are triggered the first time a shapefile is accessed via the `cartopy` [Feature interface](https://scitools.org.uk/cartopy/docs/latest/reference/feature.html).\n",
"\n",
"Notes on the code:\n",
" - The data file contents are accessed using `requests`.\n",
" - The data file is opened remotely (without downloading) as a data stream using the `BytesIO` [operator](https://docs.python.org/3/library/io.html#io.BytesIO) of the `io` module in conjunction with the `xarray.open_dataset()` [function](https://docs.xarray.dev/en/stable/generated/xarray.open_dataset.html).\n",
" - Select a map projection from the `cartopy` [map projection options](https://scitools.org.uk/cartopy/docs/latest/reference/projections.html). I used the [Plate Carree](https://scitools.org.uk/cartopy/docs/latest/reference/projections.html#platecarree) equidistant cylindrical (equirectangular) projection in this example.\n",
" - The argument `crs=ccrs.PlateCarree()` must be used with the `set_extent()` [function](https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html#cartopy.mpl.geoaxes.GeoAxes.set_extent) because the `bounding_box_corners` are entered in geographic coordinates (latitude and longitude).\n",
" - The `matplotlib.pyplot.pcolormesh()` [function](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.pcolormesh.html) is the preferred option for plotting 2D satellite data, such as estimated PM2.5. However, the `X` or `Y` argument arrays (i.e., latitude and longitude) for `pcolormesh()` are not allowed to be masked (have missing values). So the `xarray.DataArray.fillna()` [function](https://docs.xarray.dev/en/stable/generated/xarray.DataArray.fillna.html) is used to fill missing (masked) cells in the latitude and longitude arrays with a fill value of -999.99. This tricks `pcolormesh()` into plotting the estimated PM2.5, ignoring the fill values because they are outside of the valid data range.\n",
" - The argument `transform=ccrs.PlateCarree()` must be used with the `matplotlib.pyplot.pcolormesh()` [function](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.pcolormesh.html) because the plotted data are in geographic coordinates (latitude and longitude).\n",
"\n",
"**Don't forget: any files you want to keep must be downloaded from the Colab instance to Google Drive or your local computer!** Use the controls in the `Files` menu panel on the left side of the Colab window to mount your Google Drive to the Colab instance or to download files to your local computer."
],
"metadata": {
"id": "24DuVT4UpmMn"
}
},
{
"cell_type": "code",
"source": [
"# Access data file from STAR FTP website using requests\n",
"url = 'https://www.star.nesdis.noaa.gov/pub/smcd/hzhang/GOES-16/NRT/CONUS/pm25gwr/20240816/'\n",
"file_name = 'pm25_gwr_aod_exp50_2024081616.nc'\n",
"response = requests.get(url + file_name)\n",
"\n",
"# Open file contents without downloading (automatically closes file when done)\n",
"# Read file as a binary file-like object using BytesIO\n",
"with BytesIO(response.content) as remote_file:\n",
" # Open remote file using xarray\n",
" with xr.open_dataset(remote_file, engine='h5netcdf') as ds:\n",
"\n",
" # Set up figure and map projection with matplotlib & cartopy\n",
" fig = plt.figure(figsize=(8, 10))\n",
" ax = plt.axes(projection=ccrs.PlateCarree())\n",
"\n",
" # Set geographic domain of map [W_lon, E_lon, S_lat, N_lat]\n",
" # °E longitude > 0 > °W longitude, °N latitude > 0 > °S latitude\n",
" ax.set_extent([-132, -60, 14, 52.5], crs=ccrs.PlateCarree())\n",
"\n",
" # Add map borders/shading using cartopy\n",
" ax.add_feature(cfeature.LAND)\n",
" ax.add_feature(cfeature.OCEAN)\n",
" ax.add_feature(cfeature.LAKES)\n",
" ax.add_feature(cfeature.COASTLINE, linewidth=0.5)\n",
" ax.add_feature(cfeature.BORDERS, linewidth=0.5)\n",
" ax.add_feature(cfeature.STATES, linewidth=0.5)\n",
"\n",
" # Create colormap for PM2.5 (based on AQI) using matplotlib\n",
" aqi_cmap = mpl.colors.ListedColormap(['limegreen', 'yellow', 'orange',\n",
" 'red', 'darkorchid', 'firebrick'])\n",
" bounds = [0, 9.1, 35.5, 55.5, 125.5, 225.5, 500]\n",
" cmap_norm = mpl.colors.BoundaryNorm(bounds, aqi_cmap.N)\n",
"\n",
" # Read in latitude, longitude, satellite-estimated PM2.5 data\n",
" lon_ge, lat_ge, pm25sat_ge = (ds.lon_ge.fillna(-999.99), \n",
" ds.lat_ge.fillna(-999.99), ds.pm25sat_ge)\n",
" lon_gw, lat_gw, pm25sat_gw = (ds.lon_gw.fillna(-999.99), \n",
" ds.lat_gw.fillna(-999.99), ds.pm25sat_gw)\n",
"\n",
" # Plot GOES-East data\n",
" plot1 = ax.pcolormesh(lon_ge, lat_ge, pm25sat_ge, cmap=aqi_cmap,\n",
" norm=cmap_norm, transform=ccrs.PlateCarree())\n",
"\n",
" # Plot GOES-West data\n",
" plot2 = ax.pcolormesh(lon_gw, lat_gw, pm25sat_gw, cmap=aqi_cmap,\n",
" norm=cmap_norm, transform=ccrs.PlateCarree())\n",
"\n",
" # Add colorbar\n",
" cb = fig.colorbar(plot1, orientation='horizontal', fraction=0.2, pad=0.02,\n",
" shrink=0.5, ticks=bounds, spacing='uniform', extend='max')\n",
" cb.set_label(label='PM$_{2.5}$ Concentration ($\\mu$g m$\\mathregular{^{-3}}$)',\n",
" size=8, weight='bold')\n",
" cb.ax.set_xticklabels(['0', '9.1', '35.5', '55.5', '125.5', '225.5', '500'])\n",
" cb.ax.tick_params(labelsize=8)\n",
"\n",
" # Show plot\n",
" plt.show()\n",
"\n",
" # Save image file (find in \"Files\" folder on Google Colab)\n",
" save_name = 'pm25_hourly_satellite_estimated_' + file_name.split('_')[4][:10]\n",
" fig.savefig(save_name, dpi=300, bbox_inches='tight')\n",
"\n",
" # Close plot\n",
" plt.close()"
],
"metadata": {
"id": "emFcUwdYpJnY"
},
"execution_count": null,
"outputs": []
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment