Problem with reading NASA PACE data

[giswqs]

Cross posted NASA/EMIT-Data-Resources

Thank you for the excellent notebook examples from the repo. @bingqing-liu and I are building a Python package called HyperCoast for visualizing NASA Hyperspectral data interactively with a few lines of code. Here is a notebook example for visualizing EMIT data interactively. However, we run into issues with reading PACE data. The latitude and longitude coordinates are in the format of (latutude, longitude). We need to turn them into one dimension like the EMIT dataset shown in the example below. The goal is the turn the xarray dataset into a rasterio dataset so that we can plot it on the interactive map.

Any advice? @alexgleith @betolink @BriannaLind

PACE

Code snippet

import xarray as xr

# download from https://github.com/opengeos/datasets/releases/download/netcdf/PACE_OCI.20240423T184658.L2.OC_AOP.V1_0_0.NRT.nc
filepath = 'data/PACE_OCI.20240423T184658.L2.OC_AOP.V1_0_0.NRT.nc'
ds = xr.open_dataset(filepath, group="geophysical_data")
ds = ds.swap_dims({'number_of_lines': 'latitude', 'pixels_per_line': 'longitude'})
wvl = xr.open_dataset(filepath, group="sensor_band_parameters")
loc = xr.open_dataset(filepath, group="navigation_data")
lat = loc.latitude
lon = loc.longitude 
wavelengths = wvl.wavelength_3d
Rrs = ds.Rrs 

dataset = xr.Dataset(
    {
        "Rrs": (("latitude", "longitude", "wavelengths"), Rrs.data)
    },
    coords={
        "latitude": (("latitude", "longitude"), lat.data),
        "longitude": (("latitude", "longitude"), lon.data),
        "wavelengths": ("wavelengths", wavelengths.data)
    }
)
dataset

Note the two dimensions (latitude, longitude) in the latitude and longitude coordinates

EMIT

Code snippet

import hypercoast
# download from https://github.com/opengeos/datasets/releases/download/netcdf/EMIT_L2A_RFL_001_20240404T161230_2409511_009.nc
filepath = "data/EMIT_L2A_RFL_001_20240404T161230_2409511_009.nc"
dataset = hypercoast.read_emit(filepath)
dataset

EMIT demo

[itcarroll]

If I look at the original Level-2 files from EMIT, it has a similar structure as PACE/OCI level-2 products.

results = earthaccess.search_data(short_name="EMITL2ARFL", count=1)
paths = earthaccess.open(results)
xr.open_dataset(paths[0], group="location")

<xarray.Dataset> Size: 101MB
Dimensions:  (downtrack: 1280, crosstrack: 1242, ortho_y: 1890, ortho_x: 2062)
Dimensions without coordinates: downtrack, crosstrack, ortho_y, ortho_x
Data variables:
    lon      (downtrack, crosstrack) float64 13MB ...
    lat      (downtrack, crosstrack) float64 13MB ...
    elev     (downtrack, crosstrack) float64 13MB ...
    glt_x    (ortho_y, ortho_x) float64 31MB ...
    glt_y    (ortho_y, ortho_x) float64 31MB ...

Do you know where to find how the EMIT product you have from opengeos/datasets was processed to a rectilinear grid?

[giswqs]

I adopted the source code from this notebook to generate the gridded data

https://github.com/nasa/EMIT-Data-Resources/blob/main/python/tutorials/Exploring_EMIT_L2A_Reflectance.ipynb

[giswqs]

Relevant function: emit_xarray

https://hypercoast.org/emit/#hypercoast.emit.emit_xarray

[alexgleith]

These are a bit out of date, but they are my examples for loading EMIT: https://github.com/auspatious/hyperspectral-notebooks

To be opinionated, we should be consistent in calling the variables latitude, longitude and wavelength without the pluralising of wavelength.

Then you can select with something like data.sel(wavelength=420, method="nearest") and it reads like it makes sense! And it's consistent with longitude or x ...

Sorry, but I'm not much help here with the PACE data... I haven't worked with it.

This looks like the cause of the issue, though?

dataset = xr.Dataset(
    {
        "Rrs": (("latitude", "longitude", "wavelengths"), Rrs.data)
    },
    coords={
        "latitude": (("latitude", "longitude"), lat.data),
        "longitude": (("latitude", "longitude"), lon.data),
        "wavelengths": ("wavelengths", wavelengths.data)
    }
)
dataset

Is there any reason you don't do:

dataset = xr.Dataset(
    {
        "Rrs": (("latitude", "longitude", "wavelengths"), Rrs.data)
    },
    coords={
        "latitude": ("latitude", lat.data),
        "longitude": ("latitude", lon.data),
        "wavelengths": ("wavelengths", wavelengths.data)
    }
)
dataset

[itcarroll]

That lat.data array has two dimensions, so "latitude": ("latitude", lat.data) doesn't make sense and would raise an exception. In these files, latitude and longitude are what you would call "auxiliary coordinate variables" using the CF Conventions. They are not coordinates that XArray can map to a pandas index.

Here is a figure showing the location of each pixel in a PACE/OCI L2 swath, that I hope will help this conversation. Notice that the longitude of each pixel is not the same for all pixels in a row. Same for the latitude.

NASA's level 1 and level 2 data are typically geolocated like this and not projected onto a regular grid (which is what makes both the "x" and "y" coordinates one-dimensional). We will have a tutorial up on that same website that projects L2 data soon. For the time being, you may want to look at the global L3 mapped files.

[BzGEO]

Great work as always, @giswqs! FWIW, and for comparison with the EMIT, Hyperion, and select AVIRIS scenes already available in GEE, and in the framework of a PACE Early Adopters project, I have loaded a large chunk of the PACE OCI land-related datasets (surface reflectance + vegetation indices) into GEE and am making them available via the following repo (https://bit.ly/gee_repo_pace_oci) and its dependencies. (I also made a GitHub repo for said GEE code: https://github.com/BzGEO/pace_oci_toolkit/.)

[giswqs]

@BzGEO That's very cool. Thanks for sharing.