Update prep_aria to support handling of NISAR tropo and SET layers - #1514
Update prep_aria to support handling of NISAR tropo and SET layers#1514sssangha wants to merge 2 commits into
Conversation
Reviewer's GuideImplements 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_ariasequenceDiagram
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
File-Level Changes
Tips and commandsInteracting with Sourcery
Customizing Your ExperienceAccess your dashboard to:
Getting Help
|
There was a problem hiding this comment.
Hey - I've found 2 issues, and left some high level feedback:
- Both
get_nisar_dates()andinvert_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>Help me be more useful! Please click 👍 or 👎 on each comment and I'll use the feedback to improve your reviews.
| return date_list, date_utc_dict | ||
|
|
||
|
|
||
| def invert_diff_corrections(outfile, corrStack, box=None, xstep=1, ystep=1, mli_method='nearest'): |
There was a problem hiding this comment.
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 AThen 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] = dataThis 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_timeseriesUsage:
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.


Description of proposed changes
Reminders
Summary by Sourcery
Add support in prep_aria for ingesting NISAR differential correction layers and converting them into epoch-wise timeseries datasets.
New Features:
Enhancements: