Add ISCE3 processor support for data loading (prep_isce3.py)#1510
Draft
tzl21 wants to merge 28 commits into
Draft
Add ISCE3 processor support for data loading (prep_isce3.py)#1510tzl21 wants to merge 28 commits into
prep_isce3.py)#1510tzl21 wants to merge 28 commits into
Conversation
… tif to gdal reader - prep_isce3: merge geometry without ref_int_file cropping first, compute multilook from full-res shape vs target shape, then crop. ALOOKS/RLOOKS now reflect actual multilook factors (2,4 not 14,26). - isce3_utils: add crop_geometry_to_reference() to clip merged geometry to reference interferogram extent via gdalwarp. - readfile: add isce3 to processor lists that route to read_gdal, avoiding raw binary read of GeoTIFF files.
…aseline - prep_isce3: read full-resolution pixel size (dx=5, dy=10) from static_layers.h5, then compute ALOOKS/RLOOKS as target pixel size ratio (20/10=2, 20/5=4). Merge geometry cropped to interferogram extent directly via gdalwarp -te/-tr. - isce3_utils: remove unused multilook/crop helper functions. - stackDict: get_perp_baseline uses .get() with default 0 to avoid KeyError when baseline metadata is missing from RSC.
…adata with fallbacks
…d Bperp average (m) formats
… (waterMask auto-align)
…rom missing conncomp data
Contributor
Reviewer's GuideAdds ISCE3/Dolphin processor support to MintPy by introducing a new prep_isce3 pipeline (CLI + core utilities) that discovers burst-level static_layers HDF5 files, wraps them as GDAL VRTs, merges geometry to the interferogram grid, generates metadata/.rsc files for stacks, and wires this into the existing load_data dispatch and reading logic for ISCE3 GeoTIFF interferograms. Sequence diagram for ISCE3 processor preparation and dispatchsequenceDiagram
actor User
participant load_data_prepare_metadata as load_data.prepare_metadata
participant prep_isce3_cli_main as mintpy.cli.prep_isce3.main
participant prep_isce3_core as prep_isce3.prep_isce3
participant isce3_utils_extract_merge_geometry as isce3_utils.extract_merge_geometry
participant isce3_utils_read_baseline as isce3_utils.read_baseline_timeseries_isce3
participant prepare_stack_isce3 as prep_isce3.prepare_stack_isce3
User->>load_data_prepare_metadata: set processor = "isce3"
load_data_prepare_metadata->>prep_isce3_cli_main: main(iargs)
prep_isce3_cli_main->>prep_isce3_core: prep_isce3(inps)
alt metaFile is auto or missing
prep_isce3_core->>isce3_utils_generate_burst_xml_from_static: generate_burst_xml_from_static(static_layers.h5, burst.xml)
prep_isce3_core->>isce3_utils_extract_isce3_metadata: extract_isce3_metadata(burst.xml)
else metaFile provided
prep_isce3_core->>isce3_utils_extract_isce3_metadata: extract_isce3_metadata(meta_file)
end
prep_isce3_core->>isce3_utils_extract_merge_geometry: extract_merge_geometry(geom_dir, out_dir, geom_types, ref_int_file)
isce3_utils_extract_merge_geometry-->>prep_isce3_core: merged geometry files & updated LENGTH/WIDTH
prep_isce3_core->>isce3_utils_read_baseline: read_baseline_timeseries_isce3(baseline_dir)
isce3_utils_read_baseline-->>prep_isce3_core: baseline_dict
loop for each obs_files pattern
prep_isce3_core->>prepare_stack_isce3: prepare_stack_isce3(obs_file, metadata, baseline_dict)
prepare_stack_isce3-->>prep_isce3_core: .rsc files per interferogram
end
prep_isce3_core-->>load_data_prepare_metadata: prepared geometry + stack metadata
Flow diagram for ISCE3 geometry extraction and merging pipelineflowchart TD
A[extract_merge_geometry] --> B[discover burst_ids in geom_dir and extra_dirs]
B --> C[merge_geometry_files]
C --> C1[build_vrt_from_h5 or extract_h5_geometry]
C1 --> C2[gdal.Warp merge per geometry type]
C2 --> C3[compute_incidence_angle from los_east/los_north]
C2 --> C4[compute_azimuth_angle from los_east/los_north]
C3 --> D[write merged incidenceAngle.tif]
C4 --> E[write merged azimuthAngle.tif]
C2 --> F[write merged height.tif and other geometry]
F --> G[update metadata LENGTH/WIDTH via readfile.read_attribute]
G --> H[return merged geometry dict to prepare_geometry_isce3]
File-Level Changes
Tips and commandsInteracting with Sourcery
Customizing Your ExperienceAccess your dashboard to:
Getting Help
|
![sourcery-ai[bot]](https://avatars.githubusercontent.com/in/48477?s=60&v=4){kind=small}
Contributor
There was a problem hiding this comment.
Hey - I've found 2 security issues, 4 other issues, and left some high level feedback:
Security issues:
- Detected subprocess function 'run' without a static string. If this data can be controlled by a malicious actor, it may be an instance of command injection. Audit the use of this call to ensure it is not controllable by an external resource. You may consider using 'shlex.escape()'. (link)
- The native Python
xmllibrary is vulnerable to XML External Entity (XXE) attacks. These attacks can leak confidential data and "XML bombs" can cause denial of service. Do not use this library to parse untrusted input. Instead the Python documentation recommends usingdefusedxml. (link)
General comments:
- In
stackDict._warp_water_maskand other places where temporary files are created, avoidtempfile.mktempand useNamedTemporaryFileormkstempto prevent race conditions and potential security issues. - The increase of retry attempts in
tropo_pyaps3.dload_grib_filesfrom 3 to 30 may lead to very long hangs in failure cases; consider making this configurable or adding a timeout/backoff strategy instead of a hard‑coded high retry count. - Several new utilities (e.g.,
_setup_gdal_proj_dataandread_isce3_geotiff) introduce global side effects (environment variable changes) and assumptions about GDAL availability; consider localizing these behaviors, handling missing GDAL binaries more gracefully, and avoiding hard‑coded dtype mappings (e.g., treating all integer GeoTIFF bands asint16).
Prompt for AI Agents
Please address the comments from this code review:
## Overall Comments
- In `stackDict._warp_water_mask` and other places where temporary files are created, avoid `tempfile.mktemp` and use `NamedTemporaryFile` or `mkstemp` to prevent race conditions and potential security issues.
- The increase of retry attempts in `tropo_pyaps3.dload_grib_files` from 3 to 30 may lead to very long hangs in failure cases; consider making this configurable or adding a timeout/backoff strategy instead of a hard‑coded high retry count.
- Several new utilities (e.g., `_setup_gdal_proj_data` and `read_isce3_geotiff`) introduce global side effects (environment variable changes) and assumptions about GDAL availability; consider localizing these behaviors, handling missing GDAL binaries more gracefully, and avoiding hard‑coded dtype mappings (e.g., treating all integer GeoTIFF bands as `int16`).
## Individual Comments
### Comment 1
<location path="src/mintpy/utils/readfile.py" line_range="1442" />
<code_context>
+ meta['DATA_TYPE'] = 'float32'
+ elif 'int' in dtype:
+ meta['DATA_TYPE'] = 'int16'
+ meta['NO_DATA_VALUE'] = str(band.GetNoDataValue())
+
+ # File type and processor
</code_context>
<issue_to_address>
**issue:** Guard against None when reading GDAL nodata, to avoid propagating 'None' as a string.
band.GetNoDataValue() can return None when nodata is unset. Passing this through str() stores the literal 'None', which will break callers expecting a numeric nodata (e.g., float(meta['NO_DATA_VALUE'])). Please guard for None and only set NO_DATA_VALUE when you have a numeric value, or map None to an empty string in line with the rest of the metadata handling.
</issue_to_address>
### Comment 2
<location path="src/mintpy/objects/stackDict.py" line_range="646-655" />
<code_context>
+ # not apply -s_srs correctly for files missing embedded SRS).
+ import tempfile
+ import subprocess
+ tmp_f = tempfile.mktemp(suffix='.tif')
+ cmd = [
+ 'gdalwarp', '-overwrite', '-q',
+ '-s_srs', src_srs, '-t_srs', ref_srs,
+ '-te', str(xmin), str(ymin), str(xmax), str(ymax),
+ '-tr', str(dx), str(dy),
+ '-r', 'near',
+ '-of', 'GTiff', '-co', 'COMPRESS=LZW',
+ src_file, tmp_f,
+ ]
+ subprocess.run(cmd, check=True)
+ tmp_ds = gdal.Open(tmp_f)
+ result = tmp_ds.GetRasterBand(1).ReadAsArray()
+ tmp_ds = None
+ src_ds = None
+ os.remove(tmp_f)
+ return result
+
</code_context>
<issue_to_address>
**🚨 issue (security):** Avoid tempfile.mktemp due to race conditions; use NamedTemporaryFile/mkstemp instead.
Because gdalwarp writes to this path, using tempfile.mktemp() exposes a race condition where another process could create or tamper with the file first. Instead, create the file securely with tempfile.NamedTemporaryFile(delete=False, suffix='.tif') or tempfile.mkstemp(), pass that filename to gdalwarp, and ensure you remove the file in a finally block so it’s cleaned up even if gdalwarp fails.
</issue_to_address>
### Comment 3
<location path="src/mintpy/utils/isce3_utils.py" line_range="267" />
<code_context>
+ return baseline_dict
+
+
+def extract_h5_geometry(
+ h5_file: Union[str, Path],
+ output_dir: Path,
</code_context>
<issue_to_address>
**issue (complexity):** Consider refactoring shared HDF5/geometry helpers, extracting source collection/warp/LOS post-processing from `merge_geometry_files`, and moving orbit/heading math into a dedicated utility to keep this module focused and easier to maintain.
The new module is taking on many responsibilities in one place and has duplicated geometry/HDF5 logic. You can keep all behavior but reduce complexity by factoring out shared helpers and isolating concerns.
### 1. Factor out shared HDF5 → geo metadata logic
`extract_h5_geometry` and `build_vrt_from_h5` both compute EPSG, geotransform, and nodata in similar ways. Extract that into small helpers so each path focuses on its unique behavior:
```python
def _read_projection_and_coords(data_group,
projection_name="projection",
x_coords_name="x_coordinates",
y_coords_name="y_coordinates"):
# EPSG
epsg = 4326
if projection_name in data_group:
proj_ds = data_group[projection_name]
if 'epsg_code' in proj_ds.attrs:
epsg = int(proj_ds.attrs['epsg_code'])
elif 'spatial_ref' in proj_ds.attrs:
wkt = proj_ds.attrs['spatial_ref'].decode('utf-8')
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
code = srs.GetAuthorityCode(None)
if code:
epsg = int(code)
elif proj_ds.shape == () and np.issubdtype(proj_ds.dtype, np.integer):
epsg = int(proj_ds[()])
# geotransform
x_coords = data_group[x_coords_name][:]
y_coords = data_group[y_coords_name][:]
dx = abs(x_coords[1] - x_coords[0])
dy = abs(y_coords[1] - y_coords[0])
left = x_coords[0] - dx / 2
top = y_coords[0] + dy / 2
geotransform = (left, dx, 0.0, top, 0.0, -dy)
return epsg, geotransform
```
```python
def _derive_nodata(dataset):
if '_FillValue' in dataset.attrs:
return dataset.attrs['_FillValue'].item()
if np.issubdtype(dataset.dtype, np.floating):
return np.nan
if np.issubdtype(dataset.dtype, np.integer):
return np.iinfo(dataset.dtype).max
return None
```
Then both `extract_h5_geometry` and `build_vrt_from_h5` call these helpers:
```python
with h5py.File(h5_file, 'r') as h5f:
data_group = h5f['/data']
epsg, geotransform = _read_projection_and_coords(data_group)
...
nodata = _derive_nodata(dataset)
```
This keeps behavior identical but removes duplication and divergent logic.
### 2. Simplify `merge_geometry_files` control flow
Right now `merge_geometry_files` owns:
- temp dir lifecycle
- per-burst collection
- VRT vs GeoTIFF fallback
- warp options/bounds/validation
- LOS post-processing (incidence/azimuth)
You can make it easier to scan by pushing “collect sources” into a separate function and keeping the VRT vs GeoTIFF decision outside:
```python
def collect_geometry_sources(
burst_ids: List[str],
geom_dir: Path,
temp_dir: Path,
geom_types: List[str],
use_vrt: bool,
) -> tuple[dict, dict]:
geometry_files = defaultdict(list)
nodata_dict = {}
for burst_id in burst_ids:
burst_path = geom_dir / burst_id
if not burst_path.exists():
continue
for h5_file in sorted(burst_path.glob("static_layers*.h5")):
extracted = (
build_vrt_from_h5(h5_file, temp_dir, geom_types)
if use_vrt
else extract_h5_geometry(h5_file, temp_dir, geom_types)
)
for gtype, info in extracted.items():
if info['file_list']:
geometry_files[gtype].append(info['file_list'])
if gtype not in nodata_dict and info['nodata'] is not None:
nodata_dict[gtype] = info['nodata']
return geometry_files, nodata_dict
```
```python
# in merge_geometry_files
try:
geometry_files, nodata_dict = collect_geometry_sources(
burst_ids, geom_dir, temp_dir, geom_types, use_vrt=True
)
except Exception as e:
print(f'WARNING: fast VRT-based merge failed: {e}')
print(' Falling back to legacy per-burst GeoTIFF extraction ...')
geometry_files, nodata_dict = collect_geometry_sources(
burst_ids, geom_dir, temp_dir, geom_types, use_vrt=False
)
```
You can similarly extract a `warp_geometry_type(...)` helper that encapsulates the `gdal.Warp` invocation + validation and a `postprocess_los_geometry(...)` that calls `compute_incidence_angle` and `compute_azimuth_angle`. That keeps `merge_geometry_files` as a short orchestration function.
### 3. Isolate orbit/heading math from metadata extraction
`_orbit_interp_hermite` and `_compute_heading` are quite heavy and are currently interwoven into `extract_required_attributes`. Moving them into a small, separate orbit utility makes `extract_required_attributes` easier to understand and test:
```python
# orbit_utils.py
def interpolate_state(metadata, sensing_mid_s):
t_array = np.array(metadata['time'])
pos_x = CubicHermiteSpline(t_array, metadata['position_x'], metadata['velocity_x'])
pos_y = CubicHermiteSpline(t_array, metadata['position_y'], metadata['velocity_y'])
pos_z = CubicHermiteSpline(t_array, metadata['position_z'], metadata['velocity_z'])
position = np.array([pos_x(sensing_mid_s), pos_y(sensing_mid_s), pos_z(sensing_mid_s)])
velocity = np.array([
pos_x.derivative()(sensing_mid_s),
pos_y.derivative()(sensing_mid_s),
pos_z.derivative()(sensing_mid_s),
])
return position, velocity
```
```python
# in extract_required_attributes
if has_orbit and isce3_available:
try:
sensingMid = _to_seconds(metadata['sensing_mid'], metadata['reference_epoch'])
position, velocity = interpolate_state(metadata, sensingMid)
ellipsoid = isce3.core.Ellipsoid()
llh = ellipsoid.xyz_to_lon_lat(position)
heading = _compute_heading((position, velocity))
...
except Exception as e:
...
```
This keeps all existing functionality but separates orbit math from general metadata munging, making each part more focused and easier to maintain.
If you later decide the Hermite implementation is overkill, you can swap `interpolate_state` to use `isce3.core.Orbit` without touching the higher-level metadata code.
</issue_to_address>
### Comment 4
<location path="src/mintpy/prep_isce3.py" line_range="30" />
<code_context>
+ return metadata
+
+
+def prepare_geometry_isce3(geom_dir, out_dir, geom_files=None, metadata=None,
+ processor='tops', update_mode=True, ref_int_file=None,
+ target_shape=None, geom_dirs=None):
</code_context>
<issue_to_address>
**issue (complexity):** Consider extracting small helper functions for geometry preparation, stack metadata handling, and meta-file generation so each function has a single, focused responsibility and is easier to understand and maintain.
You can keep all functionality and significantly reduce cognitive load by extracting a few focused helpers. Here are concrete, small changes that keep behavior identical:
---
### 1. Split `prepare_geometry_isce3` into helpers
The function currently mixes: burst/HDF5 discovery, pixel size inference, geometry merging, looks computation, `.rsc` writing, and metadata dimension updates.
You can tease out the burst + full‑res pixel size logic and the `.rsc` writing into helpers.
**a) Extract HDF5 discovery + pixel size logic**
```python
def _infer_full_res_pixel_size(geom_dir, geom_dirs=None):
"""Return (burst_full_dx, burst_full_dy, num_bursts, first_h5_path)."""
from pathlib import Path
import h5py
burst_full_dx = None
burst_full_dy = None
num_bursts = 0
first_h5_path = None
geom_path = Path(geom_dir)
burst_subdirs = sorted(
d for d in geom_path.iterdir()
if d.is_dir() and list(d.glob('static_layers*.h5'))
)
h5_files_root = sorted(geom_path.glob('static_layers*.h5'))
if burst_subdirs:
first_h5_path = sorted(burst_subdirs[0].glob('static_layers*.h5'))[0]
num_bursts = len(burst_subdirs)
elif h5_files_root:
first_h5_path = h5_files_root[0]
num_bursts = 1
if geom_dirs and len(geom_dirs) > 1:
for edir in geom_dirs[1:]:
epath = Path(edir)
if not epath.is_dir():
continue
if list(epath.glob('static_layers*.h5')):
num_bursts += 1
else:
for sub in epath.iterdir():
if sub.is_dir() and list(sub.glob('static_layers*.h5')):
num_bursts += 1
if first_h5_path:
with h5py.File(first_h5_path, 'r') as h5:
x_coords = h5['/data/x_coordinates'][:]
y_coords = h5['/data/y_coordinates'][:]
burst_full_dx = abs(x_coords[1] - x_coords[0])
burst_full_dy = abs(y_coords[1] - y_coords[0])
return burst_full_dx, burst_full_dy, num_bursts, first_h5_path
```
Then the top of `prepare_geometry_isce3` becomes:
```python
burst_full_dx, burst_full_dy, num_bursts, first_h5_path = _infer_full_res_pixel_size(
geom_dir, geom_dirs=geom_dirs
)
if first_h5_path:
print(f'Number of bursts: {num_bursts}')
print(f'Full-resolution pixel size: dx={burst_full_dx}, dy={burst_full_dy}')
else:
print('WARNING: could not read full-res pixel size')
```
This moves all the directory/shape handling into one place.
**b) Extract `.rsc` writing for geometry outputs**
```python
def _write_geometry_rsc_files(geometry_dict, metadata, update_mode=True):
"""Write .rsc files for all geometry products."""
for geom_name, geom_path in geometry_dict.items():
if not geom_path or not os.path.isfile(geom_path):
continue
geom_path = str(geom_path)
rsc_file = geom_path + '.rsc'
geom_meta = metadata.copy() if metadata else {}
geom_meta.update(readfile.read_attribute(geom_path))
writefile.write_roipac_rsc(
geom_meta, rsc_file, update_mode=update_mode, print_msg=True
)
```
Then replace your current `.rsc` loop with:
```python
_write_geometry_rsc_files(geometry_dict, metadata, update_mode=update_mode)
```
**c) Optional: extract looks computation**
```python
def _compute_looks_from_resolutions(burst_full_dx, burst_full_dy, ref_int_file):
"""Return (lks_y, lks_x) based on full-res vs interferogram pixel size."""
from osgeo import gdal
if not (burst_full_dx and burst_full_dy and ref_int_file):
return 1, 1
ref_ds = gdal.Open(str(ref_int_file))
if not ref_ds:
return 1, 1
ref_gt = ref_ds.GetGeoTransform()
ref_ds = None
ref_dx = abs(ref_gt[1])
ref_dy = abs(ref_gt[5])
lks_x = max(1, int(round(ref_dx / burst_full_dx)))
lks_y = max(1, int(round(ref_dy / burst_full_dy)))
return lks_y, lks_x
```
And in `prepare_geometry_isce3`:
```python
if metadata is not None:
lks_y, lks_x = _compute_looks_from_resolutions(
burst_full_dx, burst_full_dy, ref_int_file
)
metadata['ALOOKS'] = str(lks_y)
metadata['RLOOKS'] = str(lks_x)
print(f'Set ALOOKS={lks_y}, RLOOKS={lks_x} '
f'(full-res dx={burst_full_dx}, dy={burst_full_dy})')
```
---
### 2. Simplify date + metadata enrichment in `prepare_stack_isce3`
The loop currently parses dates, validates ranges, merges attributes, adds baseline, and writes `.rsc` inside a single `try` with a broad `except Exception`.
You can isolate filename parsing and interferogram metadata building.
**a) Extract filename date parsing**
```python
def _parse_dates_from_filename(fbase):
"""Return (d1, d2) as YYYYMMDD strings, or None if invalid."""
date_nums = re.findall(r'\d{8}', fbase)
if len(date_nums) < 2:
return None
dates = sorted(set(date_nums))[:2]
try:
if not all(19000101 <= int(d) <= 20991231 for d in dates):
return None
except ValueError:
return None
return dates[0], dates[1]
```
Use it in the loop:
```python
for i, tif_file in enumerate(tif_files):
fbase = os.path.basename(tif_file)
dates_pair = _parse_dates_from_filename(fbase)
if not dates_pair:
prog_bar.update(i+1, suffix=f'skipped {i+1}/{num_file}')
continue
d1, d2 = dates_pair
dates = [d1, d2]
# rest of logic...
```
**b) Extract interferogram metadata building**
```python
def _build_ifgram_metadata(base_meta, tif_file, dates, baseline_dict):
"""Return metadata dict for one interferogram."""
ifg_meta = base_meta.copy()
ifg_meta.update(readfile.read_attribute(tif_file))
ifg_meta = add_ifgram_metadata(ifg_meta, dates, baseline_dict)
return ifg_meta
```
And use inside the loop:
```python
try:
# dates_parsing...
num_valid += 1
prog_bar.update(i+1, suffix=f'{dates[0]}_{dates[1]} {i+1}/{num_file}')
ifg_meta = _build_ifgram_metadata(meta, tif_file, dates, baseline_dict)
rsc_file = tif_file + '.rsc'
writefile.write_roipac_rsc(
ifg_meta, rsc_file, update_mode=update_mode, print_msg=False
)
except Exception as e:
prog_bar.update(i+1, suffix=f'error {i+1}/{num_file} ({type(e).__name__})')
continue
```
This keeps the same behavior but separates concerns and logs exception types.
---
### 3. Extract “auto meta file” handling out of `prep_isce3`
The “generate burst XML from static_layers.h5” block can become a helper, so `prep_isce3` is mostly orchestration.
```python
def _ensure_meta_file_from_static_layers(inps):
"""Ensure inps.meta_file is set, generating XML from static_layers if needed."""
if inps.meta_file and inps.meta_file != 'auto':
return inps.meta_file
print('No meta file provided. Generating burst XML from static_layers.h5...')
geom_path = Path(inps.geom_dir)
burst_dirs = sorted(
d for d in geom_path.iterdir()
if d.is_dir() and list(d.glob('static_layers*.h5'))
)
h5_files_in_root = sorted(geom_path.glob('static_layers*.h5'))
if not burst_dirs and not h5_files_in_root:
raise FileNotFoundError(f'No static_layers HDF5 found in {inps.geom_dir}')
if burst_dirs:
first_burst_dir = burst_dirs[0]
h5_list = sorted(first_burst_dir.glob('static_layers*.h5'))
first_h5 = h5_list[0]
burst_id = first_burst_dir.name
else:
first_h5 = h5_files_in_root[0]
burst_id = os.path.basename(geom_path)
xml_out = Path(inps.out_dir) / f'{burst_id}.burst.xml'
xml_out.parent.mkdir(parents=True, exist_ok=True)
isce3_utils.generate_burst_xml_from_static(str(first_h5), str(xml_out))
inps.meta_file = str(xml_out)
print(f'Generated meta file: {inps.meta_file}')
return inps.meta_file
```
Then the top of `prep_isce3` becomes:
```python
def prep_isce3(inps):
"""Prepare ISCE3/Dolphin metadata files."""
_ensure_meta_file_from_static_layers(inps)
metadata = {}
if inps.meta_file:
metadata = isce3_utils.extract_isce3_metadata(
inps.meta_file, update_mode=inps.update_mode
)
# target_shape / ref_int logic...
# geometry preparation...
# baseline + stack preparation...
```
This keeps `prep_isce3` as a clearer sequence of steps.
---
These extractions are small, preserve all current behavior (including logging and optional geom_dirs), and make each function more single‑purpose, which should address the complexity concerns without undoing your additions.
</issue_to_address>
### Comment 5
<location path="src/mintpy/objects/stackDict.py" line_range="656" />
<code_context>
subprocess.run(cmd, check=True)
</code_context>
<issue_to_address>
**security (python.lang.security.audit.dangerous-subprocess-use-audit):** Detected subprocess function 'run' without a static string. If this data can be controlled by a malicious actor, it may be an instance of command injection. Audit the use of this call to ensure it is not controllable by an external resource. You may consider using 'shlex.escape()'.
*Source: opengrep*
</issue_to_address>
### Comment 6
<location path="src/mintpy/utils/isce3_utils.py" line_range="58" />
<code_context>
tree = ET.parse(meta_file)
</code_context>
<issue_to_address>
**security (python.lang.security.use-defused-xml-parse):** The native Python `xml` library is vulnerable to XML External Entity (XXE) attacks. These attacks can leak confidential data and "XML bombs" can cause denial of service. Do not use this library to parse untrusted input. Instead the Python documentation recommends using `defusedxml`.
```suggestion
tree = defusedxml.etree.ElementTree.parse(meta_file)
```
*Source: opengrep*
</issue_to_address>Help me be more useful! Please click 👍 or 👎 on each comment and I'll use the feedback to improve your reviews.
tzl21
marked this pull request as draft
July 22, 2026 01:40
tzl21
marked this pull request as ready for review
July 22, 2026 01:48
tzl21
marked this pull request as draft
July 22, 2026 01:48
tzl21
marked this pull request as ready for review
July 22, 2026 01:49
Contributor
There was a problem hiding this comment.
Hey - I've found 2 security issues, 6 other issues, and left some high level feedback:
Security issues:
- Detected subprocess function 'run' without a static string. If this data can be controlled by a malicious actor, it may be an instance of command injection. Audit the use of this call to ensure it is not controllable by an external resource. You may consider using 'shlex.escape()'. (link)
- The native Python
xmllibrary is vulnerable to XML External Entity (XXE) attacks. These attacks can leak confidential data and "XML bombs" can cause denial of service. Do not use this library to parse untrusted input. Instead the Python documentation recommends usingdefusedxml. (link)
General comments:
- In stackDict._warp_water_mask(), if gdalwarp fails or subprocess.run raises,
resultis never set but is still returned after the finally block, which will cause a NameError; consider initializingresultbefore the try and only returning it on success (or re-raising on failure). - isce3_utils._setup_gdal_proj_data() modifies PROJ/GDAL-related environment variables at import time, which may be surprising for callers; you might want to gate this behind a function call or only run it lazily when a GDAL/OSR operation actually fails.
- The new _warp_water_mask() shell-out to the gdalwarp CLI assumes it is on PATH and does not surface stderr on failure; consider using the GDAL Python API (gdal.Warp) or at least capturing and logging stderr to make diagnosing CRS/reprojection issues easier.
Prompt for AI Agents
Please address the comments from this code review:
## Overall Comments
- In stackDict._warp_water_mask(), if gdalwarp fails or subprocess.run raises, `result` is never set but is still returned after the finally block, which will cause a NameError; consider initializing `result` before the try and only returning it on success (or re-raising on failure).
- isce3_utils._setup_gdal_proj_data() modifies PROJ/GDAL-related environment variables at import time, which may be surprising for callers; you might want to gate this behind a function call or only run it lazily when a GDAL/OSR operation actually fails.
- The new _warp_water_mask() shell-out to the gdalwarp CLI assumes it is on PATH and does not surface stderr on failure; consider using the GDAL Python API (gdal.Warp) or at least capturing and logging stderr to make diagnosing CRS/reprojection issues easier.
## Individual Comments
### Comment 1
<location path="src/mintpy/objects/stackDict.py" line_range="612-621" />
<code_context>
+ def _warp_water_mask(self, dsName, target_length, target_width):
</code_context>
<issue_to_address>
**issue (bug_risk):** Handle failures in gdalwarp path more robustly and avoid returning an uninitialized result.
In `_warp_water_mask`, if `subprocess.run(cmd, check=True)` or `gdal.Open(tmp_f)` fails, `result` is never assigned, yet the function still does `return result`, causing an `UnboundLocalError`. Also, `ref_ds = gdal.Open(ref_file)` and `src_ds = gdal.Open(src_file)` may be `None`, so later calls to `GetGeoTransform()` / `GetProjection()` can raise attribute errors. Please add explicit checks that all `gdal.Open` calls succeed and raise a clear error otherwise, handle `subprocess.CalledProcessError` (and GDAL errors) explicitly, and only return `result` from inside the block where it is guaranteed to be set rather than at the end of the function.
</issue_to_address>
### Comment 2
<location path="src/mintpy/utils/readfile.py" line_range="1435-1442" />
<code_context>
+ meta['EPSG'] = epsg
+
+ # Data type
+ band = ds.GetRasterBand(1)
+ dtype = gdal.GetDataTypeName(band.DataType).lower()
+ if 'float' in dtype:
+ meta['DATA_TYPE'] = 'float32'
+ elif 'int' in dtype:
+ meta['DATA_TYPE'] = 'int16'
+ ndv = band.GetNoDataValue()
</code_context>
<issue_to_address>
**suggestion (bug_risk):** Refine data-type mapping for ISCE3 GeoTIFFs to avoid misclassifying integer bit depths.
The current `dtype` check (`'float' in dtype` → `float32`, `'int' in dtype` → `int16`) is too coarse and will mislabel types like `Int32`, `UInt16`, or `Byte` as `int16`, which can break downstream tools that depend on correct bit depth and scaling. Please map GDAL types explicitly (e.g., via `band.DataType` or the full `dtype` string) so that `int16`, `int32`, and `uint*` are distinguished, and choose a conservative default for unknown types.
```suggestion
# Data type
band = ds.GetRasterBand(1)
gdal_type = band.DataType
gdal_type_name = gdal.GetDataTypeName(gdal_type).lower()
data_type_map = {
gdal.GDT_Byte: 'uint8',
gdal.GDT_UInt16: 'uint16',
gdal.GDT_Int16: 'int16',
gdal.GDT_UInt32: 'uint32',
gdal.GDT_Int32: 'int32',
gdal.GDT_Float32: 'float32',
gdal.GDT_Float64: 'float64',
}
# Prefer explicit mapping; fall back to GDAL type name for unknown/complex types
meta['DATA_TYPE'] = data_type_map.get(gdal_type, gdal_type_name)
ndv = band.GetNoDataValue()
```
</issue_to_address>
### Comment 3
<location path="src/mintpy/utils/isce3_utils.py" line_range="692-389" />
<code_context>
+def compute_azimuth_angle(los_east_file: Path, los_north_file: Path, output_file: Path, nodata: float = None):
</code_context>
<issue_to_address>
**suggestion:** Mask nodata using the explicit nodata value rather than only NaNs when computing azimuth angle.
In `compute_azimuth_angle`, the nodata mask only checks `np.isnan(east) | np.isnan(north)`, so rasters that use a finite numeric nodata (e.g., -9999) will be treated as valid and included in the azimuth calculation. Please also mask on the explicit `nodata` value when it is finite (e.g., `(east == nodata) | (north == nodata)` in addition to the NaN check), and apply the same approach in `compute_incidence_angle`.
</issue_to_address>
### Comment 4
<location path="src/mintpy/utils/isce3_utils.py" line_range="735-389" />
<code_context>
+def compute_incidence_angle(los_east_file: Path, los_north_file: Path, output_file: Path, nodata: float = None):
</code_context>
<issue_to_address>
**suggestion:** Apply nodata masking that matches the actual nodata encoding of LOS rasters in incidence-angle computation.
`compute_incidence_angle` currently only masks NaN values for `east`/`north`. When LOS rasters use a finite nodata value, that flag will propagate into `up_sq` and `inc_angle` and yield bogus angles instead of nodata. Please extend the mask to also cover `(east == nodata) | (north == nodata)` when `nodata` is finite, consistent with `merge_geometry_files`.
Suggested implementation:
```python
def compute_incidence_angle(los_east_file: Path, los_north_file: Path, output_file: Path, nodata: float = None):
"""Compute incidence angle from LOS east and north components.
Parameters
----------
los_east_file, los_north_file : Path
Paths to GeoTIFF files of LOS vector components.
output_file : Path
Output path for incidenceAngle.tif.
nodata : float, optional
No-data value to use.
```
```python
# Mask invalid LOS samples. When a finite nodata flag is used in the LOS rasters,
# treat those pixels as nodata in the incidence-angle computation as well.
if nodata is None or np.isnan(nodata):
mask = np.isnan(east) | np.isnan(north)
else:
mask = (east == nodata) | (north == nodata) | np.isnan(east) | np.isnan(north)
```
I only see the function signature and docstring, so I assumed there is an existing line that defines `mask` as `np.isnan(east) | np.isnan(north)` inside `compute_incidence_angle`. If the variable names or masking logic differ, update the `SEARCH`/`REPLACE` block accordingly:
1. Locate where the current NaN-only mask is computed for `east`/`north` inside `compute_incidence_angle`.
2. Replace that mask expression with the conditional logic shown above so that when `nodata` is a finite value, pixels equal to `nodata` in either LOS raster are also masked.
3. Ensure that wherever the mask is later applied to `inc_angle` (e.g., `inc_angle[mask] = nodata` or similar), the same `nodata` value is written as the band nodata using `band.SetNoDataValue(nodata)`, consistent with `merge_geometry_files`.
</issue_to_address>
### Comment 5
<location path="src/mintpy/utils/isce3_utils.py" line_range="267" />
<code_context>
+ return baseline_dict
+
+
+def extract_h5_geometry(
+ h5_file: Union[str, Path],
+ output_dir: Path,
</code_context>
<issue_to_address>
**issue (complexity):** Consider refactoring shared HDF5 geometry handling and LOS GDAL IO into small helper functions so the main routines focus on their core logic and are easier to follow.
You’ve packed a lot of logic into this module; one concrete, low‑risk way to reduce complexity (without changing behavior or reverting features) is to factor out some of the duplicated patterns so the “big” functions are easier to follow.
### 1. Factor common HDF5 geometry logic
`extract_h5_geometry` and `build_vrt_from_h5` both:
- define the same default `dataset_mapping`,
- derive EPSG and geotransform from the same `/data` group,
- compute nodata in the same way.
You can pull that into small internal helpers so both functions become narrower “orchestrators” rather than each re‑implementing the same details.
For example:
```python
# shared defaults and helpers near the top of the file
_DEFAULT_GEOM_DATASET_MAPPING = {
"height.tif": "z",
"layover_shadow_mask.tif": "layover_shadow_mask",
"local_incidence_angle.tif": "local_incidence_angle",
"los_east.tif": "los_east",
"los_north.tif": "los_north",
}
def _get_epsg_and_geotransform(data_group,
x_coords_name="x_coordinates",
y_coords_name="y_coordinates",
projection_name="projection"):
epsg = 4326
if projection_name in data_group:
proj_ds = data_group[projection_name]
if "epsg_code" in proj_ds.attrs:
epsg = int(proj_ds.attrs["epsg_code"])
elif proj_ds.shape == () and np.issubdtype(proj_ds.dtype, np.integer):
epsg = int(proj_ds[()])
x_coords = data_group[x_coords_name][:] if x_coords_name in data_group else None
y_coords = data_group[y_coords_name][:] if y_coords_name in data_group else None
geotransform = None
if x_coords is not None and y_coords is not None and len(x_coords) > 1 and len(y_coords) > 1:
dx = abs(x_coords[1] - x_coords[0])
dy = abs(y_coords[1] - y_coords[0])
left = x_coords[0] - dx / 2
top = y_coords[0] + dy / 2
geotransform = (left, dx, 0.0, top, 0.0, -dy)
return epsg, geotransform
def _get_nodata_for_dataset(dataset):
if "_FillValue" in dataset.attrs:
return dataset.attrs["_FillValue"].item()
if np.issubdtype(dataset.dtype, np.floating):
return np.nan
if np.issubdtype(dataset.dtype, np.integer):
return np.iinfo(dataset.dtype).max
return None
```
Then `extract_h5_geometry` and `build_vrt_from_h5` can reuse these:
```python
def extract_h5_geometry(..., dataset_mapping: Optional[Dict[str, str]] = None, ...):
if dataset_mapping is None:
dataset_mapping = _DEFAULT_GEOM_DATASET_MAPPING
extracted = defaultdict(lambda: {"file_list": None, "nodata": None})
h5_file = Path(h5_file)
if not h5_file.exists():
return extracted
try:
with h5py.File(h5_file, "r") as h5f:
data_group = h5f["/data"]
epsg, geotransform = _get_epsg_and_geotransform(
data_group, x_coords_name, y_coords_name, projection_name
)
...
for geom_type in geom_types:
ds_path = dataset_mapping.get(geom_type)
...
nodata = _get_nodata_for_dataset(dataset)
...
```
```python
def build_vrt_from_h5(..., dataset_mapping: Optional[Dict[str, str]] = None, ...):
if dataset_mapping is None:
dataset_mapping = _DEFAULT_GEOM_DATASET_MAPPING
extracted = defaultdict(lambda: {"file_list": None, "nodata": None})
h5_file = Path(h5_file).resolve()
if not h5_file.exists():
return extracted
with h5py.File(h5_file, "r") as h5f:
data_group = h5f["/data"]
epsg, geotransform = _get_epsg_and_geotransform(data_group)
...
for geom_type in geom_types:
ds_name = dataset_mapping.get(geom_type)
...
nodata = _get_nodata_for_dataset(dataset)
...
```
This keeps all current behavior (including fallbacks and defaults) but makes it much clearer what each function is doing, and any future bugfixes to EPSG/geotransform/nodata happen in one place.
### 2. Reduce duplication in LOS angle computations
`compute_incidence_angle` and `compute_azimuth_angle` duplicate the same GDAL IO boilerplate. A small pair of helpers can simplify both and make their math stand out:
```python
def _read_los_components(los_east_file: Path, los_north_file: Path):
ds_east = gdal.Open(str(los_east_file))
ds_north = gdal.Open(str(los_north_file))
east = ds_east.GetRasterBand(1).ReadAsArray()
north = ds_north.GetRasterBand(1).ReadAsArray()
return ds_east, ds_north, east, north
def _write_single_band_tif(template_ds, data: np.ndarray, output_file: Path, nodata: float = None):
driver = gdal.GetDriverByName("GTiff")
ds_out = driver.Create(
str(output_file),
template_ds.RasterXSize,
template_ds.RasterYSize,
1,
gdal.GDT_Float32,
options=["COMPRESS=LZW"],
)
ds_out.SetGeoTransform(template_ds.GetGeoTransform())
ds_out.SetProjection(template_ds.GetProjection())
band = ds_out.GetRasterBand(1)
band.WriteArray(data)
if nodata is not None:
band.SetNoDataValue(nodata)
band.FlushCache()
ds_out = None
```
Then:
```python
def compute_azimuth_angle(...):
ds_east, ds_north, east, north = _read_los_components(los_east_file, los_north_file)
az_angle = -1 * np.rad2deg(np.arctan2(east, north)) % 360.0
if nodata is not None:
mask = np.isnan(east) | np.isnan(north)
az_angle[mask] = nodata
_write_single_band_tif(ds_east, az_angle, output_file, nodata=nodata)
```
```python
def compute_incidence_angle(...):
ds_east, ds_north, east, north = _read_los_components(los_east_file, los_north_file)
up_sq = 1.0 - east**2 - north**2
up_sq = np.clip(up_sq, 0, 1)
up = np.sqrt(up_sq)
inc_angle = np.rad2deg(np.arccos(up))
if nodata is not None:
mask = np.isnan(east) | np.isnan(north)
inc_angle[mask] = nodata
_write_single_band_tif(ds_east, inc_angle, output_file, nodata=nodata)
```
This doesn’t change the mathematical logic or file formats, but it shortens each public function and clearly separates “what” (angle formula) from “how” (GDAL IO), making the module easier to reason about.
</issue_to_address>
### Comment 6
<location path="src/mintpy/prep_isce3.py" line_range="200" />
<code_context>
+
+
+#########################################################################
+def prep_isce3(inps):
+ """Prepare ISCE3/Dolphin metadata files."""
+ # If no meta file is provided or it is 'auto', generate one from static_layers
</code_context>
<issue_to_address>
**issue (complexity):** Consider extracting helpers for meta-file generation, target-shape detection, and static_layers discovery to simplify prep_isce3’s control flow and remove duplicated logic.
You can reduce the control‑flow / coupling in `prep_isce3` without changing behavior by extracting a couple of focused helpers and reusing the “find static_layers/burst info” logic.
### 1. Extract meta‑file generation into a helper
The “auto meta file” branch at the top of `prep_isce3` is a big chunk of logic that obscures the main orchestration. Pull it out into a dedicated helper that returns the resolved meta file path:
```python
def _ensure_meta_file(inps) -> str:
"""Return path to meta_file; auto-generate XML from static_layers if needed."""
if inps.meta_file and inps.meta_file != "auto":
return inps.meta_file
print("No meta file provided. Generating burst XML from static_layers.h5...")
geom_path = Path(inps.geom_dir)
burst_dirs = sorted(
d for d in geom_path.iterdir()
if d.is_dir() and list(d.glob("static_layers*.h5"))
)
h5_files_in_root = sorted(geom_path.glob("static_layers*.h5"))
if not burst_dirs and not h5_files_in_root:
raise FileNotFoundError(f"No static_layers HDF5 found in {inps.geom_dir}")
if burst_dirs:
first_burst_dir = burst_dirs[0]
first_h5 = sorted(first_burst_dir.glob("static_layers*.h5"))[0]
burst_id = first_burst_dir.name
else:
first_h5 = h5_files_in_root[0]
burst_id = os.path.basename(geom_path)
xml_out = Path(inps.out_dir) / f"{burst_id}.burst.xml"
xml_out.parent.mkdir(parents=True, exist_ok=True)
isce3_utils.generate_burst_xml_from_static(str(first_h5), str(xml_out))
print(f"Generated meta file: {xml_out}")
return str(xml_out)
```
Then `prep_isce3` becomes:
```python
def prep_isce3(inps):
"""Prepare ISCE3/Dolphin metadata files."""
inps.meta_file = _ensure_meta_file(inps)
metadata = isce3_utils.extract_isce3_metadata(
inps.meta_file,
update_mode=inps.update_mode,
)
target_shape, ref_int = _determine_target_shape(inps.obs_files)
if inps.geom_dir:
metadata = prepare_geometry_isce3(
geom_dir=inps.geom_dir,
out_dir=inps.out_dir,
geom_files=inps.geom_files,
metadata=metadata,
processor=inps.processor,
update_mode=inps.update_mode,
ref_int_file=ref_int,
target_shape=target_shape,
geom_dirs=getattr(inps, "geom_dirs", None),
)
baseline_dict = {}
if inps.baseline_dir:
baseline_dict = isce3_utils.read_baseline_timeseries_isce3(
inps.baseline_dir,
processor=inps.processor,
)
if inps.obs_files:
for obs_file in inps.obs_files:
prepare_stack_isce3(
obs_file,
metadata=metadata,
baseline_dict=baseline_dict,
update_mode=inps.update_mode,
)
print("Done.")
```
### 2. Extract target‑shape detection
Similarly, the “first interferogram” logic can be isolated and unit‑tested more easily:
```python
def _determine_target_shape(obs_files):
"""Return (target_shape, ref_int_file) from the first obs_files pattern, if any."""
if not obs_files:
return None, None
int_list = glob.glob(obs_files[0])
if not int_list:
return None, None
ref_int = int_list[0]
int_atr = readfile.read_attribute(ref_int)
target_shape = (int(int_atr["LENGTH"]), int(int_atr["WIDTH"]))
print(f"Target shape from interferogram: {target_shape}")
return target_shape, ref_int
```
### 3. Reuse “static_layers discovery” logic
The “find first static_layers HDF5 + burst id” pattern appears both in `prep_isce3` and in `prepare_geometry_isce3` (where you compute `first_h5_path`, `num_bursts`, etc.). You can centralize this into a single helper (in this module or in `isce3_utils`) and call it from both places to reduce duplication and make the burst‑handling semantics explicit.
For example:
```python
def find_static_layers(geom_dir, extra_dirs=None):
"""Return (first_h5_path, num_bursts)."""
geom_path = Path(geom_dir)
burst_subdirs = sorted(
d for d in geom_path.iterdir()
if d.is_dir() and list(d.glob("static_layers*.h5"))
)
h5_files_root = sorted(geom_path.glob("static_layers*.h5"))
num_bursts = len(burst_subdirs) if burst_subdirs else (1 if h5_files_root else 0)
if extra_dirs:
for edir in extra_dirs:
epath = Path(edir)
# ... same counting logic here ...
first_h5_path = (
sorted(burst_subdirs[0].glob("static_layers*.h5"))[0]
if burst_subdirs else (h5_files_root[0] if h5_files_root else None)
)
return first_h5_path, num_bursts
```
Then `prepare_geometry_isce3`’s “Step 0” becomes:
```python
first_h5_path, num_bursts = find_static_layers(geom_dir, extra_dirs=geom_dirs[1:])
if first_h5_path:
# existing h5py logic...
```
And `_ensure_meta_file` can reuse the same helper when determining `first_h5` and `burst_id`.
These small extractions keep all behavior intact but make the orchestration in `prep_isce3` much easier to scan and test.
</issue_to_address>
### Comment 7
<location path="src/mintpy/objects/stackDict.py" line_range="659" />
<code_context>
subprocess.run(cmd, check=True)
</code_context>
<issue_to_address>
**security (python.lang.security.audit.dangerous-subprocess-use-audit):** Detected subprocess function 'run' without a static string. If this data can be controlled by a malicious actor, it may be an instance of command injection. Audit the use of this call to ensure it is not controllable by an external resource. You may consider using 'shlex.escape()'.
*Source: opengrep*
</issue_to_address>
### Comment 8
<location path="src/mintpy/utils/isce3_utils.py" line_range="58" />
<code_context>
tree = ET.parse(meta_file)
</code_context>
<issue_to_address>
**security (python.lang.security.use-defused-xml-parse):** The native Python `xml` library is vulnerable to XML External Entity (XXE) attacks. These attacks can leak confidential data and "XML bombs" can cause denial of service. Do not use this library to parse untrusted input. Instead the Python documentation recommends using `defusedxml`.
```suggestion
tree = defusedxml.etree.ElementTree.parse(meta_file)
```
*Source: opengrep*
</issue_to_address>Help me be more useful! Please click 👍 or 👎 on each comment and I'll use the feedback to improve your reviews.
Comment on lines
+612
to
+621
| def _warp_water_mask(self, dsName, target_length, target_width): | ||
| """Reproject water mask to match reference geometry grid via GDAL Warp. | ||
|
|
||
| Parameters: dsName - str, dataset name (waterMask) | ||
| target_length - int, target rows | ||
| target_width - int, target columns | ||
| Returns: data - np.ndarray (target_length, target_width) | ||
| """ | ||
| src_file = self.datasetDict[dsName] | ||
| ref_file = self.datasetDict[self.dsNames[0]] |
Contributor
There was a problem hiding this comment.
issue (bug_risk): Handle failures in gdalwarp path more robustly and avoid returning an uninitialized result.
In _warp_water_mask, if subprocess.run(cmd, check=True) or gdal.Open(tmp_f) fails, result is never assigned, yet the function still does return result, causing an UnboundLocalError. Also, ref_ds = gdal.Open(ref_file) and src_ds = gdal.Open(src_file) may be None, so later calls to GetGeoTransform() / GetProjection() can raise attribute errors. Please add explicit checks that all gdal.Open calls succeed and raise a clear error otherwise, handle subprocess.CalledProcessError (and GDAL errors) explicitly, and only return result from inside the block where it is guaranteed to be set rather than at the end of the function.
| if nodata is not None: | ||
| band.SetNoDataValue(nodata) | ||
| band.FlushCache() | ||
| ds_out = None |
Contributor
There was a problem hiding this comment.
suggestion: Mask nodata using the explicit nodata value rather than only NaNs when computing azimuth angle.
In compute_azimuth_angle, the nodata mask only checks np.isnan(east) | np.isnan(north), so rasters that use a finite numeric nodata (e.g., -9999) will be treated as valid and included in the azimuth calculation. Please also mask on the explicit nodata value when it is finite (e.g., (east == nodata) | (north == nodata) in addition to the NaN check), and apply the same approach in compute_incidence_angle.
| if nodata is not None: | ||
| band.SetNoDataValue(nodata) | ||
| band.FlushCache() | ||
| ds_out = None |
Contributor
There was a problem hiding this comment.
suggestion: Apply nodata masking that matches the actual nodata encoding of LOS rasters in incidence-angle computation.
compute_incidence_angle currently only masks NaN values for east/north. When LOS rasters use a finite nodata value, that flag will propagate into up_sq and inc_angle and yield bogus angles instead of nodata. Please extend the mask to also cover (east == nodata) | (north == nodata) when nodata is finite, consistent with merge_geometry_files.
Suggested implementation:
def compute_incidence_angle(los_east_file: Path, los_north_file: Path, output_file: Path, nodata: float = None):
"""Compute incidence angle from LOS east and north components.
Parameters
----------
los_east_file, los_north_file : Path
Paths to GeoTIFF files of LOS vector components.
output_file : Path
Output path for incidenceAngle.tif.
nodata : float, optional
No-data value to use. # Mask invalid LOS samples. When a finite nodata flag is used in the LOS rasters,
# treat those pixels as nodata in the incidence-angle computation as well.
if nodata is None or np.isnan(nodata):
mask = np.isnan(east) | np.isnan(north)
else:
mask = (east == nodata) | (north == nodata) | np.isnan(east) | np.isnan(north)I only see the function signature and docstring, so I assumed there is an existing line that defines mask as np.isnan(east) | np.isnan(north) inside compute_incidence_angle. If the variable names or masking logic differ, update the SEARCH/REPLACE block accordingly:
- Locate where the current NaN-only mask is computed for
east/northinsidecompute_incidence_angle. - Replace that mask expression with the conditional logic shown above so that when
nodatais a finite value, pixels equal tonodatain either LOS raster are also masked. - Ensure that wherever the mask is later applied to
inc_angle(e.g.,inc_angle[mask] = nodataor similar), the samenodatavalue is written as the band nodata usingband.SetNoDataValue(nodata), consistent withmerge_geometry_files.
| return baseline_dict | ||
|
|
||
|
|
||
| def extract_h5_geometry( |
Contributor
There was a problem hiding this comment.
issue (complexity): Consider refactoring shared HDF5 geometry handling and LOS GDAL IO into small helper functions so the main routines focus on their core logic and are easier to follow.
You’ve packed a lot of logic into this module; one concrete, low‑risk way to reduce complexity (without changing behavior or reverting features) is to factor out some of the duplicated patterns so the “big” functions are easier to follow.
1. Factor common HDF5 geometry logic
extract_h5_geometry and build_vrt_from_h5 both:
- define the same default
dataset_mapping, - derive EPSG and geotransform from the same
/datagroup, - compute nodata in the same way.
You can pull that into small internal helpers so both functions become narrower “orchestrators” rather than each re‑implementing the same details.
For example:
# shared defaults and helpers near the top of the file
_DEFAULT_GEOM_DATASET_MAPPING = {
"height.tif": "z",
"layover_shadow_mask.tif": "layover_shadow_mask",
"local_incidence_angle.tif": "local_incidence_angle",
"los_east.tif": "los_east",
"los_north.tif": "los_north",
}
def _get_epsg_and_geotransform(data_group,
x_coords_name="x_coordinates",
y_coords_name="y_coordinates",
projection_name="projection"):
epsg = 4326
if projection_name in data_group:
proj_ds = data_group[projection_name]
if "epsg_code" in proj_ds.attrs:
epsg = int(proj_ds.attrs["epsg_code"])
elif proj_ds.shape == () and np.issubdtype(proj_ds.dtype, np.integer):
epsg = int(proj_ds[()])
x_coords = data_group[x_coords_name][:] if x_coords_name in data_group else None
y_coords = data_group[y_coords_name][:] if y_coords_name in data_group else None
geotransform = None
if x_coords is not None and y_coords is not None and len(x_coords) > 1 and len(y_coords) > 1:
dx = abs(x_coords[1] - x_coords[0])
dy = abs(y_coords[1] - y_coords[0])
left = x_coords[0] - dx / 2
top = y_coords[0] + dy / 2
geotransform = (left, dx, 0.0, top, 0.0, -dy)
return epsg, geotransform
def _get_nodata_for_dataset(dataset):
if "_FillValue" in dataset.attrs:
return dataset.attrs["_FillValue"].item()
if np.issubdtype(dataset.dtype, np.floating):
return np.nan
if np.issubdtype(dataset.dtype, np.integer):
return np.iinfo(dataset.dtype).max
return NoneThen extract_h5_geometry and build_vrt_from_h5 can reuse these:
def extract_h5_geometry(..., dataset_mapping: Optional[Dict[str, str]] = None, ...):
if dataset_mapping is None:
dataset_mapping = _DEFAULT_GEOM_DATASET_MAPPING
extracted = defaultdict(lambda: {"file_list": None, "nodata": None})
h5_file = Path(h5_file)
if not h5_file.exists():
return extracted
try:
with h5py.File(h5_file, "r") as h5f:
data_group = h5f["/data"]
epsg, geotransform = _get_epsg_and_geotransform(
data_group, x_coords_name, y_coords_name, projection_name
)
...
for geom_type in geom_types:
ds_path = dataset_mapping.get(geom_type)
...
nodata = _get_nodata_for_dataset(dataset)
...def build_vrt_from_h5(..., dataset_mapping: Optional[Dict[str, str]] = None, ...):
if dataset_mapping is None:
dataset_mapping = _DEFAULT_GEOM_DATASET_MAPPING
extracted = defaultdict(lambda: {"file_list": None, "nodata": None})
h5_file = Path(h5_file).resolve()
if not h5_file.exists():
return extracted
with h5py.File(h5_file, "r") as h5f:
data_group = h5f["/data"]
epsg, geotransform = _get_epsg_and_geotransform(data_group)
...
for geom_type in geom_types:
ds_name = dataset_mapping.get(geom_type)
...
nodata = _get_nodata_for_dataset(dataset)
...This keeps all current behavior (including fallbacks and defaults) but makes it much clearer what each function is doing, and any future bugfixes to EPSG/geotransform/nodata happen in one place.
2. Reduce duplication in LOS angle computations
compute_incidence_angle and compute_azimuth_angle duplicate the same GDAL IO boilerplate. A small pair of helpers can simplify both and make their math stand out:
def _read_los_components(los_east_file: Path, los_north_file: Path):
ds_east = gdal.Open(str(los_east_file))
ds_north = gdal.Open(str(los_north_file))
east = ds_east.GetRasterBand(1).ReadAsArray()
north = ds_north.GetRasterBand(1).ReadAsArray()
return ds_east, ds_north, east, north
def _write_single_band_tif(template_ds, data: np.ndarray, output_file: Path, nodata: float = None):
driver = gdal.GetDriverByName("GTiff")
ds_out = driver.Create(
str(output_file),
template_ds.RasterXSize,
template_ds.RasterYSize,
1,
gdal.GDT_Float32,
options=["COMPRESS=LZW"],
)
ds_out.SetGeoTransform(template_ds.GetGeoTransform())
ds_out.SetProjection(template_ds.GetProjection())
band = ds_out.GetRasterBand(1)
band.WriteArray(data)
if nodata is not None:
band.SetNoDataValue(nodata)
band.FlushCache()
ds_out = NoneThen:
def compute_azimuth_angle(...):
ds_east, ds_north, east, north = _read_los_components(los_east_file, los_north_file)
az_angle = -1 * np.rad2deg(np.arctan2(east, north)) % 360.0
if nodata is not None:
mask = np.isnan(east) | np.isnan(north)
az_angle[mask] = nodata
_write_single_band_tif(ds_east, az_angle, output_file, nodata=nodata)def compute_incidence_angle(...):
ds_east, ds_north, east, north = _read_los_components(los_east_file, los_north_file)
up_sq = 1.0 - east**2 - north**2
up_sq = np.clip(up_sq, 0, 1)
up = np.sqrt(up_sq)
inc_angle = np.rad2deg(np.arccos(up))
if nodata is not None:
mask = np.isnan(east) | np.isnan(north)
inc_angle[mask] = nodata
_write_single_band_tif(ds_east, inc_angle, output_file, nodata=nodata)This doesn’t change the mathematical logic or file formats, but it shortens each public function and clearly separates “what” (angle formula) from “how” (GDAL IO), making the module easier to reason about.
|
|
||
|
|
||
| ######################################################################### | ||
| def prep_isce3(inps): |
Contributor
There was a problem hiding this comment.
issue (complexity): Consider extracting helpers for meta-file generation, target-shape detection, and static_layers discovery to simplify prep_isce3’s control flow and remove duplicated logic.
You can reduce the control‑flow / coupling in prep_isce3 without changing behavior by extracting a couple of focused helpers and reusing the “find static_layers/burst info” logic.
1. Extract meta‑file generation into a helper
The “auto meta file” branch at the top of prep_isce3 is a big chunk of logic that obscures the main orchestration. Pull it out into a dedicated helper that returns the resolved meta file path:
def _ensure_meta_file(inps) -> str:
"""Return path to meta_file; auto-generate XML from static_layers if needed."""
if inps.meta_file and inps.meta_file != "auto":
return inps.meta_file
print("No meta file provided. Generating burst XML from static_layers.h5...")
geom_path = Path(inps.geom_dir)
burst_dirs = sorted(
d for d in geom_path.iterdir()
if d.is_dir() and list(d.glob("static_layers*.h5"))
)
h5_files_in_root = sorted(geom_path.glob("static_layers*.h5"))
if not burst_dirs and not h5_files_in_root:
raise FileNotFoundError(f"No static_layers HDF5 found in {inps.geom_dir}")
if burst_dirs:
first_burst_dir = burst_dirs[0]
first_h5 = sorted(first_burst_dir.glob("static_layers*.h5"))[0]
burst_id = first_burst_dir.name
else:
first_h5 = h5_files_in_root[0]
burst_id = os.path.basename(geom_path)
xml_out = Path(inps.out_dir) / f"{burst_id}.burst.xml"
xml_out.parent.mkdir(parents=True, exist_ok=True)
isce3_utils.generate_burst_xml_from_static(str(first_h5), str(xml_out))
print(f"Generated meta file: {xml_out}")
return str(xml_out)Then prep_isce3 becomes:
def prep_isce3(inps):
"""Prepare ISCE3/Dolphin metadata files."""
inps.meta_file = _ensure_meta_file(inps)
metadata = isce3_utils.extract_isce3_metadata(
inps.meta_file,
update_mode=inps.update_mode,
)
target_shape, ref_int = _determine_target_shape(inps.obs_files)
if inps.geom_dir:
metadata = prepare_geometry_isce3(
geom_dir=inps.geom_dir,
out_dir=inps.out_dir,
geom_files=inps.geom_files,
metadata=metadata,
processor=inps.processor,
update_mode=inps.update_mode,
ref_int_file=ref_int,
target_shape=target_shape,
geom_dirs=getattr(inps, "geom_dirs", None),
)
baseline_dict = {}
if inps.baseline_dir:
baseline_dict = isce3_utils.read_baseline_timeseries_isce3(
inps.baseline_dir,
processor=inps.processor,
)
if inps.obs_files:
for obs_file in inps.obs_files:
prepare_stack_isce3(
obs_file,
metadata=metadata,
baseline_dict=baseline_dict,
update_mode=inps.update_mode,
)
print("Done.")2. Extract target‑shape detection
Similarly, the “first interferogram” logic can be isolated and unit‑tested more easily:
def _determine_target_shape(obs_files):
"""Return (target_shape, ref_int_file) from the first obs_files pattern, if any."""
if not obs_files:
return None, None
int_list = glob.glob(obs_files[0])
if not int_list:
return None, None
ref_int = int_list[0]
int_atr = readfile.read_attribute(ref_int)
target_shape = (int(int_atr["LENGTH"]), int(int_atr["WIDTH"]))
print(f"Target shape from interferogram: {target_shape}")
return target_shape, ref_int3. Reuse “static_layers discovery” logic
The “find first static_layers HDF5 + burst id” pattern appears both in prep_isce3 and in prepare_geometry_isce3 (where you compute first_h5_path, num_bursts, etc.). You can centralize this into a single helper (in this module or in isce3_utils) and call it from both places to reduce duplication and make the burst‑handling semantics explicit.
For example:
def find_static_layers(geom_dir, extra_dirs=None):
"""Return (first_h5_path, num_bursts)."""
geom_path = Path(geom_dir)
burst_subdirs = sorted(
d for d in geom_path.iterdir()
if d.is_dir() and list(d.glob("static_layers*.h5"))
)
h5_files_root = sorted(geom_path.glob("static_layers*.h5"))
num_bursts = len(burst_subdirs) if burst_subdirs else (1 if h5_files_root else 0)
if extra_dirs:
for edir in extra_dirs:
epath = Path(edir)
# ... same counting logic here ...
first_h5_path = (
sorted(burst_subdirs[0].glob("static_layers*.h5"))[0]
if burst_subdirs else (h5_files_root[0] if h5_files_root else None)
)
return first_h5_path, num_burstsThen prepare_geometry_isce3’s “Step 0” becomes:
first_h5_path, num_bursts = find_static_layers(geom_dir, extra_dirs=geom_dirs[1:])
if first_h5_path:
# existing h5py logic...And _ensure_meta_file can reuse the same helper when determining first_h5 and burst_id.
These small extractions keep all behavior intact but make the orchestration in prep_isce3 much easier to scan and test.
| '-of', 'GTiff', '-co', 'COMPRESS=LZW', | ||
| src_file, tmp_f, | ||
| ] | ||
| subprocess.run(cmd, check=True) |
Contributor
There was a problem hiding this comment.
security (python.lang.security.audit.dangerous-subprocess-use-audit): Detected subprocess function 'run' without a static string. If this data can be controlled by a malicious actor, it may be an instance of command injection. Audit the use of this call to ensure it is not controllable by an external resource. You may consider using 'shlex.escape()'.
Source: opengrep
| Common metadata dictionary with keys required by MintPy. | ||
| """ | ||
| # Parse XML file | ||
| tree = ET.parse(meta_file) |
Contributor
There was a problem hiding this comment.
security (python.lang.security.use-defused-xml-parse): The native Python xml library is vulnerable to XML External Entity (XXE) attacks. These attacks can leak confidential data and "XML bombs" can cause denial of service. Do not use this library to parse untrusted input. Instead the Python documentation recommends using defusedxml.
Suggested change
| tree = ET.parse(meta_file) | |
| tree = defusedxml.etree.ElementTree.parse(meta_file) |
Source: opengrep
Contributor
There was a problem hiding this comment.
Hey - I've found 2 security issues, 4 other issues, and left some high level feedback:
Security issues:
- Detected subprocess function 'run' without a static string. If this data can be controlled by a malicious actor, it may be an instance of command injection. Audit the use of this call to ensure it is not controllable by an external resource. You may consider using 'shlex.escape()'. (link)
- The native Python
xmllibrary is vulnerable to XML External Entity (XXE) attacks. These attacks can leak confidential data and "XML bombs" can cause denial of service. Do not use this library to parse untrusted input. Instead the Python documentation recommends usingdefusedxml. (link)
General comments:
- In
stackDict._warp_water_mask,resultis only defined inside thetryblock; ifgdalwarpfails (e.g.,subprocess.runraises), the function will hit the finalreturn resultwithresultundefined—consider returningNoneor raising on failure instead of referencing an uninitialized variable. - The Hermite interpolation path in
_orbit_interp_hermiteassumesCubicHermiteSplineis available, but when SciPy is missing it's set toNoneand still used; it would be safer to explicitly fall back to a simpler interpolation (or raise a clear error) whenCubicHermiteSplineis not importable. - Multiple modules now import
osgeo.gdal/osrat the top level (e.g.,readfile.py,stackDict.py,isce3_utils.py), which makes GDAL a hard runtime dependency for any MintPy usage; to avoid surprising failures and reduce startup cost, consider moving these imports inside the specific functions that need them or guarding them with clear error messages when missing.
Prompt for AI Agents
Please address the comments from this code review:
## Overall Comments
- In `stackDict._warp_water_mask`, `result` is only defined inside the `try` block; if `gdalwarp` fails (e.g., `subprocess.run` raises), the function will hit the final `return result` with `result` undefined—consider returning `None` or raising on failure instead of referencing an uninitialized variable.
- The Hermite interpolation path in `_orbit_interp_hermite` assumes `CubicHermiteSpline` is available, but when SciPy is missing it's set to `None` and still used; it would be safer to explicitly fall back to a simpler interpolation (or raise a clear error) when `CubicHermiteSpline` is not importable.
- Multiple modules now import `osgeo.gdal`/`osr` at the top level (e.g., `readfile.py`, `stackDict.py`, `isce3_utils.py`), which makes GDAL a hard runtime dependency for any MintPy usage; to avoid surprising failures and reduce startup cost, consider moving these imports inside the specific functions that need them or guarding them with clear error messages when missing.
## Individual Comments
### Comment 1
<location path="src/mintpy/objects/stackDict.py" line_range="649-658" />
<code_context>
+ import subprocess
+ fd, tmp_f = tempfile.mkstemp(suffix='.tif')
+ os.close(fd)
+ try:
+ cmd = [
+ 'gdalwarp', '-overwrite', '-q',
+ '-s_srs', src_srs, '-t_srs', ref_srs,
+ '-te', str(xmin), str(ymin), str(xmax), str(ymax),
+ '-tr', str(dx), str(dy),
+ '-r', 'near',
+ '-of', 'GTiff', '-co', 'COMPRESS=LZW',
+ src_file, tmp_f,
+ ]
+ subprocess.run(cmd, check=True)
+ tmp_ds = gdal.Open(tmp_f)
+ result = tmp_ds.GetRasterBand(1).ReadAsArray()
+ tmp_ds = None
+ src_ds = None
+ return result
+ finally:
+ if os.path.exists(tmp_f):
</code_context>
<issue_to_address>
**issue (bug_risk):** Avoid returning an undefined `result` when gdalwarp or GDAL open fails in `_warp_water_mask`.
If `subprocess.run` or `gdal.Open(tmp_f)` raises, the `return` in the `try` block is skipped, `finally` runs, and the trailing `return result` executes with `result` never assigned, causing an `UnboundLocalError`. Initialize `result = None` before the `try` and only return it if set, or re-raise/handle the original exception instead of returning `result` on failure paths.
</issue_to_address>
### Comment 2
<location path="src/mintpy/utils/isce3_utils.py" line_range="716-719" />
<code_context>
+ # az = -arctan2(east, north) * 180/pi (mod 360)
+ az_angle = -1 * np.rad2deg(np.arctan2(east, north)) % 360.0
+
+ # Mask nodata
+ if nodata is not None:
+ mask = np.isnan(east) | np.isnan(north)
+ az_angle[mask] = nodata
+
+ driver = gdal.GetDriverByName('GTiff')
</code_context>
<issue_to_address>
**issue (bug_risk):** Nodata masking in `compute_azimuth_angle` ignores non-NaN nodata values and may write angles over nodata pixels.
Currently nodata is only detected via `np.isnan`, so finite nodata sentinels from the GeoTIFF (e.g., large float/int values) will be treated as valid data and preserved in the azimuth result. Please also mask `east == nodata` and `north == nodata` when `nodata` is finite, consistent with how `dstNodata` is handled in the warp step.
</issue_to_address>
### Comment 3
<location path="src/mintpy/prep_isce3.py" line_range="30" />
<code_context>
+ return metadata
+
+
+def prepare_geometry_isce3(geom_dir, out_dir, geom_files=None, metadata=None,
+ processor='tops', update_mode=True, ref_int_file=None,
+ target_shape=None, geom_dirs=None):
</code_context>
<issue_to_address>
**issue (complexity):** Consider extracting small helper functions for discovery, parsing, and setup logic so the main orchestration functions remain linear and easier to read without changing behavior.
You can keep the current behavior and reduce orchestration complexity by extracting a few tight helpers. That will make `prepare_geometry_isce3`, `prepare_stack_isce3`, and `prep_isce3` easier to follow without changing logic.
### 1. Split geometry orchestration (`prepare_geometry_isce3`)
The early burst / HDF5 discovery and full-res pixel-size reading can be moved into a helper. That removes nested branching from the main flow:
```python
def _discover_static_layers(geom_dir: str, geom_dirs=None):
geom_path = Path(geom_dir)
burst_subdirs = sorted(
d for d in geom_path.iterdir()
if d.is_dir() and list(d.glob('static_layers*.h5'))
)
h5_files_root = sorted(geom_path.glob('static_layers*.h5'))
first_h5_path = None
if burst_subdirs:
first_h5_path = sorted(burst_subdirs[0].glob('static_layers*.h5'))[0]
num_bursts = len(burst_subdirs)
elif h5_files_root:
first_h5_path = h5_files_root[0]
num_bursts = 1
else:
num_bursts = 0
if geom_dirs and len(geom_dirs) > 1:
for edir in geom_dirs[1:]:
epath = Path(edir)
if not epath.is_dir():
continue
if list(epath.glob('static_layers*.h5')):
num_bursts += 1
else:
for sub in epath.iterdir():
if sub.is_dir() and list(sub.glob('static_layers*.h5')):
num_bursts += 1
return first_h5_path, num_bursts
```
```python
def _read_full_res_pixel_size(first_h5_path):
import h5py
with h5py.File(first_h5_path, 'r') as h5:
x_coords = h5['/data/x_coordinates'][:]
y_coords = h5['/data/y_coordinates'][:]
return abs(x_coords[1] - x_coords[0]), abs(y_coords[1] - y_coords[0])
```
Then `prepare_geometry_isce3` becomes more linear:
```python
def prepare_geometry_isce3(...):
...
first_h5_path, num_bursts = _discover_static_layers(geom_dir, geom_dirs)
burst_full_dx = burst_full_dy = None
if first_h5_path:
try:
burst_full_dx, burst_full_dy = _read_full_res_pixel_size(first_h5_path)
print(f'Number of bursts: {num_bursts}')
print(f'Full-resolution pixel size: dx={burst_full_dx}, dy={burst_full_dy}')
except Exception as e:
print(f'WARNING: could not read full-res pixel size: {e}')
geometry_dict = isce3_utils.extract_merge_geometry(...)
if metadata is not None:
_populate_looks_from_ref(metadata, ref_int_file, burst_full_dx, burst_full_dy)
_write_geometry_rsc_files(geometry_dict, metadata, update_mode)
_update_metadata_dims_from_height(geometry_dict, metadata)
return metadata
```
Where the other helpers encapsulate separate responsibilities:
```python
def _populate_looks_from_ref(metadata, ref_int_file, burst_full_dx, burst_full_dy):
lks_y = lks_x = 1
if burst_full_dx is not None and burst_full_dy is not None and ref_int_file:
from osgeo import gdal
ref_ds = gdal.Open(str(ref_int_file))
if ref_ds:
ref_gt = ref_ds.GetGeoTransform()
ref_dx = abs(ref_gt[1])
ref_dy = abs(ref_gt[5])
ref_ds = None
lks_x = max(1, int(round(ref_dx / burst_full_dx)))
lks_y = max(1, int(round(ref_dy / burst_full_dy)))
metadata['ALOOKS'] = str(lks_y)
metadata['RLOOKS'] = str(lks_x)
print(f'Set ALOOKS={lks_y}, RLOOKS={lks_x} '
f'(full-res dx={burst_full_dx}, dy={burst_full_dy})')
```
```python
def _write_geometry_rsc_files(geometry_dict, metadata, update_mode):
for _, geom_path in geometry_dict.items():
if geom_path and os.path.isfile(geom_path):
geom_path = str(geom_path)
rsc_file = geom_path + '.rsc'
geom_meta = metadata.copy() if metadata else {}
geom_meta.update(readfile.read_attribute(geom_path))
writefile.write_roipac_rsc(
geom_meta, rsc_file,
update_mode=update_mode,
print_msg=True,
)
```
```python
def _update_metadata_dims_from_height(geometry_dict, metadata):
if metadata is None:
return
height_file = geometry_dict.get('height.tif')
if height_file and os.path.isfile(height_file):
geom_atr = readfile.read_attribute(str(height_file))
metadata['LENGTH'] = geom_atr['LENGTH']
metadata['WIDTH'] = geom_atr['WIDTH']
print(f"Updated metadata with LENGTH={metadata['LENGTH']}, "
f"WIDTH={metadata['WIDTH']}")
```
This keeps all behavior but flattens `prepare_geometry_isce3` into clearly named steps.
### 2. Isolate date parsing / validation in `prepare_stack_isce3`
Right now file globbing, date parsing, metadata merge, and `.rsc` writing are interleaved, plus a broad `try/except`. Pulling date parsing into a helper makes the loop more readable and testable:
```python
def _extract_valid_dates_from_name(filename: str):
date_nums = re.findall(r'\d{8}', filename)
if len(date_nums) < 2:
return None
dates = sorted(set(date_nums))[:2]
if not all(19000101 <= int(d) <= 20991231 for d in dates):
return None
return dates
```
```python
def prepare_stack_isce3(...):
...
for i, tif_file in enumerate(tif_files):
fbase = os.path.basename(tif_file)
dates = _extract_valid_dates_from_name(fbase)
if dates is None:
prog_bar.update(i+1, suffix=f'skipped {i+1}/{num_file}')
continue
try:
num_valid += 1
prog_bar.update(i+1, suffix=f'{dates[0]}_{dates[1]} {i+1}/{num_file}')
ifg_meta = meta.copy()
ifg_meta.update(readfile.read_attribute(tif_file))
ifg_meta = add_ifgram_metadata(ifg_meta, dates, baseline_dict)
rsc_file = tif_file + '.rsc'
writefile.write_roipac_rsc(
ifg_meta, rsc_file,
update_mode=update_mode,
print_msg=False,
)
except Exception:
prog_bar.update(i+1, suffix=f'error {i+1}/{num_file}')
continue
```
This keeps the `try/except` only around I/O/metadata writing (where you actually expect failures) and makes the “enumerate valid observation files” step explicit.
### 3. Factor meta-file auto-generation out of `prep_isce3`
The auto-meta logic is useful but makes `prep_isce3` dense. Moving it to a helper clarifies the main orchestration:
```python
def _ensure_meta_file(inps):
if inps.meta_file and inps.meta_file != 'auto':
return inps.meta_file
print('No meta file provided. Generating burst XML from static_layers.h5...')
geom_path = Path(inps.geom_dir)
burst_dirs = sorted(
d for d in geom_path.iterdir()
if d.is_dir() and list(d.glob('static_layers*.h5'))
)
h5_files_in_root = sorted(geom_path.glob('static_layers*.h5'))
if not burst_dirs and not h5_files_in_root:
raise FileNotFoundError(f'No static_layers HDF5 found in {inps.geom_dir}')
if burst_dirs:
first_burst_dir = burst_dirs[0]
first_h5 = sorted(first_burst_dir.glob('static_layers*.h5'))[0]
burst_id = first_burst_dir.name
else:
first_h5 = h5_files_in_root[0]
burst_id = os.path.basename(geom_path)
xml_out = Path(inps.out_dir) / f'{burst_id}.burst.xml'
xml_out.parent.mkdir(parents=True, exist_ok=True)
isce3_utils.generate_burst_xml_from_static(str(first_h5), str(xml_out))
print(f'Generated meta file: {xml_out}')
return str(xml_out)
```
```python
def prep_isce3(inps):
meta_file = _ensure_meta_file(inps)
metadata = isce3_utils.extract_isce3_metadata(
meta_file,
update_mode=inps.update_mode,
)
ref_int, target_shape = _detect_target_shape(inps.obs_files)
if inps.geom_dir:
metadata = prepare_geometry_isce3(
geom_dir=inps.geom_dir,
out_dir=inps.out_dir,
geom_files=inps.geom_files,
metadata=metadata,
processor=inps.processor,
update_mode=inps.update_mode,
ref_int_file=ref_int,
target_shape=target_shape,
geom_dirs=getattr(inps, 'geom_dirs', None),
)
...
```
With a small helper for the target shape detection:
```python
def _detect_target_shape(obs_files):
if not obs_files:
return None, None
int_list = glob.glob(obs_files[0])
if not int_list:
return None, None
ref_int = int_list[0]
int_atr = readfile.read_attribute(ref_int)
target_shape = (int(int_atr['LENGTH']), int(int_atr['WIDTH']))
print(f'Target shape from interferogram: {target_shape}')
return ref_int, target_shape
```
These small extra functions reduce the cognitive load in the main orchestration functions while preserving all current behavior and I/O patterns.
</issue_to_address>
### Comment 4
<location path="src/mintpy/objects/stackDict.py" line_range="612" />
<code_context>
return self.metadata
+ def _warp_water_mask(self, dsName, target_length, target_width):
+ """Reproject water mask to match reference geometry grid via GDAL Warp.
+
</code_context>
<issue_to_address>
**issue (complexity):** Consider refactoring the water mask warping logic to use an in-memory `gdal.Warp` and a separate helper for reference grid computation to simplify `_warp_water_mask` and remove boilerplate.
You can keep the auto-align behavior while reducing complexity and removing the external `gdalwarp` subprocess by:
1. Using `gdal.Warp` directly with an in-memory dataset.
2. Splitting “compute reference grid” from “warp water mask”.
That keeps the functionality but makes `_warp_water_mask` smaller, easier to test, and removes tempfile / subprocess boilerplate.
### 1. Extract reference grid computation
```python
def _get_ref_grid(self, target_length, target_width):
"""Get reference grid bounds/resolution/SRS from the first dataset."""
ref_file = self.datasetDict[self.dsNames[0]]
ref_ds = gdal.Open(ref_file)
ref_gt = ref_ds.GetGeoTransform()
ref_srs = ref_ds.GetProjection()
ref_ds = None
xmin = ref_gt[0]
ymax = ref_gt[3]
xmax = xmin + ref_gt[1] * target_width
ymin = ymax + ref_gt[5] * target_length
dx = abs(ref_gt[1])
dy = abs(ref_gt[5])
return {
"bounds": (xmin, ymin, xmax, ymax),
"res": (dx, dy),
"srs": ref_srs,
}
```
### 2. Simplify `_warp_water_mask` using `gdal.Warp` (no subprocess)
```python
def _warp_water_mask(self, dsName, target_length, target_width):
"""Reproject water mask to match reference geometry grid via GDAL Warp."""
src_file = self.datasetDict[dsName]
grid = self._get_ref_grid(target_length, target_width)
src_ds = gdal.Open(src_file)
src_srs = src_ds.GetProjection()
src_ds = None
if not src_srs:
print(' (source water mask has no projection, assuming EPSG:4326)')
src_srs = 'EPSG:4326'
out_ds = gdal.Warp(
destNameOrDestDS='', # in-memory dataset
srcDSOrSrcDSTab=src_file,
format='MEM',
srcSRS=src_srs,
dstSRS=grid["srs"],
outputBounds=grid["bounds"],
xRes=grid["res"][0],
yRes=grid["res"][1],
resampleAlg='near',
)
if out_ds is None:
raise RuntimeError("GDAL Warp failed for water mask")
band = out_ds.GetRasterBand(1)
data = band.ReadAsArray()
out_ds = None
return data
```
This:
- Removes `tempfile`/`subprocess`/manual cleanup.
- Keeps the same CRS/bounds/resolution logic.
- Encapsulates grid computation, so `_warp_water_mask` focuses only on warping.
- Avoids the trailing `return result` after `finally`, which could reference an undefined variable on failure.
</issue_to_address>
### Comment 5
<location path="src/mintpy/objects/stackDict.py" line_range="659" />
<code_context>
subprocess.run(cmd, check=True)
</code_context>
<issue_to_address>
**security (python.lang.security.audit.dangerous-subprocess-use-audit):** Detected subprocess function 'run' without a static string. If this data can be controlled by a malicious actor, it may be an instance of command injection. Audit the use of this call to ensure it is not controllable by an external resource. You may consider using 'shlex.escape()'.
*Source: opengrep*
</issue_to_address>
### Comment 6
<location path="src/mintpy/utils/isce3_utils.py" line_range="58" />
<code_context>
tree = ET.parse(meta_file)
</code_context>
<issue_to_address>
**security (python.lang.security.use-defused-xml-parse):** The native Python `xml` library is vulnerable to XML External Entity (XXE) attacks. These attacks can leak confidential data and "XML bombs" can cause denial of service. Do not use this library to parse untrusted input. Instead the Python documentation recommends using `defusedxml`.
```suggestion
tree = defusedxml.etree.ElementTree.parse(meta_file)
```
*Source: opengrep*
</issue_to_address>Help me be more useful! Please click 👍 or 👎 on each comment and I'll use the feedback to improve your reviews.
| '-of', 'GTiff', '-co', 'COMPRESS=LZW', | ||
| src_file, tmp_f, | ||
| ] | ||
| subprocess.run(cmd, check=True) |
Contributor
There was a problem hiding this comment.
security (python.lang.security.audit.dangerous-subprocess-use-audit): Detected subprocess function 'run' without a static string. If this data can be controlled by a malicious actor, it may be an instance of command injection. Audit the use of this call to ensure it is not controllable by an external resource. You may consider using 'shlex.escape()'.
Source: opengrep
| Common metadata dictionary with keys required by MintPy. | ||
| """ | ||
| # Parse XML file | ||
| tree = ET.parse(meta_file) |
Contributor
There was a problem hiding this comment.
security (python.lang.security.use-defused-xml-parse): The native Python xml library is vulnerable to XML External Entity (XXE) attacks. These attacks can leak confidential data and "XML bombs" can cause denial of service. Do not use this library to parse untrusted input. Instead the Python documentation recommends using defusedxml.
Suggested change
| tree = ET.parse(meta_file) | |
| tree = defusedxml.etree.ElementTree.parse(meta_file) |
Source: opengrep
…F newline, pyupgrade
…used vars, empty except
tzl21
marked this pull request as draft
July 22, 2026 02:15
tzl21
marked this pull request as ready for review
July 22, 2026 08:15
Contributor
There was a problem hiding this comment.
Hey - I've found 2 security issues, 5 other issues, and left some high level feedback:
Security issues:
- Detected subprocess function 'run' without a static string. If this data can be controlled by a malicious actor, it may be an instance of command injection. Audit the use of this call to ensure it is not controllable by an external resource. You may consider using 'shlex.escape()'. (link)
- The native Python
xmllibrary is vulnerable to XML External Entity (XXE) attacks. These attacks can leak confidential data and "XML bombs" can cause denial of service. Do not use this library to parse untrusted input. Instead the Python documentation recommends usingdefusedxml. (link)
General comments:
- In
stackDict._warp_water_mask, theresultvariable is returned after thetry/finallyblock but is only defined inside thetry; ifsubprocess.runraises, this will lead to an UnboundLocalError — consider rethrowing or returning early instead of referencingresultoutside thetry. - In
read_isce3_geotiff, theDATA_TYPEfor integer rasters is hard-coded to'int16'regardless of the actual GDAL band type; useband.DataType(e.g., map GDT_UInt16/GDT_Int32, etc.) to derive a more accurate MintPyDATA_TYPEvalue.
Prompt for AI Agents
Please address the comments from this code review:
## Overall Comments
- In `stackDict._warp_water_mask`, the `result` variable is returned after the `try/finally` block but is only defined inside the `try`; if `subprocess.run` raises, this will lead to an UnboundLocalError — consider rethrowing or returning early instead of referencing `result` outside the `try`.
- In `read_isce3_geotiff`, the `DATA_TYPE` for integer rasters is hard-coded to `'int16'` regardless of the actual GDAL band type; use `band.DataType` (e.g., map GDT_UInt16/GDT_Int32, etc.) to derive a more accurate MintPy `DATA_TYPE` value.
## Individual Comments
### Comment 1
<location path="src/mintpy/objects/stackDict.py" line_range="650-659" />
<code_context>
+ import tempfile
+ fd, tmp_f = tempfile.mkstemp(suffix='.tif')
+ os.close(fd)
+ try:
+ cmd = [
+ 'gdalwarp', '-overwrite', '-q',
+ '-s_srs', src_srs, '-t_srs', ref_srs,
+ '-te', str(xmin), str(ymin), str(xmax), str(ymax),
+ '-tr', str(dx), str(dy),
+ '-r', 'near',
+ '-of', 'GTiff', '-co', 'COMPRESS=LZW',
+ src_file, tmp_f,
+ ]
+ subprocess.run(cmd, check=True)
+ tmp_ds = gdal.Open(tmp_f)
+ result = tmp_ds.GetRasterBand(1).ReadAsArray()
+ tmp_ds = None
+ src_ds = None
+ return result
+ finally:
+ if os.path.exists(tmp_f):
</code_context>
<issue_to_address>
**issue (bug_risk):** Potential NameError in _warp_water_mask if gdalwarp invocation fails before result is assigned.
Because `result` is defined only inside the `try`, any exception from `subprocess.run(cmd, check=True)` or `gdal.Open(tmp_f)` will jump to `finally` and then hit `return result`, where `result` is undefined. To avoid this, either return `result` inside the `try` and remove the return after `finally`, or initialize `result = None` before the `try` and handle the failure explicitly (e.g., raise or return the original `data`).
</issue_to_address>
### Comment 2
<location path="src/mintpy/utils/readfile.py" line_range="1437-1443" />
<code_context>
+ meta['EPSG'] = epsg
+
+ # Data type
+ band = ds.GetRasterBand(1)
+ dtype = gdal.GetDataTypeName(band.DataType).lower()
+ if 'float' in dtype:
+ meta['DATA_TYPE'] = 'float32'
+ elif 'int' in dtype:
+ meta['DATA_TYPE'] = 'int16'
+ ndv = band.GetNoDataValue()
</code_context>
<issue_to_address>
**suggestion (bug_risk):** Integer GDAL data types are all mapped to int16, which can misrepresent int32/int8 inputs.
In `read_isce3_geotiff`, `dtype = gdal.GetDataTypeName(band.DataType).lower()` and the `elif 'int' in dtype:` branch maps all integer GDAL types to `int16`. This mislabels int32, int8, and uint16 inputs, so any code using `DATA_TYPE` may allocate insufficient precision or overflow. Please branch on the specific GDAL type names (e.g., `Int16`, `UInt16`, `Int32`) or derive bit depth directly from `band.DataType` so integer types are mapped to the correct MintPy type instead of being collapsed to `int16`.
```suggestion
# Data type
band = ds.GetRasterBand(1)
gdal_dtype_name = gdal.GetDataTypeName(band.DataType)
dtype = gdal_dtype_name.lower()
# Map GDAL data types to MintPy data types without collapsing integer precision
if dtype == 'float32':
meta['DATA_TYPE'] = 'float32'
elif dtype == 'float64':
meta['DATA_TYPE'] = 'float64'
elif dtype == 'byte':
# GDAL Byte is 8‑bit unsigned; MintPy uses int8/uint8 depending on convention.
# Use uint8 here to preserve range [0, 255].
meta['DATA_TYPE'] = 'uint8'
elif dtype == 'int16':
meta['DATA_TYPE'] = 'int16'
elif dtype == 'uint16':
meta['DATA_TYPE'] = 'uint16'
elif dtype == 'int32':
meta['DATA_TYPE'] = 'int32'
elif dtype == 'uint32':
meta['DATA_TYPE'] = 'uint32'
else:
# Fallback for complex or unknown types; store the GDAL name for downstream handling.
meta['DATA_TYPE'] = gdal_dtype_name
```
</issue_to_address>
### Comment 3
<location path="src/mintpy/utils/isce3_utils.py" line_range="312-321" />
<code_context>
+ if not h5_file.exists():
+ return extracted
+
+ with h5py.File(h5_file, 'r') as h5f:
+ data_group = h5f['/data']
+
+ # EPSG code
+ epsg = 4326
+ if 'projection' in data_group:
+ proj_ds = data_group['projection']
+ if 'epsg_code' in proj_ds.attrs:
+ epsg = int(proj_ds.attrs['epsg_code'])
+ elif proj_ds.shape == () and np.issubdtype(proj_ds.dtype, np.integer):
+ epsg = int(proj_ds[()])
+
+ # Geotransform from 1D coordinate arrays (cell centers -> top-left corner)
+ x_coords = data_group['x_coordinates'][:]
+ y_coords = data_group['y_coordinates'][:]
+ dx = abs(x_coords[1] - x_coords[0])
</code_context>
<issue_to_address>
**issue:** build_vrt_from_h5 assumes coordinate arrays have at least two elements, which can fail for degenerate grids.
`dx` and `dy` are derived from `x_coords[1] - x_coords[0]` and `y_coords[1] - y_coords[0]` without validating the array lengths. For tiles with a dimension size of 1 (e.g., single row/column), this will raise an `IndexError`. Please add a length check (e.g., require `len(x_coords) > 1 and len(y_coords) > 1`) and either use a sensible default geotransform or raise a clear error for these cases.
</issue_to_address>
### Comment 4
<location path="src/mintpy/utils/isce3_utils.py" line_range="720-722" />
<code_context>
+ # az = -arctan2(east, north) * 180/pi (mod 360)
+ az_angle = -1 * np.rad2deg(np.arctan2(east, north)) % 360.0
+
+ # Mask nodata
+ if nodata is not None:
+ mask = np.isnan(east) | np.isnan(north)
+ az_angle[mask] = nodata
+
</code_context>
<issue_to_address>
**suggestion:** Nodata masking in azimuth/incidence angle computations ignores non-NaN nodata values.
`compute_azimuth_angle` / `compute_incidence_angle` only mask NaNs, so finite sentinel nodata values (e.g., 0 or a large constant) are treated as valid and affect the angle results. When `nodata` is finite, you should also mask those values (e.g., `mask |= (east == nodata) | (north == nodata)`) or use the band’s configured nodata to avoid propagating invalid LOS pixels into the angle computation.
Suggested implementation:
```python
# Azimuth from east/north components:
# az = -arctan2(east, north) * 180/pi (mod 360)
az_angle = -1 * np.rad2deg(np.arctan2(east, north)) % 360.0
```
```python
# Azimuth from east/north components:
# az = -arctan2(east, north) * 180/pi (mod 360)
az_angle = -1 * np.rad2deg(np.arctan2(east, north)) % 360.0
# Mask nodata (NaNs and finite sentinel nodata values)
if nodata is not None:
mask = np.isnan(east) | np.isnan(north)
if np.isfinite(nodata):
mask |= (east == nodata) | (north == nodata)
az_angle[mask] = nodata
import h5py
```
1. Apply the same nodata masking logic to `compute_incidence_angle` (and any other angle-computation helpers) so that both NaNs and finite nodata sentinels in LOS components are masked before computing angles.
2. Ensure that `nodata` is passed into these functions consistently, ideally using the raster band’s configured nodata value from GDAL/metadata, so the mask correctly reflects the data source.
</issue_to_address>
### Comment 5
<location path="src/mintpy/utils/isce3_utils.py" line_range="267" />
<code_context>
+ return baseline_dict
+
+
+def extract_h5_geometry(
+ h5_file: Union[str, Path],
+ output_dir: Path,
</code_context>
<issue_to_address>
**issue (complexity):** Consider factoring shared HDF5 geometry metadata logic into a helper and splitting `extract_required_attributes` into smaller helpers to reduce duplication and make the code easier to test and maintain.
The functionality here is solid, but there are a few places where you can reduce complexity without changing behavior.
---
### 1) Unify HDF5 geometry metadata handling
`extract_h5_geometry` and `build_vrt_from_h5` both:
- read `/data` group
- derive EPSG from `projection`
- compute geotransform from `x_coordinates` / `y_coordinates`
- infer nodata from `_FillValue` / dtype
You can factor that into a small helper so the logic lives in one place, while keeping the two code paths (GeoTIFF vs VRT) intact:
```python
def _read_static_layer_info(
h5_file: Union[str, Path],
dataset_mapping: Dict[str, str],
x_coords_name: str = "x_coordinates",
y_coords_name: str = "y_coordinates",
projection_name: str = "projection",
) -> Tuple[int, Tuple[float, ...], Dict[str, Dict[str, Any]]]:
"""Return (epsg, geotransform, per-dataset info) from a static_layers HDF5."""
info: Dict[str, Dict[str, Any]] = {}
with h5py.File(h5_file, "r") as h5f:
data_group = h5f["/data"]
# EPSG
epsg = 4326
if projection_name in data_group:
proj_ds = data_group[projection_name]
if "epsg_code" in proj_ds.attrs:
epsg = int(proj_ds.attrs["epsg_code"])
elif proj_ds.shape == () and np.issubdtype(proj_ds.dtype, np.integer):
epsg = int(proj_ds[()])
# Geotransform
x_coords = data_group[x_coords_name][:] if x_coords_name in data_group else None
y_coords = data_group[y_coords_name][:] if y_coords_name in data_group else None
geotransform = None
if x_coords is not None and y_coords is not None and len(x_coords) > 1 and len(y_coords) > 1:
dx = abs(x_coords[1] - x_coords[0])
dy = abs(y_coords[1] - y_coords[0])
left = x_coords[0] - dx / 2
top = y_coords[0] + dy / 2
geotransform = (left, dx, 0.0, top, 0.0, -dy)
# Per‑dataset info (shape, dtype, nodata)
for out_name, ds_name in dataset_mapping.items():
if ds_name not in data_group:
continue
ds = data_group[ds_name]
nodata = None
if "_FillValue" in ds.attrs:
nodata = ds.attrs["_FillValue"].item()
elif np.issubdtype(ds.dtype, np.floating):
nodata = np.nan
elif np.issubdtype(ds.dtype, np.integer):
nodata = np.iinfo(ds.dtype).max
info[out_name] = {
"name": ds_name,
"shape": ds.shape,
"dtype": ds.dtype,
"nodata": nodata,
}
return epsg, geotransform, info
```
Then `extract_h5_geometry` can reuse this:
```python
def extract_h5_geometry(...):
...
epsg, geotransform, ds_info = _read_static_layer_info(
h5_file, dataset_mapping, x_coords_name, y_coords_name, projection_name
)
with h5py.File(h5_file, "r") as h5f:
data_group = h5f["/data"]
for geom_type in geom_types:
info = ds_info.get(geom_type)
if not info:
continue
dataset = data_group[info["name"]]
data = dataset[:]
nodata = info["nodata"]
...
```
And `build_vrt_from_h5` can use the same helper rather than re‑deriving EPSG / geotransform / nodata:
```python
def build_vrt_from_h5(...):
...
epsg, geotransform, ds_info = _read_static_layer_info(h5_file, dataset_mapping)
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
srs_wkt = srs.ExportToWkt()
for geom_type in geom_types:
info = ds_info.get(geom_type)
if not info:
continue
length, width = info["shape"]
gdal_dtype = _GDAL_DTYPE_MAP.get(info["dtype"].name, "Float32")
nodata = info["nodata"]
...
```
This reduces duplication and keeps the VRT vs GeoTIFF optimization entirely intact.
---
### 2) Split `extract_required_attributes` into testable helpers
`extract_required_attributes` currently:
- does key selection / fallbacks
- handles unit conversion for slant range time
- performs orbit interpolation + heading + earth radius / altitude derivation
- computes frame numbers
You can keep the public API exactly as‑is but move internals into smaller helpers to make it easier to reason about and test:
```python
def _select_basic_attrs(metadata: Dict[str, Any]) -> Dict[str, Any]:
meta: Dict[str, Any] = {}
meta['prf'] = metadata.get('prf_raw_data', metadata.get('prf', 1717.0))
meta['burstStartUTC'] = metadata.get('sensing_start',
metadata.get('zero_doppler_start_time', ''))
meta['burstStopUTC'] = metadata.get('sensing_stop',
metadata.get('zero_doppler_end_time', ''))
meta['radarWavelength'] = metadata.get('wavelength',
metadata.get('radar_wavelength', 0.05546576))
meta['startingRange'] = metadata.get('starting_range',
metadata.get('slant_range_time', 800000.0))
meta['passDirection'] = metadata.get('orbit_direction',
metadata.get('orbit_pass_direction', 'ascending'))
meta['polarization'] = metadata.get('polarization', 'VV')
meta['trackNumber'] = metadata.get('track_number', 0)
meta['orbitNumber'] = metadata.get('absolute_orbit_number', 0)
meta['sensingMid'] = metadata.get('sensing_mid', metadata.get('sensing_start', ''))
meta['azimuthTimeInterval'] = metadata.get('azimuth_time_interval',
metadata.get('azimuth_time_interval_', 0.002))
meta['rangePixelSize'] = metadata.get('range_pixel_spacing',
metadata.get('range_pixel_spacing_', 2.3))
meta['swathNumber'] = metadata.get('swathNumber', metadata.get('swath_number', 2))
meta['ascendingNodeTime'] = None
return meta
```
Unit conversion and orbit‑derived fields can also be isolated:
```python
def _normalize_starting_range(meta: Dict[str, Any], metadata: Dict[str, Any]) -> None:
val = meta['startingRange']
if isinstance(val, (int, float)):
try:
if val < 1.0: # seconds -> meters
SPEED_OF_LIGHT = 299792458.0
meta['startingRange'] = val * SPEED_OF_LIGHT / 2.0
except Exception:
meta['startingRange'] = metadata.get('starting_range', 800000.0)
def _compute_orbit_fields(meta: Dict[str, Any], metadata: Dict[str, Any]) -> None:
try:
import isce3
except ImportError:
meta.update({
'satelliteSpeed': 7545.0,
'position': [0.0, 0.0, 0.0],
'HEADING': 0.0,
'earthRadius': 6371000.0,
'altitude': 693000.0,
})
return
has_orbit = all(k in metadata for k in ['time', 'position_x', 'velocity_x'])
if not has_orbit:
meta.update({...}) # same defaults as above
return
try:
sensingMid = _to_seconds(metadata['sensing_mid'], metadata['reference_epoch'])
position, velocity_vec = _orbit_interp_hermite(metadata, sensingMid)
velocity = np.linalg.norm(velocity_vec)
ellipsoid = isce3.core.Ellipsoid()
llh = ellipsoid.xyz_to_lon_lat(position)
heading = _compute_heading((position, velocity_vec))
meta['satelliteSpeed'] = velocity
meta['position'] = position
meta['HEADING'] = heading
meta['earthRadius'] = ellipsoid.r_dir(math.radians(heading), llh[1])
meta['altitude'] = llh[2]
except Exception as e:
print(f"WARNING: orbit interpolation failed: {e}. Using default values.")
meta.update({...}) # same defaults
```
And frame numbers:
```python
def _compute_frame_numbers(meta: Dict[str, Any]) -> None:
meta['firstFrameNumber'] = 0
meta['lastFrameNumber'] = 0
if meta['ascendingNodeTime'] is None or not meta['burstStartUTC']:
return
try:
start_dt = datetime.strptime(meta['burstStartUTC'], '%Y-%m-%d %H:%M:%S.%f')
stop_dt = datetime.strptime(meta['burstStopUTC'], '%Y-%m-%d %H:%M:%S.%f')
node_dt = (datetime.strptime(meta['ascendingNodeTime'], '%Y-%m-%d %H:%M:%S.%f')
if isinstance(meta['ascendingNodeTime'], str)
else meta['ascendingNodeTime'])
time_diff_start = (start_dt - node_dt).total_seconds()
time_diff_stop = (stop_dt - node_dt).total_seconds()
meta['firstFrameNumber'] = int(0.2 * time_diff_start)
meta['lastFrameNumber'] = int(0.2 * time_diff_stop)
except Exception:
pass
```
`extract_required_attributes` can then become a thin orchestrator:
```python
def extract_required_attributes(metadata: Dict[str, Any]) -> Dict[str, Any]:
meta = _select_basic_attrs(metadata)
_normalize_starting_range(meta, metadata)
_compute_orbit_fields(meta, metadata)
_compute_frame_numbers(meta)
return meta
```
This keeps the external behavior the same but makes each concern clearer and easier to test in isolation.
</issue_to_address>
### Comment 6
<location path="src/mintpy/objects/stackDict.py" line_range="660" />
<code_context>
subprocess.run(cmd, check=True)
</code_context>
<issue_to_address>
**security (python.lang.security.audit.dangerous-subprocess-use-audit):** Detected subprocess function 'run' without a static string. If this data can be controlled by a malicious actor, it may be an instance of command injection. Audit the use of this call to ensure it is not controllable by an external resource. You may consider using 'shlex.escape()'.
*Source: opengrep*
</issue_to_address>
### Comment 7
<location path="src/mintpy/utils/isce3_utils.py" line_range="59" />
<code_context>
tree = ET.parse(meta_file)
</code_context>
<issue_to_address>
**security (python.lang.security.use-defused-xml-parse):** The native Python `xml` library is vulnerable to XML External Entity (XXE) attacks. These attacks can leak confidential data and "XML bombs" can cause denial of service. Do not use this library to parse untrusted input. Instead the Python documentation recommends using `defusedxml`.
```suggestion
tree = defusedxml.etree.ElementTree.parse(meta_file)
```
*Source: opengrep*
</issue_to_address>Help me be more useful! Please click 👍 or 👎 on each comment and I'll use the feedback to improve your reviews.
Comment on lines
+650
to
+659
| try: | ||
| cmd = [ | ||
| 'gdalwarp', '-overwrite', '-q', | ||
| '-s_srs', src_srs, '-t_srs', ref_srs, | ||
| '-te', str(xmin), str(ymin), str(xmax), str(ymax), | ||
| '-tr', str(dx), str(dy), | ||
| '-r', 'near', | ||
| '-of', 'GTiff', '-co', 'COMPRESS=LZW', | ||
| src_file, tmp_f, | ||
| ] |
Contributor
There was a problem hiding this comment.
issue (bug_risk): Potential NameError in _warp_water_mask if gdalwarp invocation fails before result is assigned.
Because result is defined only inside the try, any exception from subprocess.run(cmd, check=True) or gdal.Open(tmp_f) will jump to finally and then hit return result, where result is undefined. To avoid this, either return result inside the try and remove the return after finally, or initialize result = None before the try and handle the failure explicitly (e.g., raise or return the original data).
Comment on lines
+1437
to
+1443
| # Data type | ||
| band = ds.GetRasterBand(1) | ||
| dtype = gdal.GetDataTypeName(band.DataType).lower() | ||
| if 'float' in dtype: | ||
| meta['DATA_TYPE'] = 'float32' | ||
| elif 'int' in dtype: | ||
| meta['DATA_TYPE'] = 'int16' |
Contributor
There was a problem hiding this comment.
suggestion (bug_risk): Integer GDAL data types are all mapped to int16, which can misrepresent int32/int8 inputs.
In read_isce3_geotiff, dtype = gdal.GetDataTypeName(band.DataType).lower() and the elif 'int' in dtype: branch maps all integer GDAL types to int16. This mislabels int32, int8, and uint16 inputs, so any code using DATA_TYPE may allocate insufficient precision or overflow. Please branch on the specific GDAL type names (e.g., Int16, UInt16, Int32) or derive bit depth directly from band.DataType so integer types are mapped to the correct MintPy type instead of being collapsed to int16.
Suggested change
| # Data type | |
| band = ds.GetRasterBand(1) | |
| dtype = gdal.GetDataTypeName(band.DataType).lower() | |
| if 'float' in dtype: | |
| meta['DATA_TYPE'] = 'float32' | |
| elif 'int' in dtype: | |
| meta['DATA_TYPE'] = 'int16' | |
| # Data type | |
| band = ds.GetRasterBand(1) | |
| gdal_dtype_name = gdal.GetDataTypeName(band.DataType) | |
| dtype = gdal_dtype_name.lower() | |
| # Map GDAL data types to MintPy data types without collapsing integer precision | |
| if dtype == 'float32': | |
| meta['DATA_TYPE'] = 'float32' | |
| elif dtype == 'float64': | |
| meta['DATA_TYPE'] = 'float64' | |
| elif dtype == 'byte': | |
| # GDAL Byte is 8‑bit unsigned; MintPy uses int8/uint8 depending on convention. | |
| # Use uint8 here to preserve range [0, 255]. | |
| meta['DATA_TYPE'] = 'uint8' | |
| elif dtype == 'int16': | |
| meta['DATA_TYPE'] = 'int16' | |
| elif dtype == 'uint16': | |
| meta['DATA_TYPE'] = 'uint16' | |
| elif dtype == 'int32': | |
| meta['DATA_TYPE'] = 'int32' | |
| elif dtype == 'uint32': | |
| meta['DATA_TYPE'] = 'uint32' | |
| else: | |
| # Fallback for complex or unknown types; store the GDAL name for downstream handling. | |
| meta['DATA_TYPE'] = gdal_dtype_name |
Comment on lines
+312
to
+321
| with h5py.File(h5_file, 'r') as h5f: | ||
| data_group = h5f['/data'] | ||
|
|
||
| # Get EPSG from projection dataset | ||
| epsg = 4326 | ||
| if projection_name in data_group: | ||
| proj_ds = data_group[projection_name] | ||
| if 'epsg_code' in proj_ds.attrs: | ||
| epsg = int(proj_ds.attrs['epsg_code']) | ||
| elif 'spatial_ref' in proj_ds.attrs: |
Contributor
There was a problem hiding this comment.
issue: build_vrt_from_h5 assumes coordinate arrays have at least two elements, which can fail for degenerate grids.
dx and dy are derived from x_coords[1] - x_coords[0] and y_coords[1] - y_coords[0] without validating the array lengths. For tiles with a dimension size of 1 (e.g., single row/column), this will raise an IndexError. Please add a length check (e.g., require len(x_coords) > 1 and len(y_coords) > 1) and either use a sensible default geotransform or raise a clear error for these cases.
Comment on lines
+720
to
+722
| # Mask nodata | ||
| if nodata is not None: | ||
| mask = np.isnan(east) | np.isnan(north) |
Contributor
There was a problem hiding this comment.
suggestion: Nodata masking in azimuth/incidence angle computations ignores non-NaN nodata values.
compute_azimuth_angle / compute_incidence_angle only mask NaNs, so finite sentinel nodata values (e.g., 0 or a large constant) are treated as valid and affect the angle results. When nodata is finite, you should also mask those values (e.g., mask |= (east == nodata) | (north == nodata)) or use the band’s configured nodata to avoid propagating invalid LOS pixels into the angle computation.
Suggested implementation:
# Azimuth from east/north components:
# az = -arctan2(east, north) * 180/pi (mod 360)
az_angle = -1 * np.rad2deg(np.arctan2(east, north)) % 360.0 # Azimuth from east/north components:
# az = -arctan2(east, north) * 180/pi (mod 360)
az_angle = -1 * np.rad2deg(np.arctan2(east, north)) % 360.0
# Mask nodata (NaNs and finite sentinel nodata values)
if nodata is not None:
mask = np.isnan(east) | np.isnan(north)
if np.isfinite(nodata):
mask |= (east == nodata) | (north == nodata)
az_angle[mask] = nodata
import h5py- Apply the same nodata masking logic to
compute_incidence_angle(and any other angle-computation helpers) so that both NaNs and finite nodata sentinels in LOS components are masked before computing angles. - Ensure that
nodatais passed into these functions consistently, ideally using the raster band’s configured nodata value from GDAL/metadata, so the mask correctly reflects the data source.
| return baseline_dict | ||
|
|
||
|
|
||
| def extract_h5_geometry( |
Contributor
There was a problem hiding this comment.
issue (complexity): Consider factoring shared HDF5 geometry metadata logic into a helper and splitting extract_required_attributes into smaller helpers to reduce duplication and make the code easier to test and maintain.
The functionality here is solid, but there are a few places where you can reduce complexity without changing behavior.
1) Unify HDF5 geometry metadata handling
extract_h5_geometry and build_vrt_from_h5 both:
- read
/datagroup - derive EPSG from
projection - compute geotransform from
x_coordinates/y_coordinates - infer nodata from
_FillValue/ dtype
You can factor that into a small helper so the logic lives in one place, while keeping the two code paths (GeoTIFF vs VRT) intact:
def _read_static_layer_info(
h5_file: Union[str, Path],
dataset_mapping: Dict[str, str],
x_coords_name: str = "x_coordinates",
y_coords_name: str = "y_coordinates",
projection_name: str = "projection",
) -> Tuple[int, Tuple[float, ...], Dict[str, Dict[str, Any]]]:
"""Return (epsg, geotransform, per-dataset info) from a static_layers HDF5."""
info: Dict[str, Dict[str, Any]] = {}
with h5py.File(h5_file, "r") as h5f:
data_group = h5f["/data"]
# EPSG
epsg = 4326
if projection_name in data_group:
proj_ds = data_group[projection_name]
if "epsg_code" in proj_ds.attrs:
epsg = int(proj_ds.attrs["epsg_code"])
elif proj_ds.shape == () and np.issubdtype(proj_ds.dtype, np.integer):
epsg = int(proj_ds[()])
# Geotransform
x_coords = data_group[x_coords_name][:] if x_coords_name in data_group else None
y_coords = data_group[y_coords_name][:] if y_coords_name in data_group else None
geotransform = None
if x_coords is not None and y_coords is not None and len(x_coords) > 1 and len(y_coords) > 1:
dx = abs(x_coords[1] - x_coords[0])
dy = abs(y_coords[1] - y_coords[0])
left = x_coords[0] - dx / 2
top = y_coords[0] + dy / 2
geotransform = (left, dx, 0.0, top, 0.0, -dy)
# Per‑dataset info (shape, dtype, nodata)
for out_name, ds_name in dataset_mapping.items():
if ds_name not in data_group:
continue
ds = data_group[ds_name]
nodata = None
if "_FillValue" in ds.attrs:
nodata = ds.attrs["_FillValue"].item()
elif np.issubdtype(ds.dtype, np.floating):
nodata = np.nan
elif np.issubdtype(ds.dtype, np.integer):
nodata = np.iinfo(ds.dtype).max
info[out_name] = {
"name": ds_name,
"shape": ds.shape,
"dtype": ds.dtype,
"nodata": nodata,
}
return epsg, geotransform, infoThen extract_h5_geometry can reuse this:
def extract_h5_geometry(...):
...
epsg, geotransform, ds_info = _read_static_layer_info(
h5_file, dataset_mapping, x_coords_name, y_coords_name, projection_name
)
with h5py.File(h5_file, "r") as h5f:
data_group = h5f["/data"]
for geom_type in geom_types:
info = ds_info.get(geom_type)
if not info:
continue
dataset = data_group[info["name"]]
data = dataset[:]
nodata = info["nodata"]
...And build_vrt_from_h5 can use the same helper rather than re‑deriving EPSG / geotransform / nodata:
def build_vrt_from_h5(...):
...
epsg, geotransform, ds_info = _read_static_layer_info(h5_file, dataset_mapping)
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
srs_wkt = srs.ExportToWkt()
for geom_type in geom_types:
info = ds_info.get(geom_type)
if not info:
continue
length, width = info["shape"]
gdal_dtype = _GDAL_DTYPE_MAP.get(info["dtype"].name, "Float32")
nodata = info["nodata"]
...This reduces duplication and keeps the VRT vs GeoTIFF optimization entirely intact.
2) Split extract_required_attributes into testable helpers
extract_required_attributes currently:
- does key selection / fallbacks
- handles unit conversion for slant range time
- performs orbit interpolation + heading + earth radius / altitude derivation
- computes frame numbers
You can keep the public API exactly as‑is but move internals into smaller helpers to make it easier to reason about and test:
def _select_basic_attrs(metadata: Dict[str, Any]) -> Dict[str, Any]:
meta: Dict[str, Any] = {}
meta['prf'] = metadata.get('prf_raw_data', metadata.get('prf', 1717.0))
meta['burstStartUTC'] = metadata.get('sensing_start',
metadata.get('zero_doppler_start_time', ''))
meta['burstStopUTC'] = metadata.get('sensing_stop',
metadata.get('zero_doppler_end_time', ''))
meta['radarWavelength'] = metadata.get('wavelength',
metadata.get('radar_wavelength', 0.05546576))
meta['startingRange'] = metadata.get('starting_range',
metadata.get('slant_range_time', 800000.0))
meta['passDirection'] = metadata.get('orbit_direction',
metadata.get('orbit_pass_direction', 'ascending'))
meta['polarization'] = metadata.get('polarization', 'VV')
meta['trackNumber'] = metadata.get('track_number', 0)
meta['orbitNumber'] = metadata.get('absolute_orbit_number', 0)
meta['sensingMid'] = metadata.get('sensing_mid', metadata.get('sensing_start', ''))
meta['azimuthTimeInterval'] = metadata.get('azimuth_time_interval',
metadata.get('azimuth_time_interval_', 0.002))
meta['rangePixelSize'] = metadata.get('range_pixel_spacing',
metadata.get('range_pixel_spacing_', 2.3))
meta['swathNumber'] = metadata.get('swathNumber', metadata.get('swath_number', 2))
meta['ascendingNodeTime'] = None
return metaUnit conversion and orbit‑derived fields can also be isolated:
def _normalize_starting_range(meta: Dict[str, Any], metadata: Dict[str, Any]) -> None:
val = meta['startingRange']
if isinstance(val, (int, float)):
try:
if val < 1.0: # seconds -> meters
SPEED_OF_LIGHT = 299792458.0
meta['startingRange'] = val * SPEED_OF_LIGHT / 2.0
except Exception:
meta['startingRange'] = metadata.get('starting_range', 800000.0)
def _compute_orbit_fields(meta: Dict[str, Any], metadata: Dict[str, Any]) -> None:
try:
import isce3
except ImportError:
meta.update({
'satelliteSpeed': 7545.0,
'position': [0.0, 0.0, 0.0],
'HEADING': 0.0,
'earthRadius': 6371000.0,
'altitude': 693000.0,
})
return
has_orbit = all(k in metadata for k in ['time', 'position_x', 'velocity_x'])
if not has_orbit:
meta.update({...}) # same defaults as above
return
try:
sensingMid = _to_seconds(metadata['sensing_mid'], metadata['reference_epoch'])
position, velocity_vec = _orbit_interp_hermite(metadata, sensingMid)
velocity = np.linalg.norm(velocity_vec)
ellipsoid = isce3.core.Ellipsoid()
llh = ellipsoid.xyz_to_lon_lat(position)
heading = _compute_heading((position, velocity_vec))
meta['satelliteSpeed'] = velocity
meta['position'] = position
meta['HEADING'] = heading
meta['earthRadius'] = ellipsoid.r_dir(math.radians(heading), llh[1])
meta['altitude'] = llh[2]
except Exception as e:
print(f"WARNING: orbit interpolation failed: {e}. Using default values.")
meta.update({...}) # same defaultsAnd frame numbers:
def _compute_frame_numbers(meta: Dict[str, Any]) -> None:
meta['firstFrameNumber'] = 0
meta['lastFrameNumber'] = 0
if meta['ascendingNodeTime'] is None or not meta['burstStartUTC']:
return
try:
start_dt = datetime.strptime(meta['burstStartUTC'], '%Y-%m-%d %H:%M:%S.%f')
stop_dt = datetime.strptime(meta['burstStopUTC'], '%Y-%m-%d %H:%M:%S.%f')
node_dt = (datetime.strptime(meta['ascendingNodeTime'], '%Y-%m-%d %H:%M:%S.%f')
if isinstance(meta['ascendingNodeTime'], str)
else meta['ascendingNodeTime'])
time_diff_start = (start_dt - node_dt).total_seconds()
time_diff_stop = (stop_dt - node_dt).total_seconds()
meta['firstFrameNumber'] = int(0.2 * time_diff_start)
meta['lastFrameNumber'] = int(0.2 * time_diff_stop)
except Exception:
passextract_required_attributes can then become a thin orchestrator:
def extract_required_attributes(metadata: Dict[str, Any]) -> Dict[str, Any]:
meta = _select_basic_attrs(metadata)
_normalize_starting_range(meta, metadata)
_compute_orbit_fields(meta, metadata)
_compute_frame_numbers(meta)
return metaThis keeps the external behavior the same but makes each concern clearer and easier to test in isolation.
| '-of', 'GTiff', '-co', 'COMPRESS=LZW', | ||
| src_file, tmp_f, | ||
| ] | ||
| subprocess.run(cmd, check=True) |
Contributor
There was a problem hiding this comment.
security (python.lang.security.audit.dangerous-subprocess-use-audit): Detected subprocess function 'run' without a static string. If this data can be controlled by a malicious actor, it may be an instance of command injection. Audit the use of this call to ensure it is not controllable by an external resource. You may consider using 'shlex.escape()'.
Source: opengrep
|
|
||
| """ | ||
| # Parse XML file (Python 3.8+ disables entity expansion by default) | ||
| tree = ET.parse(meta_file) |
Contributor
There was a problem hiding this comment.
security (python.lang.security.use-defused-xml-parse): The native Python xml library is vulnerable to XML External Entity (XXE) attacks. These attacks can leak confidential data and "XML bombs" can cause denial of service. Do not use this library to parse untrusted input. Instead the Python documentation recommends using defusedxml.
Suggested change
| tree = ET.parse(meta_file) | |
| tree = defusedxml.etree.ElementTree.parse(meta_file) |
Source: opengrep
tzl21
marked this pull request as draft
July 22, 2026 08:34
…cy D212/D213 conflicts
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Add ISCE3/Dolphin processor support to MintPy. ISCE3 is the successor
to ISCE2 and produces TOPS burst-level interferograms in HDF5/GeoTIFF
format. This PR adds
prep_isce3.pyand supporting utilities, enablingMintPy to directly ingest ISCE3 outputs.
New files (3)
src/mintpy/cli/prep_isce3.pysrc/mintpy/prep_isce3.pysrc/mintpy/utils/isce3_utils.pyModified files (8)
load_data.py'isce3'toPROCESSOR_LIST; dispatch viaprep_isce3.pyobjects/stackDict.py_warp_water_mask()auto-aligns water masks viagdalwarp; robust baseline metadataobjects/stack.pyutils/readfile.pyread_isce3_geotiff()for Dolphin.int.tiffilesifgram_inversion.pyALOOKS/RLOOKSpyproject.tomlprep_isce3.pyas console_scripts entrydefaults/smallbaselineApp.cfgmintpy.load.geomSrcDirtemplate entrytropo_pyaps3.pyUsage
Typical ISCE3/Dolphin output layout
project/
├── CSLC/ # burst-level SLC products
│ └── <burst_subdir>/
│ └── static_layers.h5 # geometry, coordinates, metadata per burst
├── stitched/
│ ├── unwrapped/
│ │ └── *.unw.tif # unwrapped phase (GeoTIFF)
│ │ └── .unw.conncomp.tif # connected components
│ └── ifgrams_filtered/
│ └── .coh.tif # spatial coherence (GeoTIFF)
├── merged/
│ ├── baselines/
│ │ └── YYYYMMDD_YYYYMMDD.txt # ISCE3 baseline text files
│ └── geometry/ # prep_isce3 output (auto-generated)
│ ├── height.tif
│ ├── incidenceAngle.tif
│ └── ...
└── DEM/
└── watermask.tif
Configuration (
smallbaselineApp.cfg)Minimal CLI
prep_isce3.py -m reference/IW1.xml -g ./CSLC// -b ./merged/baselines
-f ./stitched/unwrapped/.unw.tif
Key differences from ISCE2
Item ISCE2 (processor = isce)
Reference metadata .rsc sidecar or reference .unw
Interferograms Binary + .rsc
Geometry Pre-merged, one file per layer
Baselines ISCE2 baseline text files
Water mask Matches interferogram grid
New option: mintpy.load.geomSrcDir
Path to burst-level static_layers*.h5 files. Supports glob patterns
(e.g. ./CSLC/*/) for multi-burst datasets. Internally:
across burst subdirectories.
(zero-copy, no raster duplication).
single grid, clipped to the interferogram extent from the reference
GeoTIFF.
components when available.
New option: mintpy.load.metaFile = auto
When set to auto (or omitted), metadata is extracted automatically from
the first available static_layers*.h5 file. No manual reference XML
needed. If a specific burst XML is desired (e.g. reference/IW1.xml),
pass it explicitly.
Design notes
GMTSAR, NISAR, or other processors.
Bperp average (m) column headers from ISCE3.
Checklist
Note on Codacy / pydocstyle D212/D213
Multiple Codacy issues report docstring formatting problems (D212, D213).
These two rules are mutually exclusive in pydocstyle — D212 requires
the summary to start on the first line of the docstring (
"""Summary...),while D213 requires it on the second line (
"""\n Summary...). PEP 257does not mandate either format, and both are widely used in the Python
ecosystem.
This PR adopts D212 (summary on the first line, matching the prevailing
style in the MintPy codebase). All docstrings in the new ISCE3 files are
uniformly formatted to D212. The remaining D212/D213 Codacy alerts are a
tool configuration limitation (project-level
.pydocstylenot configuredto ignore the unused rule) and do not reflect code quality issues
Summary by Sourcery
Add ISCE3/Dolphin processor support for loading geocoded interferograms, including geometry preparation, metadata extraction, and CLI integration.
New Features:
Bug Fixes:
Enhancements:
Build:
Documentation: