Exploring PACE OCI L2 Surface Reflectance

In this tutorial we will work with Level 2 data from the PACE’s Ocean Color Instrument (OCI), specificaclly the Surface Reflectance data product. For information on how this data product has been created please visit: PACE OCI Surface Reflectance (SFREFL) Algorithm Theortical Basis Document (ATBD)

The OCI instrument is a hyperspectral imaging radiometer, collecting hyperspectral measurements from 340 nm - 895 nm (Ultra Violet (UV) to near-infrared-NIR), and multi-spectral measurements from 940-2260 nm (NIR to shortwave-infrared).

This tutorial leverages code from Ocean Color OB DAAC Training Tutorials and NASA LP DAAC VITALS tutorials

Learning Objectives

This tutorial will guide you through working with PACE OCI L2 SFREL data on a cloud platform. Upon completion of this tutorial, you will be able to:

Access and search for available granules by time and location.
Explore data structure and format.
Visualize data by plotting individual spectral bands.
Apply cloud and water flags to effectively mask out unwanted pixels.
Create interactive RGB maps.
Display, and export spectral signatures for selected locations.
Calculate advanced hyperspectral indices using these L2 datasets.

Setup
Searching for collections and granules
Read PACE L2 Surface Reflectance data
Mask clouds and water from images
Display spectral signatures for selected locations
Calculate unique spectral indices

Setup

!conda install cf_xarray -y

Channels:
 - conda-forge
Platform: win-64
Collecting package metadata (repodata.json): - \ | / - \ | / - \ | / - \ | / - \ | / - \ | done
Solving environment: - \ | / - \ | / - \ | / done

# All requested packages already installed.

# Import required libraries
import csv
import math
from io import BytesIO

import cf_xarray
import earthaccess
import rasterio as rio
import xarray as xr
import cartopy

import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import holoviews as hv
import geoviews as gv
import hvplot.xarray

Authentication

earthaccess creates and leverages Earthdata Login tokens to authenticate with NASA systems. Earthdata Login tokens expire after a month. To retrieve a token from Earthdata Login, you can either enter your username and password each time you use earthaccess, or use a .netrc file. A .netrc file is a configuration file that is commonly used to store login credentials for remote systems. If you don’t have a .netrc or don’t know if you have one or not, you can use the persist argument with the login function below to create or update an existing one, then use it for authentication.

If you do not have an Earthdata Account, you can create one here.

auth = earthaccess.login(persist=True)

Searching for Collections

The PACE mission produces several collections or datasets available via the NASA Earthdata cloud archive.

To view what’s available, we can use the search_datasets function and with the keyword argument.

# Retrieve Collections
collections = earthaccess.search_datasets(keyword='PACE Surface Reflectance')
# Print Quantity of Results
print(f'Collections found: {len(collections)}')

Collections found: 74

If you print the collections object you can explore all of the json metadata.

# # Print collections
# collections

We can also create a list of the short-name, concept-id, and version of each result collection using list comprehension. These fields are important for specifying and searching for data within collections.

collections_info = [
    {
        'short_name': c['umm']['ShortName'],
        'collection_concept_id': c['meta']['concept-id'],
        'version': c['umm']['Version'],
        'entry_title': c['umm']['EntryTitle']
    }
    for c in collections
]
pd.set_option('display.max_colwidth', 150)
collections_info = pd.DataFrame(collections_info)
collections_info

	short_name	collection_concept_id	version	entry_title
0	PACE_OCI_L2_SFREFL	C3620139902-OB_CLOUD	3.1	PACE OCI Level-2 Regional Surface Reflectance Data, version 3.1
1	PACE_OCI_L3M_SFREFL	C3620140468-OB_CLOUD	3.1	PACE OCI Level-3 Global Mapped Surface Reflectance Data, version 3.1
2	PACE_OCI_L2_SFREFL_NRT	C3620139865-OB_CLOUD	3.1	PACE OCI Level-2 Regional Surface Reflectance - Near Real-time (NRT) Data, version 3.1
3	PACE_OCI_L3B_SFREFL	C3620140190-OB_CLOUD	3.1	PACE OCI Level-3 Global Binned Surface Reflectance Data, version 3.1
4	PACE_OCI_L3B_SFREFL_NRT	C3620140174-OB_CLOUD	3.1	PACE OCI Level-3 Global Binned Surface Reflectance - Near Real-time (NRT) Data, version 3.1
...	...	...	...	...
69	PACE_SPEXONE_L2_MAPOL_OCEAN_NRT	C4126583180-OB_CLOUD	4.0	PACE SPEXone Level-2 Regional Aerosol Over Ocean Optical Properties, FastMAPOL Algorithm - Near Real-time (NRT) Data, version 4.0
70	PVST_AIRSHARP	C3394587427-OB_DAAC	0	Airborne asSessment of Hyperspectral Aerosol optical depth and water-leaving Reflectance Product Performance for PACE
71	PACE-PAX_Ivanpah-Playa_Data	C3525273524-LARC_CLOUD	1	PACE-PAX Ivanpah Playa Data
72	AirSHARP_MetNav_AircraftInSitu_CIRPAS-TO_Data	C4126792971-LARC_CLOUD	1	AirSHARP CIRPAS Twin Otter In-Situ Meteorology and Navigation Data
73	ARCTAS_AircraftRemoteSensing_P3B_CAR_Data	C3228971208-LARC_CLOUD	1	ARCTAS P-3B Aircraft CAR Data

74 rows × 4 columns

The collection concept-id is the best way to search for data within a collection, as this is unique to each collection. The short-name can be used as well, however the version should be passed as well as there can be multiple versions available with the same short name. After finding the collection you want to search, you can use the concept-id to search for granules within that collection.

Searching for Granules

A granule can be thought of as a unique spatiotemporal grouping within a collection. To search for granules, we can use the search_data function from earthaccess and provide the arguments for our search. Its possible to specify search products using several criteria shown in the table below:

dataset origin and location	spatio temporal parameters	dataset metadata parameters
archive_center	bounding_box	concept_id
data_center	temporal	entry_title
daac	point	granule_name
provider	polygon	version
cloud_hosted	line	short_name

For this example we will search for the PACE OCI Level-2 Regional Surface Reflectance Data using a bounding box and temporal parameters, and add a cloud_cover parameter. Note that not all datasets have cloud cover information, so this parameter may not work for all datasets.

tspan = ("2025-03-01", "2025-03-20")
bbox = (-86.512, 31.70, -86.512, 34.938)
 
results = earthaccess.search_data(
    short_name="PACE_OCI_L2_SFREFL",
    version="3.1",
    temporal=tspan,
    bounding_box=bbox,
    cloud_cover=(0, 50)
)
print(f'Granules found: {len(results)}')

Granules found: 17

C:\Users\ebolch\AppData\Local\miniforge3\envs\lpdaac_vitals\Lib\site-packages\earthaccess\results.py:348: FutureWarning: As of version 1.0, `DataGranule.size` will be accessed as an attribute; e.g. use `DataCollection.size` **not** `DataCollection.size()`
  self["size"] = self.size()

Results is a list, so we can use an index to view a single result.

results[0]

C:\Users\ebolch\AppData\Local\miniforge3\envs\lpdaac_vitals\Lib\site-packages\earthaccess\results.py:375: FutureWarning: As of version 1.0, `DataGranule.size` will be accessed as an attribute; e.g. use `DataCollection.size` **not** `DataCollection.size()`
  Size(MB): {self.size()}
C:\Users\ebolch\AppData\Local\miniforge3\envs\lpdaac_vitals\Lib\site-packages\earthaccess\formatters.py:40: FutureWarning: As of version 1.0, `DataGranule.size` will be accessed as an attribute; e.g. use `DataCollection.size` **not** `DataCollection.size()`
  granule_size = round(granule.size(), 2)

Data: PACE_OCI.20250301T183734.L2.SFREFL.V3_1.nc

Size: 745.93 MB

Cloud Hosted: True

We can also retrieve specific metadata for a result using .keys() since this object also acts as a dictionary.

results[0].keys()
# results[0]['meta']

dict_keys(['meta', 'umm', 'size'])

Feel free to explore the metadata of the results, to check if there are fields relevant for your use case.

From here, we can also retrieve the links for data using earthaccess, which we can either download or stream.

links = [result.data_links()[0] for result in results]
# Display the first 5 links
links[:5]

['https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250301T183734.L2.SFREFL.V3_1.nc',
 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250302T191306.L2.SFREFL.V3_1.nc',
 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250303T181013.L2.SFREFL.V3_1.nc',
 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250304T184545.L2.SFREFL.V3_1.nc',
 'https://obdaac-tea.earthdatacloud.nasa.gov/ob-cumulus-prod-public/PACE_OCI.20250306T181826.L2.SFREFL.V3_1.nc']

To access data from NASA, you’ll need to provide your Earthdata Login credentials. When streaming this can best be done using the token or cookies set up by the earthaccess library. Since we’ve already logged in, we can start an fsspec session to manage our connection to a remote file, including sending credentials. This allows other libraries to work with a URL as if it is a local file.

In addition to fsspec for this example, we’ll also use BytesIO to read the full dataset into memory. This is not always the best approach, but for working with a single PACE scene which is ~750 mb, it speeds up the process due to the data structure and some caveats regarding fsspec.

fs = earthaccess.get_fsspec_https_session()

Now, open the connection to the remote file and read into memory.

with fs.open(links[0]) as f:
    file_content = BytesIO(f.read())

Reading dataset

The PACE_OCI_L2_SFREFL files are heirarchichal netCDF files, so we will use the xarray.datatree function to read in and view all of the groups. Using the datatree representation object, we can explore the data structure, which is useful for understanding the groups, variables, dimension order, and coordinates.

datatree = xr.open_datatree(file_content, decode_timedelta=False)
datatree

It will be easier to work with the data if we create a single dataset, which we can do by merging the data and coordinate variables we need from different groups. Longitude and latitude appear as data variables, in the group navigation_data, they need to be explicitly set as coordinates. Note the group geophysical_data and its data variables: rhos and l2_flags. The rhos variable are surface reflectances, and the l2_flags are quality flags as defined by the Ocean Biology Processing Group (OBPG).

ds = xr.merge(
    (
        datatree.ds,
        datatree["geophysical_data"].ds[["rhos", "l2_flags"]],
        datatree["sensor_band_parameters"].coords,
        datatree["navigation_data"].ds.set_coords(("longitude", "latitude")).coords,
    )
)
ds

We can also see which wavelengths we have surface reflectance measurements at by accessing the wavelength_3d coordinate:

# Check the wavelengths available for PACE 
ds["wavelength_3d"]

<xarray.DataArray 'wavelength_3d' (wavelength_3d: 122)> Size: 976B
array([ 346.,  351.,  356.,  361.,  366.,  371.,  375.,  380.,  385.,  390.,
        395.,  400.,  405.,  410.,  415.,  420.,  425.,  430.,  435.,  440.,
        445.,  450.,  455.,  460.,  465.,  470.,  475.,  480.,  485.,  490.,
        495.,  500.,  505.,  510.,  515.,  520.,  525.,  530.,  535.,  540.,
        545.,  550.,  555.,  560.,  565.,  570.,  575.,  580.,  586.,  615.,
        620.,  625.,  630.,  635.,  640.,  642.,  645.,  647.,  650.,  652.,
        655.,  657.,  660.,  662.,  665.,  667.,  670.,  672.,  675.,  677.,
        679.,  682.,  697.,  699.,  702.,  704.,  707.,  709.,  712.,  714.,
        719.,  724.,  729.,  734.,  739.,  742.,  744.,  747.,  749.,  752.,
        754.,  772.,  774.,  779.,  784.,  789.,  794.,  799.,  804.,  809.,
        814.,  819.,  824.,  829.,  835.,  840.,  845.,  850.,  855.,  860.,
        865.,  870.,  875.,  880.,  885.,  890.,  895., 1038., 1249., 1618.,
       2131., 2258.])
Coordinates:
  * wavelength_3d  (wavelength_3d) float64 976B 346.0 351.0 ... 2.258e+03
Attributes:
    long_name:  wavelengths
    units:      nm
    valid_min:  0
    valid_max:  20000

xarray.DataArray

'wavelength_3d'

wavelength_3d: 122

346.0 351.0 356.0 361.0 ... 1.249e+03 1.618e+03 2.131e+03 2.258e+03

array([ 346.,  351.,  356.,  361.,  366.,  371.,  375.,  380.,  385.,  390.,
        395.,  400.,  405.,  410.,  415.,  420.,  425.,  430.,  435.,  440.,
        445.,  450.,  455.,  460.,  465.,  470.,  475.,  480.,  485.,  490.,
        495.,  500.,  505.,  510.,  515.,  520.,  525.,  530.,  535.,  540.,
        545.,  550.,  555.,  560.,  565.,  570.,  575.,  580.,  586.,  615.,
        620.,  625.,  630.,  635.,  640.,  642.,  645.,  647.,  650.,  652.,
        655.,  657.,  660.,  662.,  665.,  667.,  670.,  672.,  675.,  677.,
        679.,  682.,  697.,  699.,  702.,  704.,  707.,  709.,  712.,  714.,
        719.,  724.,  729.,  734.,  739.,  742.,  744.,  747.,  749.,  752.,
        754.,  772.,  774.,  779.,  784.,  789.,  794.,  799.,  804.,  809.,
        814.,  819.,  824.,  829.,  835.,  840.,  845.,  850.,  855.,  860.,
        865.,  870.,  875.,  880.,  885.,  890.,  895., 1038., 1249., 1618.,
       2131., 2258.])

Coordinates: (1)

wavelength_3d

(wavelength_3d)

float64

346.0 351.0 ... 2.131e+03 2.258e+03

long_name :: wavelengths
units :: nm
valid_min :: 0
valid_max :: 20000

array([ 346.,  351.,  356.,  361.,  366.,  371.,  375.,  380.,  385.,  390.,
        395.,  400.,  405.,  410.,  415.,  420.,  425.,  430.,  435.,  440.,
        445.,  450.,  455.,  460.,  465.,  470.,  475.,  480.,  485.,  490.,
        495.,  500.,  505.,  510.,  515.,  520.,  525.,  530.,  535.,  540.,
        545.,  550.,  555.,  560.,  565.,  570.,  575.,  580.,  586.,  615.,
        620.,  625.,  630.,  635.,  640.,  642.,  645.,  647.,  650.,  652.,
        655.,  657.,  660.,  662.,  665.,  667.,  670.,  672.,  675.,  677.,
        679.,  682.,  697.,  699.,  702.,  704.,  707.,  709.,  712.,  714.,
        719.,  724.,  729.,  734.,  739.,  742.,  744.,  747.,  749.,  752.,
        754.,  772.,  774.,  779.,  784.,  789.,  794.,  799.,  804.,  809.,
        814.,  819.,  824.,  829.,  835.,  840.,  845.,  850.,  855.,  860.,
        865.,  870.,  875.,  880.,  885.,  890.,  895., 1038., 1249., 1618.,
       2131., 2258.])

Attributes: (4)
long_name :
wavelengths
units :
nm
valid_min :
0
valid_max :
20000

Note that “wavelength_3d” is an indexed coordinate, which allows us to subset the dataset by slicing or choosing individual wavelengths. The method=“nearest” argument lets us select one wavelength without knowning its exact value.

Hence, we can select one wavelength even if it does not represent an exact value in the “wavelength_3d” array and plot it using the method=“nearest” argument.

rhos_860 = ds["rhos"].sel({"wavelength_3d": 860}, method="nearest")

fig, ax = plt.subplots(figsize=(9, 5), subplot_kw={"projection": cartopy.crs.PlateCarree()})
ax.gridlines(draw_labels={"left": "y", "bottom": "x"}, linewidth=0.25)
ax.coastlines(linewidth=0.5)
ax.add_feature(cartopy.feature.OCEAN, edgecolor="w", linewidth=0.01)
ax.add_feature(cartopy.feature.LAND, edgecolor="w", linewidth=0.01)
rhos_860.plot(x="longitude", y="latitude", cmap="Greys_r", vmin=0, vmax=1.0)
plt.show()

rhos_860

<xarray.DataArray 'rhos' (number_of_lines: 1709, pixels_per_line: 1272)> Size: 9MB
[2173848 values with dtype=float32]
Coordinates:
    longitude      (number_of_lines, pixels_per_line) float32 9MB ...
    latitude       (number_of_lines, pixels_per_line) float32 9MB ...
    wavelength_3d  float64 8B 860.0
Dimensions without coordinates: number_of_lines, pixels_per_line
Attributes:
    long_name:  Surface reflectance
    valid_min:  -0.05
    valid_max:  1.5

Masking clouds

#Now flags 
ds["l2_flags"]

<xarray.DataArray 'l2_flags' (number_of_lines: 1709, pixels_per_line: 1272)> Size: 9MB
[2173848 values with dtype=int32]
Coordinates:
    longitude  (number_of_lines, pixels_per_line) float32 9MB ...
    latitude   (number_of_lines, pixels_per_line) float32 9MB ...
Dimensions without coordinates: number_of_lines, pixels_per_line
Attributes:
    long_name:      Level-2 Processing Flags
    valid_min:      -2147483648
    valid_max:      2147483647
    flag_masks:     [          1           2           4           8         ...
    flag_meanings:  ATMFAIL LAND PRODWARN HIGLINT HILT HISATZEN COASTZ SPARE ...

We will use the cf_xarray library to create a mask including only land and excluding clouds and ice.

# Create and apply Mask
ds["l2_flags"].cf.is_flag_variable
cldwater_mask = (
    (ds["l2_flags"].cf == "LAND")
& ~(ds["l2_flags"].cf == "CLDICE")
)
# Apply mask, creating new dataset
rhos = ds["rhos"].where(cldwater_mask)

# Plot a band of masked data
rhos_860 = rhos.sel({"wavelength_3d": 860}, method="nearest")

fig, ax = plt.subplots(figsize=(9, 5), subplot_kw={"projection": cartopy.crs.PlateCarree()})
ax.coastlines(linewidth=0.25)
ax.gridlines(draw_labels={"left": "y", "bottom": "x"}, linewidth=0.25)
ax.add_feature(cartopy.feature.OCEAN, edgecolor="w", linewidth=0.01)
ax.add_feature(cartopy.feature.LAND, edgecolor="w", linewidth=0.01)
ax.add_feature(cartopy.feature.LAKES, edgecolor="k", linewidth=0.1)
rhos_860.plot(x="longitude", y="latitude", cmap="Greys_r", vmin=0, vmax=1.0)
plt.show()

Projecting data

From viewing the data structure, we can see that the data is in swath format based on the dimensions. This means that we will need to grid the data to work with it in a typical manner.

There are a couple ways to do this transformation, but for this example, we’ll opt for a simple approach using rioxarray. Note that the process shown here will create a grid unique to this scene, so it won’t easily merge or stack with others created in the same way for something like a time-series analysis.

We encourage you to check the following resources for additional information on data transformation for PACE datasets: SeaDAS and Python Jupyter notebook

To grid the data with rioxarray, we need to change the order of the dimensions in the array. Hyperspectral data is often stored (y,x,band) since most processing is conducted along the band dimension; however, rioxarray expects the dimensions to be in the order of (band, y, x), so we must transpose the array to grid the data using rioxarray. Note that we’ll also want to ensure we’re using nearest neighbor resampling to preserve the spectral information during the gridding process. If we use another method we could be creating artificial spectra which could lead to inaccurate unmixing or classification results.

# Transpose data
sr_src = rhos.transpose("wavelength_3d", ...)
sr_src

<xarray.DataArray 'rhos' (wavelength_3d: 122, number_of_lines: 1709,
                          pixels_per_line: 1272)> Size: 1GB
array([[[       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        ...,
        [       nan,        nan,        nan, ..., 0.08980089,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ..., 0.09087893,
                nan,        nan]],

       [[       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
...
        [       nan,        nan,        nan, ..., 0.13494207,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ..., 0.13389756,
                nan,        nan]],

       [[       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        ...,
        [       nan,        nan,        nan, ..., 0.13306819,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ..., 0.13225715,
                nan,        nan]]], shape=(122, 1709, 1272), dtype=float32)
Coordinates:
  * wavelength_3d  (wavelength_3d) float64 976B 346.0 351.0 ... 2.258e+03
    longitude      (number_of_lines, pixels_per_line) float32 9MB -99.13 ... ...
    latitude       (number_of_lines, pixels_per_line) float32 9MB 13.6 ... 36.6
Dimensions without coordinates: number_of_lines, pixels_per_line
Attributes:
    long_name:  Surface reflectance
    valid_min:  -0.05
    valid_max:  1.5

xarray.DataArray

'rhos'

wavelength_3d: 122
number_of_lines: 1709
pixels_per_line: 1272

nan nan nan nan nan nan nan ... nan nan 0.1361 0.1331 0.1323 nan nan

array([[[       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        ...,
        [       nan,        nan,        nan, ..., 0.08980089,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ..., 0.09087893,
                nan,        nan]],

       [[       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
...
        [       nan,        nan,        nan, ..., 0.13494207,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ..., 0.13389756,
                nan,        nan]],

       [[       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        ...,
        [       nan,        nan,        nan, ..., 0.13306819,
                nan,        nan],
        [       nan,        nan,        nan, ...,        nan,
                nan,        nan],
        [       nan,        nan,        nan, ..., 0.13225715,
                nan,        nan]]], shape=(122, 1709, 1272), dtype=float32)

Coordinates: (3)

wavelength_3d

(wavelength_3d)

float64

346.0 351.0 ... 2.131e+03 2.258e+03

long_name :: wavelengths
units :: nm
valid_min :: 0
valid_max :: 20000

array([ 346.,  351.,  356.,  361.,  366.,  371.,  375.,  380.,  385.,  390.,
        395.,  400.,  405.,  410.,  415.,  420.,  425.,  430.,  435.,  440.,
        445.,  450.,  455.,  460.,  465.,  470.,  475.,  480.,  485.,  490.,
        495.,  500.,  505.,  510.,  515.,  520.,  525.,  530.,  535.,  540.,
        545.,  550.,  555.,  560.,  565.,  570.,  575.,  580.,  586.,  615.,
        620.,  625.,  630.,  635.,  640.,  642.,  645.,  647.,  650.,  652.,
        655.,  657.,  660.,  662.,  665.,  667.,  670.,  672.,  675.,  677.,
        679.,  682.,  697.,  699.,  702.,  704.,  707.,  709.,  712.,  714.,
        719.,  724.,  729.,  734.,  739.,  742.,  744.,  747.,  749.,  752.,
        754.,  772.,  774.,  779.,  784.,  789.,  794.,  799.,  804.,  809.,
        814.,  819.,  824.,  829.,  835.,  840.,  845.,  850.,  855.,  860.,
        865.,  870.,  875.,  880.,  885.,  890.,  895., 1038., 1249., 1618.,
       2131., 2258.])

longitude

(number_of_lines, pixels_per_line)

float32

-99.13 -99.06 ... -76.79 -76.71

long_name :: Longitude
units :: degrees_east
standard_name :: longitude
valid_min :: -180.0
valid_max :: 180.0

array([[ -99.129036,  -99.064575,  -99.00084 , ...,  -74.639114,
         -74.573906,  -74.50793 ],
       [ -99.132484,  -99.06801 ,  -99.004265, ...,  -74.641335,
         -74.576126,  -74.51015 ],
       [ -99.13476 ,  -99.07029 ,  -99.00655 , ...,  -74.6423  ,
         -74.57708 ,  -74.51109 ],
       ...,
       [-105.207855, -105.13512 , -105.06324 , ...,  -76.86469 ,
         -76.78669 ,  -76.70763 ],
       [-105.212906, -105.13999 , -105.06777 , ...,  -76.86625 ,
         -76.788216,  -76.709015],
       [-105.2168  , -105.142365, -105.07115 , ...,  -76.86663 ,
         -76.78864 ,  -76.70912 ]], shape=(1709, 1272), dtype=float32)

latitude

(number_of_lines, pixels_per_line)

float32

13.6 13.61 13.62 ... 36.59 36.6

long_name :: Latitude
units :: degrees_north
standard_name :: latitude
valid_min :: -90.0
valid_max :: 90.0

array([[13.596701, 13.608392, 13.619964, ..., 18.716452, 18.731403,
        18.746536],
       [13.606942, 13.618638, 13.630214, ..., 18.726784, 18.741734,
        18.756865],
       [13.617398, 13.629096, 13.640673, ..., 18.737404, 18.752354,
        18.76748 ],
       ...,
       [31.077051, 31.094622, 31.111973, ..., 36.56201 , 36.57278 ,
        36.583675],
       [31.087013, 31.10461 , 31.122004, ..., 36.572365, 36.583122,
        36.59406 ],
       [31.097267, 31.11505 , 31.132324, ..., 36.582874, 36.59362 ,
        36.604626]], shape=(1709, 1272), dtype=float32)

Attributes: (3)
long_name :
Surface reflectance
valid_min :
-0.05
valid_max :
1.5

# Set spatial dimensions and swath CRS
sr_src = sr_src.rio.set_spatial_dims("pixels_per_line", "number_of_lines").rio.write_crs("epsg:4326")
sr_src

# Grid the data using rioxarray and rename dimensions
sr_dst = sr_src.rio.reproject(
    dst_crs=sr_src.rio.crs,
    src_geoloc_array=(
        sr_src.coords["longitude"],
        sr_src.coords["latitude"],
    ),
    nodata=np.nan,
    resampling=rio.warp.Resampling.nearest,
).rename({'y': 'latitude', 'x': 'longitude'})

sr_dst

From this point, we can either subset the data spatially to a smaller region of interest, or extract a series of location points to work with.

Interactive Spatial and Spectral Plots

Combining the Spatial and Spectral information into a single visualization can be a powerful tool for exploring and inspecting imaging spectroscopy data. Using the streams module from Holoviews we can link a spatial map to a plot of spectra.

We could plot a single band image as we previously have, but using a multiband image, like an RGB may help infer what targets we’re examining. Build function to select and adjust the brightness of bands to create a nice RGB image.

Select bands to represent red (650 nm), green (560 nm), and blue (470 nm) by finding the nearest to a wavelength chosen to represent that color.

rgb_da = sr_dst.sel(wavelength_3d=[650, 560, 470], method='nearest')
rgb_da

To get a better rgb image, we can balance the brightness of the selected wavelengths to make a prettier map. This will not affect the data, just the visuals. To do this we will use the function below. We can change the bright argument to increase or decrease the brightness of the scene as a whole.

def gamma_adjust(rgb_da, bright=0.3):
    """
    Adjust gamma across all bands in the RGB dataset.
    """
    array = rgb_da.data
    # Mask nan and negative values
    invalid = np.isnan(array) | (array < 0)
    valid = ~invalid

    # Calculate gamma based on the mean of valid values
    mean_valid = np.nanmean(array[valid])
    gamma = math.log(bright) / math.log(mean_valid)
    # Apply scaling and clip
    scaled = np.full_like(array, np.nan)
    scaled[valid] = np.power(array[valid], gamma)
    rgb_da.data = np.clip(scaled, 0, 1)
    return rgb_da

rgb_da = gamma_adjust(rgb_da, bright=0.3)

rgb_da.hvplot.rgb(rasterize=True, dynamic=False, x='longitude', y='latitude', bands='wavelength_3d', frame_height=400, geo=True, crs='EPSG:4326', tiles='OSM', title="Stretched RGB PACE Image")

Now that we have an RGB dataset, we can use that to create a spatial plot, and data selected by clicking on our rgb map, which can be inputs for a function to return values from the full dataset at the selected latitude and longitude location using the cell below. To visualize the spectral and spatial data side-by-side, we use the Point Draw tool from the holoviews library.

Click in the RGB image to add spectra to the plot. You can also click and hold the mouse button then drag previously placed points. To remove a point click and hold the mouse button down, then press the backspace key.

# Starting Point from center pixel
x_start = -85.3879
y_start = 31.3446

# Set Point Limit
POINT_LIMIT = 10

# Set up Color Cycling
color_cycle = hv.Cycle('Category20')

first_point = ([x_start], [y_start], [0])

points = gv.Points(first_point, kdims=['longitude','latitude'], vdims='id', crs=cartopy.crs.PlateCarree())

points_stream = hv.streams.PointDraw(
    data=points.columns(),
    source=points,
    drag=True,
    num_objects=POINT_LIMIT,
    styles={'fill_color': color_cycle.values[:POINT_LIMIT], 'line_color': 'gray'}
)

# RGB Plot without Basemap

rgb_map = rgb_da.hvplot.rgb(x='longitude', y='latitude', bands='wavelength_3d',
                            frame_height=480, frame_width=480,
                            crs=cartopy.crs.PlateCarree(),
                            title="Stretched RGB PACE Image")

# Coastlines
coastlines = gv.feature.coastline(projection=cartopy.crs.PlateCarree()).opts(line_color='black', line_width=1)

posxy = hv.streams.PointerXY(source=rgb_map, x=x_start, y=y_start)
clickxy = hv.streams.Tap(source=rgb_map, x=x_start, y=y_start)

# Function to build spectral plot of clicked location to show on hover stream plot
def click_spectra(data):
    coordinates = [c for c in zip(data['longitude'], data['latitude'])]
    
    plots = {}
    for i, coords in enumerate(coordinates):
        x, y = coords
        selected = sr_dst.sel(longitude=x, latitude=y, method="nearest")
        plots[i] = (
            selected.hvplot.line(
                y="rhos",
                x="wavelength_3d",
                xlim=(340, 895), #you can comment or change the range here to see the multispectral bands too
                label=f"{i}"
            )
        )
        points_stream.data["id"][i] = i
    return hv.NdOverlay(plots).opts(hv.opts.Curve(color=color_cycle))

def hover_spectra(x,y):
    return sr_dst.sel(longitude=x,latitude=y,method='nearest').hvplot.line(y='rhos',x='wavelength_3d',
                                                                           color='black', frame_width=480)
# Define the Dynamic Maps
click_dmap = hv.DynamicMap(click_spectra, streams=[points_stream])
hover_dmap = hv.DynamicMap(hover_spectra, streams=[posxy])

# Plot the Map and Dynamic Map side by side
hv.Layout(hover_dmap*click_dmap + rgb_map * coastlines* points).cols(2).opts(
    hv.opts.Points(active_tools=['point_draw'], size=10, tools=['hover'], color='white', line_color='gray'),
    hv.opts.Overlay(show_legend=False, show_title=False, fontscale=1.5, frame_height=480)
)

We can take these selected points and the corresponding reflectance spectra and save them as a .csv for later use.

data = points_stream.data
wavelengths = sr_dst.wavelength_3d.values

rows = [["id", "longitude", "latitude"] + [str(i) for i in wavelengths]]
 
for p in zip(data['longitude'], data['latitude'], data['id']):
    x, y, i = p
    spectra = sr_dst.sel(longitude=x, latitude=y, method="nearest").values
    row = [i, x, y] + list(spectra)
    rows.append(row)

with open('../data/pace_interactive_plot_data.csv', 'w', newline='') as f:
    writer = csv.writer(f)
    writer.writerows(rows)

# In case you cannot save the file, check your directory
# import os
# print(os.getcwd())

Calculate unique indices

Hyperspectral-enabled Vegetation Indices (VI)

We’ll take the Chlorophyll Index Red Edge (CIRE) as an example of a hyperspectral-enabled VI. CIRE uses bands from the red edge and the NIR to get at relative canopy chlorophyll content.

CIRE: (ρ 800 /ρ 705)−1

Carotenoid Content Index (Car): [(1/ρ495)−(1/ρ705)] * ρ800

Information on Hyperspectral enabled indices by PACE here.

## Chlorophyll Index Red Edge (CIRE)
# Select bands
sr_800 = sr_dst.sel({"wavelength_3d": 800}, method="nearest")
sr_705 = sr_dst.sel({"wavelength_3d": 705}, method="nearest")
#Calculate
cire = (sr_800 / sr_705) - 1
cire.attrs["long_name"] = "Chlorophyll Index Red Edge (CIRE)"

cire

<xarray.DataArray 'rhos' (latitude: 1524, longitude: 2037)> Size: 12MB
array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]],
      shape=(1524, 2037), dtype=float32)
Coordinates:
  * latitude     (latitude) float64 12kB 36.62 36.6 36.59 ... 13.63 13.62 13.6
  * longitude    (longitude) float64 16kB -105.2 -105.2 -105.2 ... -74.46 -74.44
    spatial_ref  int64 8B 0
Attributes:
    long_name:   Chlorophyll Index Red Edge (CIRE)
    valid_min:   -0.05
    valid_max:   1.5
    _FillValue:  nan

xarray.DataArray

'rhos'

latitude: 1524
longitude: 2037

nan nan nan nan nan nan nan nan ... nan nan nan nan nan nan nan nan

array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]],
      shape=(1524, 2037), dtype=float32)

Coordinates: (3)
- latitude
  (latitude)
  float64
  36.62 36.6 36.59 ... 13.62 13.6
  axis :
  Y
  long_name :
  latitude
  standard_name :
  latitude
  units :
  degrees_north
```
array([36.61871 , 36.603598, 36.588485, ..., 13.632375, 13.617263, 13.60215 ],
      shape=(1524,))
```
- longitude
  (longitude)
  float64
  -105.2 -105.2 ... -74.46 -74.44
  axis :
  X
  long_name :
  longitude
  standard_name :
  longitude
  units :
  degrees_east
```
array([-105.213132, -105.198019, -105.182906, ...,  -74.474009,  -74.458897,
        -74.443784], shape=(2037,))
```
- spatial_ref
  ()
  int64
  0
  crs_wkt :
  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
  semi_major_axis :
  6378137.0
  semi_minor_axis :
  6356752.314245179
  inverse_flattening :
  298.257223563
  reference_ellipsoid_name :
  WGS 84
  longitude_of_prime_meridian :
  0.0
  prime_meridian_name :
  Greenwich
  geographic_crs_name :
  WGS 84
  horizontal_datum_name :
  World Geodetic System 1984
  grid_mapping_name :
  latitude_longitude
  spatial_ref :
  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
  GeoTransform :
  -105.22068786621094 0.015112646133606075 0.0 36.62626647949219 0.0 -0.015112646133606075
```
array(0)
```
Attributes: (4)
long_name :
Chlorophyll Index Red Edge (CIRE)
valid_min :
-0.05
valid_max :
1.5
_FillValue :
nan

#Plot
map1 = cire.hvplot.image(rasterize=True, dynamic=False, x='longitude', y='latitude', frame_height=480, geo=True, cmap="viridis", tiles='CartoDark', title="CIRE from PACE OCI")
map1

## Carotenoid Content Index (Car) 
# Select additional bands
sr_495 = sr_dst.sel({"wavelength_3d": 495}, method="nearest")

#Calculate
car = ((1 / sr_495)- (1 / sr_705)) * sr_800
car.attrs["long_name"] = "Caretenoid Content Index (Car)"

map2 = car.hvplot.image(rasterize=True, dynamic=False, x='longitude', y='latitude', frame_height=480, geo=True, cmap="plasma", tiles='CartoDark', title="CAR from PACE OCI")
map2

# Compare CIRE and CAR side to side
map1.opts(frame_height=400, frame_width=480)+map2.opts(frame_height=400, frame_width=480)

Enjoy PACE data

Calculate your own vegetation index

Select relevant bands
Calculate index
Plot data

Contact Info:

Email: LPDAAC@usgs.gov
Voice: +1-866-573-3222
Organization: Land Processes Distributed Active Archive Center (LP DAAC)¹
Website: https://www.earthdata.nasa.gov/centers/lp-daac

¹Work performed under USGS contract 140G0126D0001 for NASA contract NNG14HH33I.

Learning Objectives

Contents

Setup

Authentication

Searching for Collections

Searching for Granules

Reading dataset

Masking clouds

Projecting data

Interactive Spatial and Spectral Plots

Calculate unique indices

Hyperspectral-enabled Vegetation Indices (VI)

Enjoy PACE data

Calculate your own vegetation index

Contact Info: