feat: add JPSS CrIS FSR Level 1 spectral radiance DataFrameSource (PAINFUL)

[NickGeneva]

Description

Add JPSS_CRIS DataFrameSource for JPSS Cross-track Infrared Sounder (CrIS) Full Spectral Resolution (FSR) Level 1 spectral radiance observations from NOAA Open Data on AWS.

This complements the existing JPSS_ATMS microwave sounder source by adding the infrared hyperspectral sounder from the same JPSS satellite constellation (NOAA-20, NOAA-21, Suomi NPP).

Data source details

Property	Value
Source type	DataFrameSource
Remote store	https://registry.opendata.aws/noaa-jpss/
Format	HDF5 (paired SDR + GEO files per granule)
Spectral channels	2223 channels across 3 bands (LWIR 717, MWIR 869, SWIR 637)
Temporal resolution	~32s per granule, continuous polar-orbiting swaths
Date range	2023-09-06 to present
Region	Global
Authentication	Anonymous (public S3 buckets)
Satellites	NOAA-20 (n20), NOAA-21 (n21), Suomi NPP (npp)

Data licensing

License: NOAA Open Data Dissemination (NODD) — public domain
URL: https://registry.opendata.aws/noaa-jpss/

NOAA satellite data on AWS is freely available for commercial and non-commercial use with no restrictions.

Dependencies added

No new dependencies needed — uses h5py, s3fs, pandas, pyarrow, numpy which are already in core deps.

Changes

• New data source file: earth2studio/data/jpss_cris.py (~1315 lines)
• New lexicon: JPSSCrISLexicon in earth2studio/lexicon/jpss.py (2223 channel mappings)
• New vocab entry: crisfsr added to E2STUDIO_VOCAB
• Unit tests: test/data/test_jpss_cris.py with mock and network tests
• Documentation: added to datasources_dataframe.rst
• CHANGELOG.md updated

Key implementation details

• Async S3 listing + parallel download of HDF5 granule pairs (SDR + GEO)
• Thread-pool based HDF5 decoding (h5py is not async-compatible)
• Quality filtering via CrIS QF3 bit flags
• Granule-level temporal subsampling via subsample parameter
• Time tolerance support for asymmetric observation windows
• Hamming apodization (3-tap [0.23, 0.54, 0.23] kernel) with optional bypass via apodize=False
• Brightness temperature conversion via inverse Planck function using CRTM SpcCoeff constants

CrIS channel mapping — validation against UFS/GSI

The channel_index column uses the GSI sensor_chan numbering convention (1–2211 contiguous for science channels) to enable direct comparison with UFSObsSat CrIS observations.

Band	Array indices	Content	sensor_chan
LWIR	0–1	Guard (648.75–649.375 cm⁻¹)	0 (sentinel)
LWIR	2–714	Science (650.0–1095.0 cm⁻¹)	1–713
LWIR	715–716	Guard (1095.625–1096.25 cm⁻¹)	0 (sentinel)
MWIR	0–1	Guard (1208.75–1209.375 cm⁻¹)	0 (sentinel)
MWIR	2–866	Science (1210.0–1750.0 cm⁻¹)	714–1578
MWIR	867–868	Guard (1750.625–1751.25 cm⁻¹)	0 (sentinel)
SWIR	0–1	Guard (2153.75–2154.375 cm⁻¹)	0 (sentinel)
SWIR	2–634	Science (2155.0–2550.0 cm⁻¹)	1579–2211
SWIR	635–636	Guard (2550.625–2551.25 cm⁻¹)	0 (sentinel)

(Guard bands dropped after apodization)

Apodization verification

The Hamming apodization implementation was verified against the CrIS SDR ATBD (JPSS 474-00032, Section 3.7.2, Equation 54). The 3-tap spectral convolution [0.23, 0.54, 0.23] is mathematically equivalent to interferogram-domain Hamming multiplication — confirmed numerically to machine precision (~1e-14) via DCT-I round-trip testing. Boundary handling at band edges uses reflect padding, which differs from the ATBD's one-sided formula, but this only affects guard channels that are trimmed in the output.

Verification sources consulted

• GSI Fortran: read_cris.f90, setuprad.f90, crtm_interface.f90 — confirmed BUFR CRCHNM channel numbers match CRTM sensor_channel via direct integer equality (no offset/transform)
• CrIS SDR Userguide: — SDR output contains 717/869/637 channels starting at band edge (guard channels included in output) Table 4 and Section 4.1
• CrIS SDR Algorithm Theoretical Basis Document Equation 52-56, Figure 3.17 "He also showed that the optimum value of a was a function of the number of points in the spectrum; however, the optimum value of a converged to a = 0.23 for more than 100 points. "

Checklist

• I am familiar with the Contributing Guidelines.
• New or existing tests cover these changes.
• The documentation is up to date with these changes.
• The CHANGELOG.md is up to date with these changes.
• An issue is linked to this pull request.
• Assess and address Greptile feedback (AI code review bot).

Dependencies

None — all required packages are already in core dependencies.

[NickGeneva]

Sanity-Check Validation

Source: JPSS_CRIS — CrIS FSR Level 1 spectral radiance (2223 channels) compared against UFS
Satellites: NOAA-20 (n20), NOAA-21 (n21), Suomi NPP (npp)
Format: HDF5 paired SDR + GEO granules from NOAA S3

There a total 2k channels in the CRIS raw data, but only 100ish in UFS. So these are comparing those channel with apodization

Select scatter plots

Plot Script

import os
from math import ceil

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
from loguru import logger
from tqdm import tqdm

logger.remove()
logger.add(lambda msg: tqdm.write(msg, end=""), colorize=True)

from earth2studio.data import JPSS_CRIS, UFSObsSat, fetch_dataframe

os.makedirs("outputs", exist_ok=True)

# ----- Configuration --------------------------------------------------------
analysis_time = np.array([np.datetime64("2024-01-01T00:00")])
CHANNELS_PER_IMAGE = 10

jpss_source = JPSS_CRIS(
    satellites=["n20"],
    time_tolerance=np.timedelta64(1, "h"),
    apodize=True,
)
ufs_source = UFSObsSat(
    satellites=["n20"],
    time_tolerance=np.timedelta64(1, "h"),
)

# ----- Fetch ----------------------------------------------------------------
jpss_df = fetch_dataframe(
    jpss_source, time=analysis_time, variable=np.array(["crisfsr"])
)
logger.info(f"JPSS CrIS: {len(jpss_df):,} rows")

ufs_df = fetch_dataframe(ufs_source, time=analysis_time, variable=np.array(["crisfsr"]))
logger.info(f"UFS CrIS:  {len(ufs_df):,} rows")

# ----- Filter NaN / invalid values ------------------------------------------
n_before = len(jpss_df)
jpss_df = jpss_df[jpss_df["observation"].notna()].reset_index(drop=True)
logger.info(
    f"JPSS after NaN filter: {len(jpss_df):,} rows "
    f"(removed {n_before - len(jpss_df):,} NaN values)"
)

# ----- Find common channels -------------------------------------------------
jpss_channels = set(jpss_df["channel_index"].unique())
ufs_channels = sorted(ufs_df["channel_index"].unique())
common_channels = sorted(jpss_channels & set(ufs_channels))
logger.info(
    f"JPSS channels: {len(jpss_channels)}, UFS channels: {len(ufs_channels)}, "
    f"common: {len(common_channels)}"
)

if len(common_channels) == 0:
    logger.error("No common channels between JPSS and UFS — cannot plot")
    raise SystemExit(1)

# Split into groups of CHANNELS_PER_IMAGE
n_images = ceil(len(common_channels) / CHANNELS_PER_IMAGE)
channel_groups = [
    common_channels[i * CHANNELS_PER_IMAGE : (i + 1) * CHANNELS_PER_IMAGE]
    for i in range(n_images)
]

logger.info(
    f"Plotting {len(common_channels)} channels across {n_images} images "
    f"({CHANNELS_PER_IMAGE} channels each)"
)

# ----- Per-channel obs counts (diagnostic) ----------------------------------
for sc in common_channels:
    n_jpss = int((jpss_df["channel_index"] == sc).sum())
    n_ufs = int((ufs_df["channel_index"] == sc).sum())
    logger.info(f"  sc {sc}: JPSS={n_jpss:,}  UFS={n_ufs:,}")


# ----- Wavenumber helper (for labels only) ----------------------------------
def _sensor_chan_to_wn(sc: int) -> float:
    """Convert GSI sensor_chan to wavenumber (cm-1)."""
    if sc <= 713:
        return 650.0 + (sc - 1) * 0.625
    elif sc <= 1578:
        return 1210.0 + (sc - 714) * 0.625
    else:
        return 2155.0 + (sc - 1579) * 0.625


def _band_label(sc: int) -> str:
    wn = _sensor_chan_to_wn(sc)
    if sc <= 713:
        return f"LWIR sc{sc}\n{wn:.1f} cm$^{{-1}}$"
    elif sc <= 1578:
        return f"MWIR sc{sc}\n{wn:.1f} cm$^{{-1}}$"
    else:
        return f"SWIR sc{sc}\n{wn:.1f} cm$^{{-1}}$"


# ----- Plot each image ------------------------------------------------------
rng = np.random.default_rng(42)
projection = ccrs.Robinson()

for img_idx, channels in enumerate(channel_groups):
    n_rows = len(channels)
    plt.close("all")
    fig, axes = plt.subplots(
        n_rows,
        2,
        subplot_kw={"projection": projection},
        figsize=(16, 3.0 * n_rows),
    )
    if n_rows == 1:
        axes = axes[np.newaxis, :]

    for row, sc in enumerate(channels):
        jpss_ch = jpss_df[jpss_df["channel_index"] == sc]
        ufs_ch = ufs_df[ufs_df["channel_index"] == sc]

        # Shared colour limits across both sources (1st-99th percentile)
        all_obs = np.concatenate(
            [
                jpss_ch["observation"].dropna().values,
                ufs_ch["observation"].dropna().values,
            ]
        )
        if len(all_obs) > 0:
            vmin = float(np.nanpercentile(all_obs, 1))
            vmax = float(np.nanpercentile(all_obs, 99))
        else:
            vmin, vmax = 180.0, 320.0

        for col, (label, df_ch, pt_size) in enumerate(
            [
                ("JPSS (Tb)", jpss_ch, 0.3),
                ("UFS/GSI (Tb)", ufs_ch, 4.0),
            ]
        ):
            ax = axes[row, col]
            ax.set_global()
            ax.coastlines(linewidth=0.4)
            ax.gridlines(linewidth=0.2, alpha=0.4)

            if len(df_ch) > 0:
                lons = df_ch["lon"].values
                lats = df_ch["lat"].values
                obs = df_ch["observation"].values

                # Subsample for plotting performance
                max_points = 50000
                if len(obs) > max_points:
                    idx = rng.choice(len(obs), size=max_points, replace=False)
                    lons, lats, obs = lons[idx], lats[idx], obs[idx]

                im = ax.scatter(
                    lons,
                    lats,
                    c=obs,
                    s=pt_size,
                    alpha=0.6,
                    cmap="turbo",
                    vmin=vmin,
                    vmax=vmax,
                    transform=ccrs.PlateCarree(),
                    rasterized=True,
                )
                cbar = fig.colorbar(
                    im,
                    ax=ax,
                    orientation="vertical",
                    shrink=0.7,
                    pad=0.02,
                    aspect=12,
                )
                cbar.set_label("Brightness Temperature (K)", fontsize=6)
                cbar.ax.tick_params(labelsize=5)
            else:
                ax.scatter([], [], c=[], cmap="turbo", transform=ccrs.PlateCarree())

            if col == 0:
                ax.text(
                    -0.08,
                    0.5,
                    _band_label(sc),
                    fontsize=7,
                    va="center",
                    ha="center",
                    rotation="vertical",
                    rotation_mode="anchor",
                    transform=ax.transAxes,
                )
            if row == 0:
                ax.set_title(label, fontsize=12, fontweight="bold")

            ax.text(
                0.98,
                0.02,
                f"n={len(df_ch):,}",
                fontsize=6,
                ha="right",
                va="bottom",
                transform=ax.transAxes,
                bbox=dict(boxstyle="round,pad=0.2", fc="white", alpha=0.7, lw=0.3),
            )

    sc_lo = channels[0]
    sc_hi = channels[-1]
    fig.suptitle(
        f"JPSS vs UFS CrIS FSR Brightness Temperature (n20) — "
        f"{str(analysis_time[0])[:16]} UTC (±1 h)\n"
        f"Image {img_idx + 1}/{n_images}: sensor_chan {sc_lo}–{sc_hi}",
        fontsize=14,
        y=1.0,
    )
    plt.tight_layout()
    fname = f"outputs/jpss_vs_ufs_cris_{img_idx + 1:02d}.jpg"
    plt.savefig(fname, dpi=150, bbox_inches="tight")
    logger.info(f"Saved {fname}")

logger.info(f"Done — {n_images} images saved to outputs/")

Histograms

Data is from 2024-01-01T00:00 within a +/- 1hr range. I believe the delta is due to the GIS pipeline filtering / sub sampling removing typical "cooler" values on the outer scan angles.

Without the hamming functions, the histograms are way off, so it has a pretty large effect.

Plot Script

```python

import os
from math import ceil

import matplotlib.pyplot as plt
import numpy as np
from loguru import logger
from tqdm import tqdm

logger.remove()
logger.add(lambda msg: tqdm.write(msg, end=""), colorize=True)

from earth2studio.data import JPSS_CRIS, UFSObsSat, fetch_dataframe # noqa: E402

os.makedirs("outputs", exist_ok=True)

----- Configuration --------------------------------------------------------

analysis_time = np.array([np.datetime64("2024-01-01T00:00")])
GRID_ROWS = 8
GRID_COLS = 8
CHANNELS_PER_PAGE = GRID_ROWS * GRID_COLS # 64
HIST_BINS = 80
HIST_ALPHA = 0.55

----- Data sources ---------------------------------------------------------

jpss_source = JPSS_CRIS(
satellites=["n20"],
time_tolerance=np.timedelta64(1, "h"),
apodize=True,
)
ufs_source = UFSObsSat(
satellites=["n20"],
time_tolerance=np.timedelta64(1, "h"),
)

----- Fetch ----------------------------------------------------------------

logger.info("Fetching JPSS CrIS data ...")
jpss_df = fetch_dataframe(
jpss_source, time=analysis_time, variable=np.array(["crisfsr"])
)
logger.info(f"JPSS CrIS: {len(jpss_df):,} rows")

logger.info("Fetching UFS CrIS data ...")
ufs_df = fetch_dataframe(ufs_source, time=analysis_time, variable=np.array(["crisfsr"]))
logger.info(f"UFS CrIS: {len(ufs_df):,} rows")

----- Filter NaN / invalid --------------------------------------------------

jpss_df = jpss_df[jpss_df["observation"].notna()].reset_index(drop=True)
ufs_df = ufs_df[ufs_df["observation"].notna()].reset_index(drop=True)

----- Find common channels (plot every UFS channel) -------------------------

jpss_channels = set(jpss_df["channel_index"].unique())
ufs_channels = sorted(ufs_df["channel_index"].unique())
common_channels = sorted(jpss_channels & set(ufs_channels))

logger.info(
f"JPSS channels: {len(jpss_channels)}, UFS channels: {len(ufs_channels)}, "
f"common: {len(common_channels)}"
)

if len(common_channels) == 0:
logger.error("No common channels between JPSS and UFS — cannot plot")
raise SystemExit(1)

----- Helpers ---------------------------------------------------------------

def _sensor_chan_to_wn(sc: int) -> float:
"""Convert GSI sensor_chan to wavenumber (cm-1)."""
if sc <= 713:
return 650.0 + (sc - 1) * 0.625
elif sc <= 1578:
return 1210.0 + (sc - 714) * 0.625
else:
return 2155.0 + (sc - 1579) * 0.625

def _band_name(sc: int) -> str:
if sc <= 713:
return "LWIR"
elif sc <= 1578:
return "MWIR"
return "SWIR"

----- Split into pages of 64 ------------------------------------------------

n_pages = ceil(len(common_channels) / CHANNELS_PER_PAGE)
channel_pages = [
common_channels[i * CHANNELS_PER_PAGE : (i + 1) * CHANNELS_PER_PAGE]
for i in range(n_pages)
]

logger.info(
f"Plotting {len(common_channels)} channels across {n_pages} page(s) "
f"({GRID_ROWS}x{GRID_COLS} = {CHANNELS_PER_PAGE} channels per page)"
)

----- Plot -------------------------------------------------------------------

for page_idx, channels in enumerate(channel_pages):
plt.close("all")
fig, axes = plt.subplots(
GRID_ROWS,
GRID_COLS,
figsize=(28, 24),
constrained_layout=True,
)
axes_flat = axes.flatten()

for i, sc in enumerate(tqdm(channels, desc=f"Page {page_idx + 1}/{n_pages}")):
    ax = axes_flat[i]

    jpss_bt = jpss_df.loc[jpss_df["channel_index"] == sc, "observation"].values
    ufs_bt = ufs_df.loc[ufs_df["channel_index"] == sc, "observation"].values

    # Compute shared bin range
    all_bt = np.concatenate([jpss_bt, ufs_bt])
    lo, hi = np.nanpercentile(all_bt, [0.5, 99.5])
    margin = max((hi - lo) * 0.05, 0.5)
    bins = np.linspace(lo - margin, hi + margin, HIST_BINS + 1)

    ax.hist(
        jpss_bt,
        bins=bins,
        alpha=HIST_ALPHA,
        color="steelblue",
        label=f"JPSS ({len(jpss_bt):,})",
        density=True,
    )
    ax.hist(
        ufs_bt,
        bins=bins,
        alpha=HIST_ALPHA,
        color="orangered",
        label=f"UFS ({len(ufs_bt):,})",
        density=True,
    )

    wn = _sensor_chan_to_wn(sc)
    band = _band_name(sc)
    ax.set_title(f"sc {sc}  ({band} {wn:.1f} cm$^{{-1}}$)", fontsize=8)
    ax.tick_params(labelsize=6)
    ax.legend(fontsize=5, loc="upper right")

    # Show mean delta in corner
    mean_jpss = np.nanmean(jpss_bt) if len(jpss_bt) > 0 else np.nan
    mean_ufs = np.nanmean(ufs_bt) if len(ufs_bt) > 0 else np.nan
    delta = mean_jpss - mean_ufs
    ax.text(
        0.03,
        0.95,
        f"$\\Delta$={delta:+.2f} K",
        transform=ax.transAxes,
        fontsize=6,
        verticalalignment="top",
        color="black",
        bbox=dict(facecolor="white", alpha=0.7, edgecolor="none", pad=1),
    )

# Hide unused subplots
for j in range(len(channels), GRID_ROWS * GRID_COLS):
    axes_flat[j].set_visible(False)

sc_lo = channels[0]
sc_hi = channels[-1]
fig.suptitle(
    f"CrIS FSR BT Distribution — JPSS vs UFS (n20, 1h tolerance)\n"
    f"sensor_chan {sc_lo}–{sc_hi}  |  Page {page_idx + 1}/{n_pages}",
    fontsize=14,
    fontweight="bold",
)

out_path = f"outputs/cris_bt_histograms_page{page_idx + 1:02d}.png"
fig.savefig(out_path, dpi=150, bbox_inches="tight")
logger.info(f"Saved {out_path}")

logger.info("Done.")

</details>

[greptile-apps]

Greptile Summary

This PR adds JPSS_CRIS, a new DataFrameSource for JPSS CrIS FSR Level-1 brightness temperatures, closely mirroring the existing JPSS_ATMS implementation. The channel-mapping, Hamming apodization, Planck inversion, and guard-channel handling are all carefully verified and well-tested; prior review concerns about satellite start dates, dedup keys, and QF metadata semantics have been addressed.

Confidence Score: 5/5

Safe to merge; all remaining findings are P2 with no runtime impact on the primary brightness-temperature data path.

Prior P0/P1 concerns (satellite start dates, dedup key, QF filtering semantics) are resolved. The three remaining findings are all P2: a wrong wavenumber range in a reference-only lexicon attribute, unused dict entries, and a nibble-overlap in the quality metadata field that only matters if QF3 values exceed 15.

earth2studio/lexicon/jpss.py (CRIS_BAND_RANGES upper bounds); earth2studio/data/jpss_cris.py (QF3 packing, _SDR_QF_KEYS dead entries)

Important Files Changed

Filename	Overview
earth2studio/data/jpss_cris.py	New 1291-line DataFrameSource for JPSS CrIS FSR L1 brightness temperatures; sound architecture mirroring JPSS_ATMS, with a few minor correctness nits in QF packing and dead code in _SDR_QF_KEYS
earth2studio/lexicon/jpss.py	Adds JPSSCrISLexicon; CRIS_BAND_RANGES upper bounds are off by 2 channels (states 4 high-end guards in comment, actual layout is 2)
test/data/test_jpss_cris.py	Comprehensive offline mock tests plus network slow-tests; covers apodized/unapodized paths, subsampling, schema field selection, validation, and filename parsing
earth2studio/data/init.py	One-line export of JPSS_CRIS — correct
earth2studio/lexicon/init.py	Exports JPSSCrISLexicon alongside existing JPSS lexicons — correct

<sub>Reviews (2): Last reviewed commit: "Merge branch 'main' into ngeneva/cris" | Re-trigger Greptile</sub>

[NickGeneva]

@greptile-apps

[NickGeneva]

/blossom-ci