Skip to content

Update prep_aria to support handling of NISAR tropo and SET layers - #1514

Open
sssangha wants to merge 2 commits into
mainfrom
nisar_tropo_set
Open

Update prep_aria to support handling of NISAR tropo and SET layers#1514
sssangha wants to merge 2 commits into
mainfrom
nisar_tropo_set

Conversation

@sssangha

@sssangha sssangha commented Jul 31, 2026

Copy link
Copy Markdown
Contributor

Description of proposed changes

Reminders

  • Fix #xxxx
  • Pass Pre-commit check (green)
  • Pass Codacy code review (green)
  • Pass Circle CI test (green)
  • Make sure that your code follows our style. Use the other functions/files as a basis.
  • If modifying functionality, describe changes to function behavior and arguments in a comment below the function declaration.
  • If adding new functionality, add a detailed description to the documentation and/or an example.

Summary by Sourcery

Add support in prep_aria for ingesting NISAR differential correction layers and converting them into epoch-wise timeseries datasets.

New Features:

  • Introduce inversion of NISAR differential troposphere and solid Earth tide layers into per-epoch displacement/path-delay timeseries.
  • Add extraction of unique NISAR acquisition dates and sensing times from VRT metadata when building correction timeseries.

Enhancements:

  • Extend correction layer type detection to recognize additional troposphere model names and generic tropo layers from metadata.
  • Update ARIA correction loading to choose between direct timeseries writing and differential inversion based on platform metadata and to standardize per-layer HDF5 output paths.

@sourcery-ai

sourcery-ai Bot commented Jul 31, 2026

Copy link
Copy Markdown
Contributor

Reviewer's Guide

Implements NISAR-specific handling of differential correction layers in prep_aria by inverting pairwise corrections into epoch-wise time series, extending tropo layer detection, and branching the ARIA loading path to handle NISAR metadata and per-layer HDF5 outputs.

Sequence diagram for NISAR-aware correction layer loading in load_aria

sequenceDiagram
    participant load_aria
    participant get_correction_layer
    participant get_number_of_epochs
    participant get_nisar_dates
    participant writefile_layout_hdf5 as writefile.layout_hdf5
    participant write_timeseries
    participant invert_diff_corrections

    load_aria->>get_correction_layer: get_correction_layer(layer)
    get_correction_layer-->>load_aria: layer_name, layer_type

    load_aria->>load_aria: is_nisar = meta[PLATFORM].upper().startswith(NISAR)

    alt is_nisar
        load_aria->>get_nisar_dates: get_nisar_dates(layer)
        get_nisar_dates-->>load_aria: date_list
        load_aria->>load_aria: num_dates = len(date_list)
    else not NISAR
        load_aria->>get_number_of_epochs: get_number_of_epochs(layer)
        get_number_of_epochs-->>load_aria: num_dates
    end

    load_aria->>load_aria: define ds_name_dict
    load_aria->>writefile_layout_hdf5: layout_hdf5(out_file, ds_name_dict, metadata)

    alt is_nisar
        load_aria->>invert_diff_corrections: invert_diff_corrections(out_file, corrStack=layer, box, xstep, ystep, mli_method)
        invert_diff_corrections-->>load_aria: out_file
    else not NISAR
        load_aria->>write_timeseries: write_timeseries(out_file, corrStack=layer, box, xstep, ystep)
        write_timeseries-->>load_aria: out_file
    end
Loading

File-Level Changes

Change Details Files
Add utilities to extract NISAR acquisition dates and invert differential correction stacks into epoch-wise time series with correct physical sign conventions.
  • Introduce get_nisar_dates() to parse band metadata, extract unique acquisition dates and UTC sensing times from NISAR VRTs, and return sorted date lists plus a date-to-UTC lookup.
  • Implement invert_diff_corrections() that reads differential correction pairs, converts phase to range using layer-specific sign conventions, builds and pseudo-inverts a design matrix to solve for per-epoch values, handles multilooking and NaN masking, and writes a timeseries HDF5 including date and sensingMid arrays and MODIFICATION_TIME metadata.
src/mintpy/prep_aria.py
Extend correction layer type detection and ARIA loading logic to recognize NISAR tropo/SET layers and route them through the new inversion path while standardizing timeseries HDF5 creation.
  • Augment get_correction_layer() to classify additional tropospheric model names and generic tropo-like layer names as tropo, using a new tropo_models list and a substring check.
  • Modify load_aria() to detect NISAR platforms via metadata, derive the number of epochs from NISAR differential layers using get_nisar_dates(), compute per-layer HDF5 output file paths, reuse ds_name_dict for layout creation, and choose between invert_diff_corrections() for NISAR layers and write_timeseries() for non-NISAR layers while aligning run_or_skip handling with the new out_file naming.
src/mintpy/prep_aria.py

Tips and commands

Interacting with Sourcery

  • Trigger a new review: Comment @sourcery-ai review on the pull request.
  • Continue discussions: Reply directly to Sourcery's review comments.
  • Generate a GitHub issue from a review comment: Ask Sourcery to create an
    issue from a review comment by replying to it. You can also reply to a
    review comment with @sourcery-ai issue to create an issue from it.
  • Generate a pull request title: Write @sourcery-ai anywhere in the pull
    request title to generate a title at any time. You can also comment
    @sourcery-ai title on the pull request to (re-)generate the title at any time.
  • Generate a pull request summary: Write @sourcery-ai summary anywhere in
    the pull request body to generate a PR summary at any time exactly where you
    want it. You can also comment @sourcery-ai summary on the pull request to
    (re-)generate the summary at any time.
  • Generate reviewer's guide: Comment @sourcery-ai guide on the pull
    request to (re-)generate the reviewer's guide at any time.
  • Resolve all Sourcery comments: Comment @sourcery-ai resolve on the
    pull request to resolve all Sourcery comments. Useful if you've already
    addressed all the comments and don't want to see them anymore.
  • Dismiss all Sourcery reviews: Comment @sourcery-ai dismiss on the pull
    request to dismiss all existing Sourcery reviews. Especially useful if you
    want to start fresh with a new review - don't forget to comment
    @sourcery-ai review to trigger a new review!

Customizing Your Experience

Access your dashboard to:

  • Enable or disable review features such as the Sourcery-generated pull request
    summary, the reviewer's guide, and others.
  • Change the review language.
  • Add, remove or edit custom review instructions.
  • Adjust other review settings.

Getting Help

@sssangha

Copy link
Copy Markdown
Contributor Author

To replicate my test case, access the latest dev version of ARIA-tools and run the following:

ariaDownload.py --mission NISAR -d 'descending' -o Url -s 20251030 -e 20251230 -b '28.80868 31.28686 -92.74528 -95.86692' -t 012

ariaTSsetup.py -f "products/*h5" -l 'troposphereTotal,ionosphere,solidEarthTide' -w TS_extract_masked -of ISCE -d Download -m Download >& ariaTSsetup.log

prep_aria.py -s TS_extract_masked/stack -d TS_extract_masked/DEM/glo_90.dem -w TS_extract_masked/mask/esa_world_cover_2021.msk -i TS_extract_masked/incidenceAngle/*.vrt -a TS_extract_masked/azimuthAngle/*.vrt --iono 'TS_extract_masked/stack/ionoStack.vrt' --set TS_extract_masked/stack/setStack.vrt --tropo TS_extract_masked/stack/tropoStack.vrt >& prep_aria.log

Here is the comparison of the MintPy internal pyaps workflow derived tropo delays (from ERA5) on the left vs the prep_aria.py derived tropo delays directly from the NISAR GUNW extracted layers:
Screenshot 2026-07-30 at 10 49 37 PM

Here is the comparison of the MintPy internal pysolid workflow derived SET on the left vs the prep_aria.py derived SET directly from the NISAR GUNW extracted layers:
Screenshot 2026-07-30 at 10 51 20 PM

@sssangha
sssangha marked this pull request as ready for review July 31, 2026 05:55

@sourcery-ai sourcery-ai Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hey - I've found 2 issues, and left some high level feedback:

  • Both get_nisar_dates() and invert_diff_corrections() assume specific metadata keys (Dates, UTCTime (HH:MM:SS.ss), Wavelength (m)) and a non-empty metadata domain list; add defensive checks and clearer error messages when these are missing or differently named to avoid obscure runtime failures on variant NISAR products.
  • In invert_diff_corrections(), the handling of NaNs replaces them with 0.0 before inversion and then zeroes the entire time series for pixels with any NaN; consider preserving a mask (e.g., using NaNs in the output) so that invalid pixels are distinguishable from valid zero-valued corrections.
Prompt for AI Agents
Please address the comments from this code review:

## Overall Comments
- Both `get_nisar_dates()` and `invert_diff_corrections()` assume specific metadata keys (`Dates`, `UTCTime (HH:MM:SS.ss)`, `Wavelength (m)`) and a non-empty metadata domain list; add defensive checks and clearer error messages when these are missing or differently named to avoid obscure runtime failures on variant NISAR products.
- In `invert_diff_corrections()`, the handling of NaNs replaces them with 0.0 before inversion and then zeroes the entire time series for pixels with any NaN; consider preserving a mask (e.g., using NaNs in the output) so that invalid pixels are distinguishable from valid zero-valued corrections.

## Individual Comments

### Comment 1
<location path="src/mintpy/prep_aria.py" line_range="689-693" />
<code_context>

     # Get type of correction
-    if layer_name in ['GMAO', 'HRES', 'HRRR', "ERA5"]:
+    tropo_models = [
+        'GMAO', 'HRES', 'HRRR', 'ERA5',
+        'troposphereTotal', 'troposphereWet', 'troposphereHydrostatic',
+    ]
+    if layer_name in tropo_models or 'trop' in layer_name.lower():
         layer_type = 'tropo'
     else:
</code_context>
<issue_to_address>
**question (bug_risk):** The generic `'trop'` substring match might over-classify some layers as tropospheric.

`if layer_name in tropo_models or 'trop' in layer_name.lower():` will mark any layer containing `trop` as `tropo`, even if it isn’t actually tropospheric. That could cause incorrect use of the `phase2range` sign convention. If only specific model names should be treated as troposphere, consider tightening the match (e.g., prefix/regex) or using just the explicit allowlist.
</issue_to_address>

### Comment 2
<location path="src/mintpy/prep_aria.py" line_range="552" />
<code_context>
+    return date_list, date_utc_dict
+
+
+def invert_diff_corrections(outfile, corrStack, box=None, xstep=1, ystep=1, mli_method='nearest'):
+    """Invert differential pairwise corrections into epoch-wise timeseries.
+
</code_context>
<issue_to_address>
**issue (complexity):** Consider refactoring the new NISAR path by splitting `invert_diff_corrections` into smaller helpers for math, HDF5 writing, and design-matrix creation, and by encapsulating the NISAR/non-NISAR branching logic in `load_aria` to keep the main flow simple.

The new NISAR logic is functionally solid but `invert_diff_corrections` is taking on too many responsibilities and duplicating parts of the existing timeseries writing path. You can reduce complexity by teasing apart the math, IO, and metadata while keeping behavior identical.

### 1. Extract math into a focused helper

Move the dense inversion + masking logic into a small helper that takes already‑prepared arrays. This keeps `invert_diff_corrections` focused on reading metadata and orchestrating:

```python
def compute_epoch_timeseries(diff_data, A, phase2range, no_data_val):
    """Invert differential pairs → epoch time series.

    diff_data: (n_pairs, length, width) phase corrections
    A:        (n_pairs, n_dates) design matrix
    """
    n_pairs, length, width = diff_data.shape
    n_dates = A.shape[1]
    num_pixels = length * width

    # Apply phase → range and mask NoData
    diff_data = diff_data.astype(np.float32)
    diff_data[diff_data == no_data_val] = np.nan
    diff_data *= phase2range

    A_sub   = A[:, 1:]
    A_pinv  = np.linalg.pinv(A_sub)
    flat    = diff_data.reshape(n_pairs, num_pixels)
    nan_mask = np.isnan(flat)
    flat[nan_mask] = 0.0

    ts_sub = A_pinv @ flat

    ts_flat = np.zeros((n_dates, num_pixels), dtype=np.float32)
    ts_flat[1:, :] = ts_sub

    # Any NaN in any pair ⇒ zero all epochs for that pixel
    has_nan_pixel = np.any(nan_mask, axis=0)
    ts_flat[:, has_nan_pixel] = 0.0

    return ts_flat.reshape(n_dates, length, width)
```

Then `invert_diff_corrections` can be simplified to:

```python
# after building A and diff_data (without phase2range / masking)
timeseries = compute_epoch_timeseries(
    diff_data=diff_data,
    A=A,
    phase2range=phase2range,
    no_data_val=no_data_val,
)
```

This isolates the inversion semantics (incl. NaN masking) in one place and makes them easier to test.

### 2. Extract a shared HDF5 writer for timeseries cubes

The HDF5 write path in `invert_diff_corrections` mirrors `write_timeseries`. You can factor a common helper that both NISAR and non‑NISAR flows call, avoiding divergence when metadata layout changes:

```python
def write_timeseries_cube(outfile, date_list, date_utc_dict, timeseries):
    """Write timeseries cube + date/sensingMid to HDF5."""
    with h5py.File(outfile, "a") as f:
        n_dates = len(date_list)
        prog_bar = ptime.progressBar(maxValue=n_dates)
        for ii, date in enumerate(date_list):
            utc = date_utc_dict.get(date, f"{date} 00:00:00")
            prog_bar.update(ii + 1, suffix=f'{date} {ii + 1}/{n_dates}')

            f["date"][ii]        = date.encode("utf-8")
            f["sensingMid"][ii]  = utc.encode("utf-8")
            f["timeseries"][ii]  = timeseries[ii]

        prog_bar.close()
        f["timeseries"].attrs['MODIFICATION_TIME'] = str(time.time())
```

Then in `invert_diff_corrections`:

```python
timeseries = compute_epoch_timeseries(...)
write_timeseries_cube(outfile, date_list, date_utc_dict, timeseries)
```

And you can similarly refactor `write_timeseries` to use `write_timeseries_cube` once it has built its `timeseries` array, so there is a single path for writing the HDF5 cube.

### 3. Separate design‑matrix construction from IO

Building `A` is currently interleaved with reading each band. You can simplify this by a small helper that only maps pair metadata → `A`:

```python
def build_pair_design_matrix(ds_cor, layer_name, date2idx):
    n_pairs = ds_cor.RasterCount
    n_dates = len(date2idx)
    A = np.zeros((n_pairs, n_dates), dtype=np.float32)

    for ii in range(n_pairs):
        bnd = ds_cor.GetRasterBand(ii + 1)
        d_sec, d_ref = bnd.GetMetadata(layer_name)["Dates"].split("_")
        A[ii, date2idx[d_sec]] =  1.0
        A[ii, date2idx[d_ref]] = -1.0

    return A
```

Then `invert_diff_corrections` becomes:

```python
date_list, date_utc_dict = get_nisar_dates(corrStack)
date2idx = {d: i for i, d in enumerate(date_list)}

A = build_pair_design_matrix(ds_cor, layer_name, date2idx)

# separate loop purely for reading data into diff_data
for ii in range(n_pairs):
    bnd = ds_cor.GetRasterBand(ii + 1)
    data = bnd.ReadAsArray(**kwargs)
    if xstep * ystep > 1:
        data = multilook_data(data, ystep, xstep, method=mli_method)
    diff_data[ii] = data
```

This reduces cognitive load by making “what indices correspond to which dates” an explicit, testable unit instead of being buried in the read loop.

### 4. Keep NISAR branching in `load_aria` small

The new NISAR branch in `load_aria` is reasonable, but you can keep it from growing by encapsulating the strategy selection:

```python
def get_correction_writer(layer, meta, inps):
    is_nisar = meta.get("PLATFORM", "").upper().startswith("NISAR")
    if is_nisar:
        return invert_diff_corrections
    else:
        return write_timeseries
```

Usage:

```python
writer = get_correction_writer(layer, meta, inps)
if run_or_skip(inps, ds_name_dict, out_file=out_file) == 'run':
    writefile.layout_hdf5(out_file, ds_name_dict, metadata=meta, ...)
    writer(out_file, corrStack=layer, box=box,
           xstep=inps.xstep, ystep=inps.ystep,
           mli_method=getattr(inps, "method", "nearest"))
```

This keeps platform‑specific branching out of the main `load_aria` loop and makes it easier to add future platform behaviors without expanding the loop itself.

These small, focused helpers keep all existing functionality (including the NISAR‑specific inversion semantics) but make the code easier to reason about, test, and evolve.
</issue_to_address>

Sourcery is free for open source - if you like our reviews please consider sharing them ✨
Help me be more useful! Please click 👍 or 👎 on each comment and I'll use the feedback to improve your reviews.

Comment thread src/mintpy/prep_aria.py Outdated
Comment thread src/mintpy/prep_aria.py
return date_list, date_utc_dict


def invert_diff_corrections(outfile, corrStack, box=None, xstep=1, ystep=1, mli_method='nearest'):

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

issue (complexity): Consider refactoring the new NISAR path by splitting invert_diff_corrections into smaller helpers for math, HDF5 writing, and design-matrix creation, and by encapsulating the NISAR/non-NISAR branching logic in load_aria to keep the main flow simple.

The new NISAR logic is functionally solid but invert_diff_corrections is taking on too many responsibilities and duplicating parts of the existing timeseries writing path. You can reduce complexity by teasing apart the math, IO, and metadata while keeping behavior identical.

1. Extract math into a focused helper

Move the dense inversion + masking logic into a small helper that takes already‑prepared arrays. This keeps invert_diff_corrections focused on reading metadata and orchestrating:

def compute_epoch_timeseries(diff_data, A, phase2range, no_data_val):
    """Invert differential pairs → epoch time series.

    diff_data: (n_pairs, length, width) phase corrections
    A:        (n_pairs, n_dates) design matrix
    """
    n_pairs, length, width = diff_data.shape
    n_dates = A.shape[1]
    num_pixels = length * width

    # Apply phase → range and mask NoData
    diff_data = diff_data.astype(np.float32)
    diff_data[diff_data == no_data_val] = np.nan
    diff_data *= phase2range

    A_sub   = A[:, 1:]
    A_pinv  = np.linalg.pinv(A_sub)
    flat    = diff_data.reshape(n_pairs, num_pixels)
    nan_mask = np.isnan(flat)
    flat[nan_mask] = 0.0

    ts_sub = A_pinv @ flat

    ts_flat = np.zeros((n_dates, num_pixels), dtype=np.float32)
    ts_flat[1:, :] = ts_sub

    # Any NaN in any pair ⇒ zero all epochs for that pixel
    has_nan_pixel = np.any(nan_mask, axis=0)
    ts_flat[:, has_nan_pixel] = 0.0

    return ts_flat.reshape(n_dates, length, width)

Then invert_diff_corrections can be simplified to:

# after building A and diff_data (without phase2range / masking)
timeseries = compute_epoch_timeseries(
    diff_data=diff_data,
    A=A,
    phase2range=phase2range,
    no_data_val=no_data_val,
)

This isolates the inversion semantics (incl. NaN masking) in one place and makes them easier to test.

2. Extract a shared HDF5 writer for timeseries cubes

The HDF5 write path in invert_diff_corrections mirrors write_timeseries. You can factor a common helper that both NISAR and non‑NISAR flows call, avoiding divergence when metadata layout changes:

def write_timeseries_cube(outfile, date_list, date_utc_dict, timeseries):
    """Write timeseries cube + date/sensingMid to HDF5."""
    with h5py.File(outfile, "a") as f:
        n_dates = len(date_list)
        prog_bar = ptime.progressBar(maxValue=n_dates)
        for ii, date in enumerate(date_list):
            utc = date_utc_dict.get(date, f"{date} 00:00:00")
            prog_bar.update(ii + 1, suffix=f'{date} {ii + 1}/{n_dates}')

            f["date"][ii]        = date.encode("utf-8")
            f["sensingMid"][ii]  = utc.encode("utf-8")
            f["timeseries"][ii]  = timeseries[ii]

        prog_bar.close()
        f["timeseries"].attrs['MODIFICATION_TIME'] = str(time.time())

Then in invert_diff_corrections:

timeseries = compute_epoch_timeseries(...)
write_timeseries_cube(outfile, date_list, date_utc_dict, timeseries)

And you can similarly refactor write_timeseries to use write_timeseries_cube once it has built its timeseries array, so there is a single path for writing the HDF5 cube.

3. Separate design‑matrix construction from IO

Building A is currently interleaved with reading each band. You can simplify this by a small helper that only maps pair metadata → A:

def build_pair_design_matrix(ds_cor, layer_name, date2idx):
    n_pairs = ds_cor.RasterCount
    n_dates = len(date2idx)
    A = np.zeros((n_pairs, n_dates), dtype=np.float32)

    for ii in range(n_pairs):
        bnd = ds_cor.GetRasterBand(ii + 1)
        d_sec, d_ref = bnd.GetMetadata(layer_name)["Dates"].split("_")
        A[ii, date2idx[d_sec]] =  1.0
        A[ii, date2idx[d_ref]] = -1.0

    return A

Then invert_diff_corrections becomes:

date_list, date_utc_dict = get_nisar_dates(corrStack)
date2idx = {d: i for i, d in enumerate(date_list)}

A = build_pair_design_matrix(ds_cor, layer_name, date2idx)

# separate loop purely for reading data into diff_data
for ii in range(n_pairs):
    bnd = ds_cor.GetRasterBand(ii + 1)
    data = bnd.ReadAsArray(**kwargs)
    if xstep * ystep > 1:
        data = multilook_data(data, ystep, xstep, method=mli_method)
    diff_data[ii] = data

This reduces cognitive load by making “what indices correspond to which dates” an explicit, testable unit instead of being buried in the read loop.

4. Keep NISAR branching in load_aria small

The new NISAR branch in load_aria is reasonable, but you can keep it from growing by encapsulating the strategy selection:

def get_correction_writer(layer, meta, inps):
    is_nisar = meta.get("PLATFORM", "").upper().startswith("NISAR")
    if is_nisar:
        return invert_diff_corrections
    else:
        return write_timeseries

Usage:

writer = get_correction_writer(layer, meta, inps)
if run_or_skip(inps, ds_name_dict, out_file=out_file) == 'run':
    writefile.layout_hdf5(out_file, ds_name_dict, metadata=meta, ...)
    writer(out_file, corrStack=layer, box=box,
           xstep=inps.xstep, ystep=inps.ystep,
           mli_method=getattr(inps, "method", "nearest"))

This keeps platform‑specific branching out of the main load_aria loop and makes it easier to add future platform behaviors without expanding the loop itself.

These small, focused helpers keep all existing functionality (including the NISAR‑specific inversion semantics) but make the code easier to reason about, test, and evolve.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant