{ "cells": [ { "cell_type": "markdown", "id": "60bcdae2", "metadata": {}, "source": [ "## Processing PACE and EMIT Data\n", "\n", "Authors: Skye Caplan1 and Erik Bolch2\n", "\n", "1. Science Systems and Applications, Inc., NASA Goddard Space Flight Center.\n", "2. KBR Inc., contractor to the USGS Earth Observation Research and Science Center, NASA Land Processes Distributed Active Archive Center.\n", "\n", "\n", "Last updated: 06/24/2026\n", "\n", "[edl]: https://urs.earthdata.nasa.gov/\n", "[oci-data-access]: https://oceancolor.gsfc.nasa.gov/resources/docs/tutorials/notebooks/oci_data_access/\n", "[emit-data-access]: https://nasa.github.io/VITALS/python/Exploring_EMIT_L2A_RFL.html\n", "[help-hub]: https://nasa.github.io/oceandata-notebooks/\n", "[vitals]: https://nasa.github.io/VITALS/\n", "\n", "
\n", "\n", "An [Earthdata Login][edl] account is required to access data from the NASA Earthdata system, including NASA PACE and EMIT data.\n", "\n", "
\n", "\n", "## Summary\n", "NASA's Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) satellite mission and Earth Surface Mineral Dust Source Investigation (EMIT) sensor aboard the International Space Station (ISS) both make hyperspectral observations of the Earth's surface. Taking into account their different spectral, spatial, and temporal coverage, these instruments' complementary capabilities provide a wealth of new opportunities when used in concert. In this series of tutorials, we will demonstrate how to find, prepare, and use data from PACE and EMIT together, illustrating both missions' strengths and how they can help fill in each other's gaps. \n", "\n", "This is the second tutorial in the PACE-EMIT series. This notebook will demonstrate how to preprocess the reflectance data downloaded from the first tutorial, `01_Colocate_PACE_EMIT_data`, so that they can be used together to investigate science questions. We will use functions from the `emit_tools` and `pace_tools` modules to orthorectify, geolocate, and merge EMIT and PACE Ocean Color Instrument (OCI) data where necessary. Those modules will then also be used to place OCI and EMIT granules on an equivalent, stackable grid. We can then compare the data from the two sensors. These techniques will also be integral to our case study analysis in the third tutorial, `03_Wildfires_with_PACE_EMIT`.\n", "\n", "## Background\n", "\n", "The goal of this tutorial series is to look at how EMIT and PACE OCI make complementary observations of the same location. Before getting to the science in notebook 3, the datasets require some preprocessing to help make working with them more intuitive. This is especially true considering we have data from two different instruments with unique orbit and observation characteristics. \n", "\n", "Recall that the reflectances downloaded in the first notebook are in raw Level-2 (L2) data format, meaning both datasets are still in swath resolution. Additionally, neither dataset is geolocated or masked for quality issues. Referring to the instrument summary table below, we also know that EMIT and PACE OCI are not at the same spatial resolution, which can make comparing spectra more complicated. Using pre-written tools for each instrument, we can geolocate and grid the data so that each pixel is associated with a physical latitude/longitude coordinate. To deal with the different pixel sizes, we can either upscale PACE OCI to the finer EMIT grid, or downscale EMIT to the coarser PACE OCI grid, ensuring the data are stackable and directly comparable.\n", "\n", "| Resolutions | EMIT | PACE OCI |\n", "| ----------- | ------- | -------- |\n", "| Spatial | 60 m | 1.2 km (at nadir)|\n", "| Temporal | Variable| 1 - 2 days|\n", "| Spectral Range | 380 - 2500 nm | 350 - 895, with 7 multispectral SWIR bands from 935-2258 nm|\n", "| Spectral Sampling | 7.5 nm | 2.5 nm, with 1.25 nm in some spectral regions |\n", "\n", "## Learning Objectives\n", "\n", "At the end of this notebook, you will know how to:\n", "- Use prebuilt tools to orthorectify, merge, and geolocate EMIT and PACE OCI data\n", "- Upscale PACE OCI to EMIT resolution\n", "- Compare spectral signatures over a region of interest (ROI)\n", "\n", "## Contents\n", "1. [Setup](#1-setup)\n", "2. [Preprocessing Instrument Data](#2-preprocessing-instrument-data) \n", "3. [Grid Alignment](#3-grid-alignment) \n", "4. [Spectral Comparisons](#4-spectral-comparisons) \n", "\n" ] }, { "cell_type": "markdown", "id": "001cffd9", "metadata": {}, "source": [ "## 1. Setup \n", "\n", "Begin by importing all of the packages used in this notebook, and defining some functions we will use in the notebook." ] }, { "cell_type": "code", "execution_count": null, "id": "bed28531", "metadata": {}, "outputs": [], "source": [ "import cartopy\n", "import skimage\n", "import numpy as np \n", "import xarray as xr \n", "import hvplot.xarray\n", "import holoviews as hv\n", "import rioxarray as rio\n", "import geopandas as gpd\n", "import cartopy.crs as ccrs \n", "import matplotlib.pyplot as plt\n", "import earthaccess \n", "from pathlib import Path\n", "NB_DIR = Path.cwd()\n", "\n", "from vitals import pace_tools as pace\n", "from vitals import emit_tools as emit\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "\n", "def plot_features(ax, gridline_alpha=1):\n", " \"\"\"\n", " Function to add gridlines and cartopy features to a figure.\n", " Args:\n", " ax - matplotlib Axes object on which to plot the features\n", " gridline_alpha - fractional opacity of gridlines on the plot. Will also \n", " make sure x, y axes are labelled with lon/lats, respectively. 0 will \n", " include axes labels but not plot gridlines\n", " \"\"\"\n", " for a in range(len(ax)):\n", " ax[a].gridlines(draw_labels={\"left\": \"y\", \"bottom\": \"x\"}, linewidth=0.25, alpha=gridline_alpha)\n", " ax[a].coastlines(linewidth=0.5)\n", " ax[a].add_feature(cartopy.feature.OCEAN, edgecolor=\"w\", linewidth=0.01)\n", " ax[a].add_feature(cartopy.feature.LAND, edgecolor=\"w\", linewidth=0.01)\n", " ax[a].add_feature(cartopy.feature.LAKES, edgecolor=\"w\", linewidth=0.01)\n", "\n", "def _file_stem(f):\n", " \"\"\"Get filename stem from a Path or fsspec file object.\"\"\"\n", " if hasattr(f, 'path'):\n", " return Path(f.path.split('/')[-1]).stem\n", " elif hasattr(f, 'info'):\n", " return Path(f.info()['name'].split('/')[-1]).stem\n", " return f.stem" ] }, { "cell_type": "markdown", "id": "8299aaec", "metadata": {}, "source": [ "In the last notebook, we found and optionally downloaded EMIT and PACE data corresponding to pre-fire and post-fire days for the Grampians National Park fire in Victoria, Australia, as well as several PACE OCI granules between these days. Since we'll be merging some adjacent scenes from EMIT and regridding data from the instruments to align, we need a standard extent for cropping and gridding. Open the `grampians_bbox.geojson` which is a bounding box generated from the grampians shapefile used in the first notebook." ] }, { "cell_type": "code", "execution_count": null, "id": "94db18d0", "metadata": {}, "outputs": [], "source": [ "# Set up extents with GDF\n", "gdf_bbox = gpd.read_file(NB_DIR / \"data/grampians_bbox.geojson\")" ] }, { "cell_type": "markdown", "id": "ab388375", "metadata": {}, "source": [ "This notebook supports a download based or streaming workflow. The download workflow is usually the best pathway if outside the cloud and working on a local computer.The streaming workflow is typically recommended if working in the AWS us-west-2 cloud, where connection speeds are very fast.\n", "\n", "The cell below provides the URLs for scenes used in this notebook, and logs in using your Earthdata Login, in case you need to re-download or want to stream the data. " ] }, { "cell_type": "code", "execution_count": null, "id": "a9137370", "metadata": {}, "outputs": [], "source": [ "# Define URLs for download\n", "emit_urls = [\"https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20241114T035759_2431902_009/EMIT_L2A_RFL_001_20241114T035759_2431902_009.nc\",\n", " \"https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20241114T035759_2431902_009/EMIT_L2A_MASK_001_20241114T035759_2431902_009.nc\",\n", " \"https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20241114T035747_2431902_008/EMIT_L2A_RFL_001_20241114T035747_2431902_008.nc\",\n", " \"https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20241114T035747_2431902_008/EMIT_L2A_MASK_001_20241114T035747_2431902_008.nc\",\n", " \"https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20250319T022818_2507801_008/EMIT_L2A_RFL_001_20250319T022818_2507801_008.nc\",\n", " \"https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20250319T022818_2507801_008/EMIT_L2A_MASK_001_20250319T022818_2507801_008.nc\",\n", " \"https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20250319T022830_2507801_009/EMIT_L2A_RFL_001_20250319T022830_2507801_009.nc\",\n", " \"https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20250319T022830_2507801_009/EMIT_L2A_MASK_001_20250319T022830_2507801_009.nc\"\n", " ]\n", "# NOTE: We can use a handful of these or only host the subsets to reduce data size if necessary (sc)\n", "pace_urls = ['https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20241114T035224.L2.SFREFL.V3_1.nc',\n", " 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20241215T033822.L2.SFREFL.V3_1.nc',\n", " 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250101T033916.L2.SFREFL.V3_1.nc',\n", " 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250118T033938.L2.SFREFL.V3_1.nc',\n", " 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250122T041511.L2.SFREFL.V3_1.nc',\n", " 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250205T040925.L2.SFREFL.V3_1.nc',\n", " 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250211T042542.L2.SFREFL.V3_1.nc',\n", " 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250212T032254.L2.SFREFL.V3_1.nc',\n", " 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250224T035613.L2.SFREFL.V3_1.nc',\n", " 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250301T033714.L2.SFREFL.V3_1.nc',\n", " 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250304T034526.L2.SFREFL.V3_1.nc',\n", " 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250305T042058.L2.SFREFL.V3_1.nc',\n", " 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250307T035339.L2.SFREFL.V3_1.nc',\n", " 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250313T041003.L2.SFREFL.V3_1.nc',\n", " 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250319T042120.L2.SFREFL.V3_1.nc',\n", " ]" ] }, { "cell_type": "markdown", "id": "298403ae", "metadata": {}, "source": [ "### Download Workflow\n", "\n", "If you did not download the data in `01_Colocate_PACE_EMIT_Data.ipynb`, you can uncomment the two cells below. If you've already downloaded the data in the first notebook, uncomment only the second cell below to set the filepaths for the scenes you've already downloaded. " ] }, { "cell_type": "code", "execution_count": null, "id": "0d36e461", "metadata": {}, "outputs": [], "source": [ "# Download Data\n", "earthaccess.login()\n", "emit_paths = earthaccess.download(emit_urls, \"data/\")\n", "oci_paths = earthaccess.download(pace_urls, \"data/\")" ] }, { "cell_type": "code", "execution_count": null, "id": "7dc6864f", "metadata": {}, "outputs": [], "source": [ "# Download Workflow\n", "\n", "emit_paths = sorted((NB_DIR / \"data\").glob(\"EMIT_L2A*.nc\"))\n", "oci_paths_all = sorted((NB_DIR / \"data\").glob(\"PACE_OCI*.nc\"))\n", "oci_paths = [oci_paths_all[0], oci_paths_all[-1]]\n", "oci_ts_paths = oci_paths_all[1:-1]" ] }, { "cell_type": "markdown", "id": "c4c7fa4d", "metadata": {}, "source": [ "### Cloud Workflow\n", "\n", "To stream the data, uncomment out the following cell to construct S3 URLs from the existing URL lists and open them using `earthaccess.open`." ] }, { "cell_type": "code", "execution_count": null, "id": "36e0e49a", "metadata": {}, "outputs": [], "source": [ "# # Cloud Workflow\n", "# earthaccess.login()\n", "\n", "# # Sort the EMIT URLs since they contain both RFL and MASK files then open with earthaccess\n", "# emit_paths = earthaccess.open(\n", "# sorted(\n", "# [u.replace('https://data.lpdaac.earthdatacloud.nasa.gov/', 's3://') for u in emit_urls],\n", "# key=lambda u: u.split('/')[-1]\n", "# ), provider = \"LPCLOUD\", open_kwargs={\"max_concurrency\": 20}\n", "# )\n", "\n", "# # Repeat the Same for the PACE Scenes\n", "# oci_paths = earthaccess.open(\n", "# sorted(\n", "# [u.replace('https://obdaac-tea.earthdatacloud.nasa.gov/', 's3://') for u in [pace_urls[0], pace_urls[-1]]],\n", "# key=lambda u: u.split('/')[-1]\n", "# ), provider= \"OB_CLOUD\", open_kwargs={\"max_concurrency\": 20}\n", "# )\n", "\n", "# oci_ts_paths = earthaccess.open(\n", "# sorted(\n", "# [u.replace('https://obdaac-tea.earthdatacloud.nasa.gov/', 's3://') for u in pace_urls[1:-1]],\n", "# key=lambda u: u.split('/')[-1]\n", "# ), provider=\"OB_CLOUD\", open_kwargs={\"max_concurrency\": 20}\n", "# )" ] }, { "cell_type": "markdown", "id": "1d5c76b9", "metadata": {}, "source": [ "## 2. Preprocessing Instrument Data \n", "\n", "Before we can use these data for science, there are a couple steps we should take to preprocess the files so that they are easier to handle in the long run. Recall from the previous tutorial that we are using Level-2 (L2) data. Raw L2 data are still in the instrument swath, or in other words, not geolocated or on any type of uniform grid. As such, our pixels are yet associated with any coordinates, and can represent different sizes on the ground. This is especially true in the case of PACE OCI, since its swath is very wide, but applies to EMIT as well. The first notebook also showed that our region of interest (ROI) spans two EMIT granules, but is only a small part of any PACE OCI scenes. To simplify our workflow and minimize the size of our datasets, we can use tools written specifically for managing these instruments' data. We'll preprocess the data from each mission separately to work with their distinct characteristics, demonstrating a few different ways of handling data.\n", "\n", "### EMIT Data\n", "\n", "For EMIT data, we want to orthorectify, geolocate, and merge the pre-fire granules and post-fire granules, respectively, so that we have one file for each date covering the full ROI. The cell below contains a function which combines all of these steps into one by calling functions from the `emit_tools` module. More information on the tools below can be found in the [EMIT-Data-Resources repository](https://github.com/nasa/EMIT-Data-Resources).\n", "\n", "Briefly, the `process_and_merge_emit` function calls on `emit_xarray` to open both the reflectance and mask files, and uses the mask to remove any unwanted pixels (clouds, bad data bands, etc.). Then, the `spatial_subset` function is called to clip the data only to the ROI, in this case a bounding box around Grampians National Park. Finally, both granules are merged so the ROI is contained in one file, and that file is saved in an output directory of the user's choice. " ] }, { "cell_type": "code", "execution_count": null, "id": "22bad5cd", "metadata": {}, "outputs": [], "source": [ "def process_and_merge_emit(rfl_paths, mask_paths, polygon_path, mask_bands, outdir):\n", " \"\"\"Preprocess a list of EMIT RFL+MASK pairs and write a single merged NetCDF.\n", " Args:\n", " rfl_paths - list of Paths to EMIT reflectance files to process\n", " mask_paths - list of Paths to EMIT mask files corresponding the reflectance\n", " paths to use in processing\n", " polygon_path - Path to polygon to use for subsettnig\n", " mask_bands - list of which EMIT mask bands to use\n", " outdir - Path to save merged output to\n", " Returns:\n", " out_path - Path to merged output\n", " \"\"\"\n", " outdir.mkdir(parents=True, exist_ok=True)\n", " out_path = outdir / f\"{_file_stem(rfl_paths[0])}_ortho_merged.nc\"\n", " if out_path.exists():\n", " print(\"Output already exists!\")\n", " return out_path\n", " gdf = gpd.read_file(polygon_path)\n", "\n", " scenes = []\n", " for rfl_path, mask_path in zip(rfl_paths, mask_paths):\n", " rfl_ds = emit.emit_xarray(rfl_path)\n", " rfl = emit.spatial_subset(rfl_ds, gdf)\n", " del rfl_ds\n", " mask_ds = emit.emit_xarray(mask_path)\n", " mask = (emit.spatial_subset(mask_ds, gdf)[\"mask\"]\n", " .sel(mask_bands=mask_bands).any(\"mask_bands\").astype(int))\n", " del mask_ds\n", " rfl[\"reflectance\"].values[mask.values == 1, :] = np.nan\n", " scenes.append(emit.ortho_xr(rfl))\n", " del rfl\n", "\n", " merged = emit.merge_emit(scenes, gdf)\n", " del scenes\n", " \n", " merged.to_netcdf(out_path)\n", " return out_path" ] }, { "cell_type": "markdown", "id": "3585250e", "metadata": {}, "source": [ "We set up the arguments for the `process_and_merge_emit` function in the cell below. Let's also print out the list of sorted EMIT files so we know which indices to use for the function call." ] }, { "cell_type": "code", "execution_count": null, "id": "d736ec97", "metadata": {}, "outputs": [], "source": [ "OUTDIR = NB_DIR / \"data\" / \"outputs\"\n", "gdf_path = NB_DIR / \"data/grampians_bbox.geojson\"\n", "mask_bands = ['Cloud flag','Cirrus flag','Water flag','Spacecraft Flag']\n", "\n", "emit_paths" ] }, { "cell_type": "markdown", "id": "d23ab108", "metadata": {}, "source": [ "We can then call the function with the correct filepaths. Note that this cell should take about a minute to run fully. The resulting files will be merged, masked, and subsetted to the ROI." ] }, { "cell_type": "code", "execution_count": null, "id": "99ec74c7", "metadata": {}, "outputs": [], "source": [ "prefire_emit_path = process_and_merge_emit(rfl_paths=emit_paths[4:6],\n", " mask_paths=emit_paths[0:2],\n", " polygon_path=gdf_path,\n", " mask_bands=mask_bands,\n", " outdir=OUTDIR)\n", "postfire_emit_path = process_and_merge_emit(rfl_paths=emit_paths[6:8],\n", " mask_paths=emit_paths[2:4],\n", " polygon_path=gdf_path,\n", " mask_bands=mask_bands,\n", " outdir=OUTDIR)" ] }, { "cell_type": "markdown", "id": "d7f91258", "metadata": {}, "source": [ "Once processing is complete, we can take a look at the contents of the merged files, and plot the data as well." ] }, { "cell_type": "code", "execution_count": null, "id": "c6d37d45", "metadata": {}, "outputs": [], "source": [ "prefire_emit_ds = xr.open_dataset(prefire_emit_path)\n", "postfire_emit_ds = xr.open_dataset(postfire_emit_path)\n", "\n", "prefire_emit_ds" ] }, { "cell_type": "code", "execution_count": null, "id": "7d9d85c6", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1,2, figsize=(9,5),sharey=True, \n", " subplot_kw={\"projection\": ccrs.PlateCarree()})\n", "fig.tight_layout(pad=0.5)\n", "im = prefire_emit_ds[\"reflectance\"].sel(wavelengths=860, \n", " method=\"nearest\").plot(ax=ax[0],vmin=0, vmax=0.6,\n", " add_colorbar=False)\n", "postfire_emit_ds[\"reflectance\"].sel(wavelengths=860, \n", " method=\"nearest\").plot(ax=ax[1],vmin=0, vmax=0.6, \n", " add_colorbar=False)\n", "\n", "plot_features(ax.flatten(), gridline_alpha=0)\n", "bbox = [141.9385, -37.74163, 142.77324, -36.79961]\n", "\n", "cax = fig.add_axes([1.0025, 0.03, 0.03, 0.94]) \n", "fig.colorbar(im, cax=cax, label=\"Surface Reflectance\")\n", "\n", "ax[0].set_title(\"Pre-Fire EMIT Reflectance at 860 nm\")\n", "ax[1].set_title(\"Post-Fire EMIT Reflectance at 860 nm\") \n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "e1961308", "metadata": {}, "source": [ "We have successfully clipped our EMIT data, and the data represent physical locations with coordinates. Let's do the same for the PACE OCI reflectances.\n", "\n", "### PACE OCI Data\n", "\n", "For PACE data, we'll need to subset, mask, and geolocate each granule. A merge is not needed since the full ROI is contained in each scene. The `pace_tools` module provides a suite of functions for working with PACE OCI data, similar to the `emit_tools` module, which the `process_pace` function calls below to accomplish the preprocessing. \n", "\n", "Briefly, the function uses the `open_l2` and `mask_ds` functions to remove pixels containing clouds or water. Then, the `grid_data` function is called with a user-defined resolution to grid the data. We will use a resolution of 0.015 degrees (~1.67 km at the equator) to account for OCI's wide swath. Once the data is gridded, `xarray`'s `.sel()` method can be used to subset the data based on the coordinates of a bounding box. More information on these tools can be found on the [Help Hub](https://nasa.github.io/oceandata-notebooks/). \n", "\n", "Note that in the following cells, we'll process both the pre-fire and post-fire data as well as the time series data across that time period. These cells may take some time to complete due to the size of the files. This cell may take ~10 minutes depending on your connection speed. Occasionally, when working in the cloud, these cells can hang and run without producing any outputs. Stop the kernel with the stop button and re-run the cell. It should pick up from the last file processed." ] }, { "cell_type": "code", "execution_count": null, "id": "a182faa1", "metadata": {}, "outputs": [], "source": [ "def process_pace(pace_path, resolution, bbox, outdir):\n", " out_path = outdir / f\"{_file_stem(pace_path)}_ortho_subset.nc\"\n", " if out_path.exists():\n", " print(\"Output already exists!\")\n", " return out_path\n", " # Load the full data here to avoid an issue with fsspec/datatree\n", " ds = pace.open_l2(pace_path).load()\n", " # Close the earthaccess connection and dataset connection to the source file after the data is loaded into memory\n", " # this seems to improve reliability of the processing loops in the following cells when working in the cloud.\n", " if hasattr(pace_path, 'close'):\n", " pace_path.close()\n", " ds.close()\n", " # Mask Data\n", " ds = pace.mask_ds(ds, flag=[\"CLDICE\",\"LAND\"], reverse=[False,True])\n", " # Project onto a rectilinear grid\n", " ds = pace.grid_data(ds, resolution=resolution)\n", " # Slice to area of interest\n", " ds_sub = ds.sel({\"latitude\":slice(bbox[3],bbox[1]), \"longitude\":slice(bbox[0],bbox[2])})\n", " del ds\n", " # Write Output\n", " ds_sub.to_netcdf(out_path)\n", " del ds_sub\n", " return out_path\n", "\n", "resolution = (0.015,0.015)" ] }, { "cell_type": "code", "execution_count": null, "id": "de3218f0", "metadata": {}, "outputs": [], "source": [ "pace_subsets = [process_pace(oci_path, resolution, gdf_bbox.total_bounds, OUTDIR) for oci_path in oci_paths]" ] }, { "cell_type": "code", "execution_count": null, "id": "bb7c57db", "metadata": {}, "outputs": [], "source": [ "pace_ts_subsets = [process_pace(oci_path , resolution, gdf_bbox.total_bounds, OUTDIR) for oci_path in oci_ts_paths]" ] }, { "cell_type": "markdown", "id": "4bd0d941", "metadata": {}, "source": [ "Let's plot the output pre- and post-fire OCI granules to see the processed data:" ] }, { "cell_type": "code", "execution_count": null, "id": "5ad6c4ed", "metadata": {}, "outputs": [], "source": [ "prefire_oci_ds = xr.open_dataset(pace_subsets[0])\n", "postfire_oci_ds = xr.open_dataset(pace_subsets[1])\n", "\n", "fig, ax = plt.subplots(1,2, figsize=(9,5),sharey=True, \n", " subplot_kw={\"projection\": ccrs.PlateCarree()})\n", "fig.tight_layout(pad=0.5)\n", "prefire_oci_ds[\"rhos\"].sel(wavelength_3d=860, \n", " method=\"nearest\").plot(ax=ax[0],vmin=0, vmax=0.6, add_colorbar=False)\n", "im = postfire_oci_ds[\"rhos\"].sel(wavelength_3d=860, \n", " method=\"nearest\").plot(ax=ax[1],vmin=0, vmax=0.6, add_colorbar=False)\n", "\n", "plot_features(ax.flatten(), gridline_alpha=0)\n", "\n", "cax = fig.add_axes([1.0025, 0.03, 0.03, 0.94])\n", "fig.colorbar(im, cax=cax, label='Surface Reflectance')\n", "\n", "ax[0].set_title(\"Pre-Fire PACE Reflectance at 860 nm\")\n", "ax[1].set_title(\"Post-Fire PACE Reflectance at 860 nm\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "311ed43f", "metadata": {}, "source": [ "The preprocessing for PACE data has worked as well! A quick glance at the above OCI plots compared to those from EMIT illustrates the difference in spatial scale covered by the two missions, but the overall reflectances at 860 nm look similar. \n", "\n", "We can compare full spectra from both instruments for a more detailed view of the data. This could be done by using both datasets at their nominal resolutions, but some analyses may benefit from or even require a PACE dataset with the same resolution as the EMIT dataset (or vice versa). For example, upscaling methods are often employed to prepare a consistent grid for analysis when the goal is to preserve fine resolution details from the higher resolution dataset. Coarsening, on the other hand, can be useful to reduce noise in high-resolution imagery, or to reduce file sizes for large scale regional and global studies. \n", "\n", "Upscaling the OCI data would be useful for the next tutorial in this series since we want to retain features from the EMIT scene while having both datasets on the same grid. We'll use the `grid_data` function from the `pace_tools` module to regrid our OCI data to a finer resolution. \n", "\n", "## 3. Grid Alignment \n", "\n", "First, let's get an idea of what the original EMIT and PACE grids look like overlaid on one another." ] }, { "cell_type": "code", "execution_count": null, "id": "c66fe267", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(13,5), subplot_kw=dict(projection=ccrs.PlateCarree()))\n", "prefire_oci_ds[\"rhos\"].sel({\"wavelength_3d\":860}, method=\"nearest\").plot(vmin=-1, vmax=0, cmap=\"Greys_r\", transform=ccrs.PlateCarree(),linewidth=2, edgecolors=\"red\", add_colorbar=False)\n", "prefire_emit_ds[\"reflectance\"].sel({\"wavelengths\":860}, method=\"nearest\").plot(vmin=0, vmax=1, cmap=\"viridis\", transform=ccrs.PlateCarree(), linewidth=0.05, edgecolors=\"k\", add_colorbar=False, alpha=1)\n", "ax.set_extent([142.615, 142.64, -37.625, -37.60])\n", "ax.gridlines(draw_labels={\"left\": \"y\", \"bottom\": \"x\"}, alpha=0.0)\n", "ax.set_title(\"Pixel Sizes\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "079dd393", "metadata": {}, "source": [ "The discrepancy in pixel size here can be fixed by regridding OCI data with EMIT's specification. Thanks to the processing steps we took in Section 2, we have access to the EMIT grid parameters in the form of its affine transform. We can feed that tranform into the `grid_data` function to coerce the PACE data into that shape. \n", "\n", "
\n", "\n", "It is important to note that sharpening the data does not increase the information content of the data.\n", "\n", "
\n", "\n", "In other words, this approach will be using existing values to fill a larger quantity of pixels covering the same spatial extent, enabling us to stack granules and compare spectra, do band math, and more." ] }, { "cell_type": "code", "execution_count": null, "id": "e9865e60", "metadata": {}, "outputs": [], "source": [ "prefire_oci_finer = pace.grid_data(src=prefire_oci_ds, \n", " dst_transform=prefire_emit_ds.rio.transform())\n", "postfire_oci_finer = pace.grid_data(src=postfire_oci_ds, \n", " dst_transform=postfire_emit_ds.rio.transform())" ] }, { "cell_type": "markdown", "id": "94420cdc", "metadata": {}, "source": [ "After the regridding is done, we can plot the pixel boundaries again to see if they align as we intend them to. " ] }, { "cell_type": "code", "execution_count": null, "id": "296c7668", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(13,5), subplot_kw=dict(projection=ccrs.PlateCarree()))\n", "prefire_oci_finer[\"rhos\"].sel({\"wavelength_3d\":860}, \n", " method=\"nearest\").plot(vmin=-1, vmax=0, cmap=\"Greys_r\", \n", " transform=ccrs.PlateCarree(), \n", " linewidth=2, edgecolors=\"red\", \n", " add_colorbar=False)\n", "prefire_emit_ds[\"reflectance\"].sel({\"wavelengths\":860}, \n", " method=\"nearest\").plot(vmin=0, vmax=1, \n", " transform=ccrs.PlateCarree(), \n", " linewidth=0.05, edgecolors=\"k\", \n", " add_colorbar=False)\n", "\n", "ax.set_extent([142.615, 142.64, -37.625, -37.60])\n", "ax.gridlines(draw_labels={\"left\": \"y\", \"bottom\": \"x\"}, alpha=0.0)\n", "ax.set_title(\"Pixel Alignment\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "24325df6", "metadata": {}, "source": [ "With the plot above, we can see that the pixels now align. Take a look at the shape of each dataset, as well - their x and y dimensions should match in size. \n", "\n", "## 4. Spectral Comparisons \n", "\n", "With our granules now preprocessed and on the same grid, we can look at what some of the full spectral information looks like for our pre- and post-fire datasets. First, to simplify the plotting, we'll make sure our variable names are the same between each dataset by renaming the PACE variables. We'll also change EMIT's `nodata` values to `np.nan` instead of -9999, which would affect the color bar. " ] }, { "cell_type": "code", "execution_count": null, "id": "67fe955d", "metadata": {}, "outputs": [], "source": [ "# Band Masking for EMIT\n", "prefire_emit_ds['reflectance'].data[:,:,prefire_emit_ds['good_wavelengths'].data==0] = np.nan\n", "postfire_emit_ds['reflectance'].data[:,:,postfire_emit_ds['good_wavelengths'].data==0] = np.nan\n", "\n", "# Rename PACE Vars\n", "postfire_oci_finer = postfire_oci_finer.rename({\n", " \"rhos\": \"reflectance\",\n", " \"wavelength_3d\": \"wavelengths\",\n", "})\n", "prefire_oci_finer = prefire_oci_finer.rename({\n", " \"rhos\": \"reflectance\",\n", " \"wavelength_3d\": \"wavelengths\",\n", "})" ] }, { "cell_type": "markdown", "id": "7beebc47", "metadata": {}, "source": [ "We'll then use `hvplot` to create a hoverplot displaying the spectral data. The cells below create true color images from the post-fire OCI and EMIT datasets, which you can then hover over to show what the spectra from each time point and sensor look like. " ] }, { "cell_type": "code", "execution_count": null, "id": "4ac323bf", "metadata": {}, "outputs": [], "source": [ "def scale_for_rgb(d_in):\n", " d_in = np.where(np.isnan(d_in), -9999, d_in)\n", " d = np.where(d_in < -9990, np.nan, d_in)\n", " d_out = np.zeros_like(d)\n", " for i in np.arange(d.shape[-1]):\n", " mi, ma = np.nanpercentile(d[:, :, i], [.5, 99.5])\n", " di = skimage.exposure.rescale_intensity(d_in[:, :, i], in_range=(mi, ma))\n", " if ~np.all(np.isnan(di)):\n", " di = skimage.exposure.equalize_hist(di)\n", " d_out[:, :, i] = skimage.exposure.rescale_intensity(\n", " di, in_range=(di.min(), di.max())\n", " )\n", " return np.where(d_in > -9990, d_out, -9999)\n", "\n", "rgb_ds_emit = postfire_emit_ds['reflectance'].sel(wavelengths=[650,560,470], method='nearest')\n", "rgb_ds_oci = postfire_oci_finer['reflectance'].sel(wavelengths=[650,560,470], method='nearest')\n", "rgb_ds_emit.data = scale_for_rgb(rgb_ds_emit.data)\n", "rgb_ds_oci.data = scale_for_rgb(rgb_ds_oci.data)" ] }, { "cell_type": "code", "execution_count": null, "id": "397726eb", "metadata": {}, "outputs": [], "source": [ "rows, cols = len(rgb_ds_emit.latitude), len(rgb_ds_emit.longitude)\n", "scene_center = rgb_ds_emit.latitude.values[int(rows/2)],rgb_ds_emit.longitude.values[int(cols/2)]\n", "\n", "datasets = {\n", " \"EMIT pre-fire\": prefire_emit_ds,\n", " \"PACE pre-fire\": prefire_oci_finer,\n", " \"EMIT post-fire\": postfire_emit_ds,\n", " \"PACE post-fire\": postfire_oci_finer,\n", "}\n", "colors = {\n", " \"EMIT pre-fire\": \"lightgreen\", \n", " \"PACE pre-fire\": \"forestgreen\",\n", " \"EMIT post-fire\":\"darkorange\", \n", " \"PACE post-fire\":\"firebrick\",\n", "}\n", "\n", "line_styles = {\n", " \"EMIT pre-fire\": \"solid\",\n", " \"EMIT post-fire\": \"solid\",\n", " \"PACE pre-fire\": \"dashed\",\n", " \"PACE post-fire\": \"dashed\",\n", "}\n", "\n", "# Set Start Position\n", "x_start = scene_center[1]\n", "y_start = scene_center[0]\n", "\n", "# Plot RGB\n", "rgb_map_emit = rgb_ds_emit.hvplot.rgb(\n", " x=\"longitude\", y=\"latitude\", bands=\"wavelengths\",\n", " crs=ccrs.PlateCarree(),\n", " title=\"Stretched RGB EMIT Image\",\n", " frame_width=240, frame_height=480\n", ")\n", "\n", "# Plot RGB\n", "rgb_map_oci = rgb_ds_oci.hvplot.rgb(\n", " x=\"longitude\", y=\"latitude\", bands=\"wavelengths\",\n", " crs=ccrs.PlateCarree(),\n", " title=\"Stretched RGB OCI Image\",\n", " frame_width=240, frame_height=480\n", ")\n", "\n", "# Pointer stream\n", "posxy = hv.streams.PointerXY(source=rgb_map_emit, x=x_start, y=y_start)\n", "\n", "# Hover spectra function\n", "def hover_spectra(x, y):\n", " plots = {}\n", " for label, ds in datasets.items():\n", " plots[label] = (\n", " ds[\"reflectance\"]\n", " .sel(longitude=x, latitude=y, method=\"nearest\")\n", " .hvplot.line(\n", " x=\"wavelengths\", y=\"reflectance\",\n", " label=label, color=colors[label],\n", " line_dash=line_styles[label],\n", " frame_width=380,frame_height=380\n", " )\n", " )\n", " return hv.NdOverlay(plots).opts(\n", " legend_position=\"bottom\",\n", " show_legend=True\n", " )\n", "\n", "hover_dmap = hv.DynamicMap(hover_spectra, streams=[posxy])\n", "\n", "# Set Layout\n", "hv.Layout(rgb_map_emit + hover_dmap + rgb_map_oci).cols(3).opts(\n", " hv.opts.Overlay(show_legend=True, fontscale=1.3, frame_height=480),\n", ")" ] }, { "cell_type": "markdown", "id": "985dcdb0", "metadata": {}, "source": [ "### Mid-fire Time Series\n", "\n", "So far, we have focused on the individual pre- and post-fire data from EMIT and OCI, but we also have OCI data from in between those days that we can explore. This area was subject to two forest fires, one of which (the Yarram Gap Road fire) started on December 16th. Let's take a slice of the area burned in that fire and look at the median spectrum over time, using each of the 12 OCI granules we gathered through the course of the fire. \n", "\n", "Looking at the plot below, we can see that as time progressed, the fire expanded over the full slice of the park we are looking at. The spectral signatures shifted, going from a classic representation of live, healthy vegetation to those which clearly show the transition to burnt material. The next two notebooks in this series will use the features of these time series spectra to characterize the fires which took place in our ROI." ] }, { "cell_type": "code", "execution_count": null, "id": "b62cfcf3", "metadata": {}, "outputs": [], "source": [ "cmap = plt.get_cmap(\"Reds\")\n", "rgbs = cmap(np.linspace(0, 1, len(pace_ts_subsets)))\n", "\n", "fig, ax = plt.subplots(figsize=(15,7))\n", "for i in range(len(pace_ts_subsets)):\n", " ts=pace_ts_subsets[i].name[9:17]\n", " ds = xr.open_dataset(pace_ts_subsets[i])\n", "\n", " # Select section to compute the median on for the Yarram fire\n", " ds_yarram = ds[\"rhos\"].sel({\"longitude\":slice(142.2, 142.3), \n", " \"latitude\":slice(-37.28, -37.4)}).median(dim=[\"latitude\", \"longitude\"])\n", " \n", " ax.plot(ds[\"wavelength_3d\"][:-6], ds_yarram[:-6], c=rgbs[i], label=f\"{ts[4:6]}/{ts[6:]}/{ts[:4]}\")\n", " ax.plot(ds[\"wavelength_3d\"][-6:], ds_yarram[-6:], c=rgbs[i], ls='--', marker='o')\n", "ax.set_title(\"Yarram mid-fire spectra from PACE OCI\")\n", "ax.set_xlabel(\"Wavelength (nm)\")\n", "ax.set_ylabel(\"Reflectance\")\n", "ax.legend(bbox_to_anchor=(1.01, 1), loc=\"upper left\")\n", "plt.show()\n", " " ] }, { "cell_type": "markdown", "id": "561b372b", "metadata": {}, "source": [ "
\n", "\n", "You have completed the second notebook in the PACE+EMIT series. We plan to add a third notebook, `03_Wildfires_with_PACE_EMIT`, which uses the data we processed here to further investigate the wildfires which occurred in Grampians National Park, Victoria, Australia in 2024-2025.\n", "\n", "
" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.14" } }, "nbformat": 4, "nbformat_minor": 5 }