diff --git a/.github/workflows/build-n-publish-to-pypi.yml b/.github/workflows/build-n-publish-to-pypi.yml index 5300f7de0..bc823cfe8 100644 --- a/.github/workflows/build-n-publish-to-pypi.yml +++ b/.github/workflows/build-n-publish-to-pypi.yml @@ -21,7 +21,7 @@ jobs: fetch-depth: 0 - name: Set up Python 3.10 - uses: actions/setup-python@v6 + uses: actions/setup-python@v7 with: python-version: "3.10" diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 88deec22a..20a4d97e7 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -27,7 +27,7 @@ repos: exclude: tests/data/ - repo: https://github.com/PyCQA/isort - rev: "9.0.0b1" + rev: "9.0.0b5" hooks: - id: isort name: sort imports diff --git a/pyproject.toml b/pyproject.toml index 22880a5cf..f5472562c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -67,6 +67,7 @@ Issues = "https://github.com/insarlab/MintPy/issues" "prep_gmtsar.py" = "mintpy.cli.prep_gmtsar:main" "prep_hyp3.py" = "mintpy.cli.prep_hyp3:main" "prep_isce.py" = "mintpy.cli.prep_isce:main" +"prep_isce3.py" = "mintpy.cli.prep_isce3:main" "prep_nisar.py" = "mintpy.cli.prep_nisar:main" "prep_roipac.py" = "mintpy.cli.prep_roipac:main" "prep_snap.py" = "mintpy.cli.prep_snap:main" diff --git a/src/mintpy/cli/prep_isce3.py b/src/mintpy/cli/prep_isce3.py new file mode 100644 index 000000000..9bcd8446c --- /dev/null +++ b/src/mintpy/cli/prep_isce3.py @@ -0,0 +1,117 @@ +#!/usr/bin/env python3 +############################################################ +# Program is part of MintPy # +# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi # +# Author: Zhang Yunjun, 2024 # +############################################################ + + +import glob +import os +import sys + +from mintpy.utils.arg_utils import create_argument_parser +from mintpy.utils.isce3_utils import GEOMETRY_FILENAMES + +######################################################################### +# Default geometry files to extract from static_layers HDF5 +# (single source of truth defined in mintpy.utils.isce3_utils) + +EXAMPLE = """example: + ## Dolphin/ISCE-3 topsStack (auto‑generate metadata) + prep_isce3.py -f "../../dolphin/unwrapped/*.unw.tif" -b ../baselines -g ../merged/geom/ + + ## with existing metadata file + prep_isce3.py -f "../../dolphin/unwrapped/*.unw.tif" -m ../reference/IW1.xml -g ../merged/geom/ + + ## force overwrite existing .rsc files + prep_isce3.py -f "../../dolphin/unwrapped/*.unw.tif" -g ../merged/geom/ --force +""" + +def create_parser(subparsers=None): + """Command line parser.""" + synopsis = 'Prepare ISCE-3 / Dolphin metadata files.' + epilog = EXAMPLE + name = __name__.split('.')[-1] + parser = create_argument_parser( + name, synopsis=synopsis, description=synopsis, epilog=epilog, subparsers=subparsers) + + # observations + parser.add_argument('-f', dest='obs_files', type=str, nargs='+', + default=['../../dolphin/unwrapped/*.unw.tif'], + help='Wildcard path pattern(s) for observation files.\n' + 'E.g.: unwrapped phase: ../../dolphin/unwrapped/*.unw.tif\n' + ' coherence: ../../dolphin/interferograms/*.int.cor.tif\n' + ' wrapped phase: ../../dolphin/interferograms/*.int.tif\n' + ' connected comp: ../../dolphin/unwrapped/*.unw.conncomp.tif') + + # metadata + parser.add_argument('-m', '--meta-file', dest='meta_file', type=str, default=None, + help='Metadata file to extract common metadata for the stack.\n' + 'E.g.: reference burst XML (e.g., reference/IW1.xml). ' + 'If not provided, one will be generated from static_layers.h5.') + + # geometry and baseline + parser.add_argument('-b', '--baseline-dir', dest='baseline_dir', type=str, default=None, + help='Directory with baseline files (e.g., ../baselines). ' + 'If omitted, baseline info will not be added to metadata.') + parser.add_argument('-g', '--geometry-dir', dest='geom_dir', type=str, required=True, + help='Directory containing burst subdirectories with static_layers*.h5 files.\n' + 'E.g.: ../merged/geom/') + parser.add_argument('--geom-files', dest='geom_files', type=str, nargs='*', + default=GEOMETRY_FILENAMES, + help='List of geometry file basenames to extract/merge. Default: %(default)s.') + + # processing flag + parser.add_argument('--force', dest='update_mode', action='store_false', + help='Force to overwrite all .rsc metadata files (disable update mode).') + + parser.add_argument('--out-dir', dest='out_dir', type=str, default=None, + help='Output directory for merged geometry files. ' + 'If not provided, defaults to (geometry_dir)/../merged_geom') + + return parser + + +def cmd_line_parse(iargs=None): + parser = create_parser() + inps = parser.parse_args(args=iargs) + + if inps.meta_file and '*' in inps.meta_file: + fnames = glob.glob(inps.meta_file) + if fnames: + inps.meta_file = fnames[0] + else: + raise FileNotFoundError(inps.meta_file) + + # Expand glob patterns in geometry directory (e.g. "../../t124*/20210104/") + inps.geom_dirs = [inps.geom_dir] + if inps.geom_dir and ('*' in inps.geom_dir or '?' in inps.geom_dir): + matches = sorted(glob.glob(inps.geom_dir)) + if matches: + inps.geom_dir = matches[0] + inps.geom_dirs = matches + + # Set default output directory if not provided + if inps.out_dir is None: + inps.out_dir = os.path.join(os.path.dirname(inps.geom_dir), 'merged_geom') + + inps.processor = 'tops' + return inps + + +######################################################################### +def main(iargs=None): + # parse + inps = cmd_line_parse(iargs) + + # import core function + from mintpy.prep_isce3 import prep_isce3 + + # run + prep_isce3(inps) + + +######################################################################### +if __name__ == '__main__': + main(sys.argv[1:]) diff --git a/src/mintpy/defaults/smallbaselineApp.cfg b/src/mintpy/defaults/smallbaselineApp.cfg index c8dca0196..b2cd2abe5 100644 --- a/src/mintpy/defaults/smallbaselineApp.cfg +++ b/src/mintpy/defaults/smallbaselineApp.cfg @@ -26,14 +26,14 @@ mintpy.compute.config = auto #[none / slurm / pbs / lsf ], auto for none (sam ## no - save 0% disk usage, fast [default] ## lzf - save ~57% disk usage, relative slow ## gzip - save ~62% disk usage, very slow [not recommend] -mintpy.load.processor = auto #[isce, aria, hyp3, gmtsar, snap, gamma, roipac, nisar], auto for isce +mintpy.load.processor = auto #[isce, isce3, aria, hyp3, gmtsar, snap, gamma, roipac, nisar], auto for isce mintpy.load.autoPath = auto #[yes / no], auto for no, use pre-defined auto path mintpy.load.updateMode = auto #[yes / no], auto for yes, skip re-loading if HDF5 files are complete mintpy.load.compression = auto #[gzip / lzf / none / default], auto for default (none/lzf for stack/geometry). mintpy.load.frequency = auto #[auto / A / B], auto for A, NISAR only ##---------for ISCE only: -mintpy.load.metaFile = auto #[path of common metadata file for the stack], i.e.: ./reference/IW1.xml, ./referenceShelve/data.dat -mintpy.load.baselineDir = auto #[path of the baseline dir], i.e.: ./baselines +mintpy.load.metaFile = auto #[path of common metadata file for the stack], i.e.: ./reference/IW1.xml, ./referenceShelve/data.dat. For isce3: reference burst XML; auto generates one from static_layers.h5 +mintpy.load.baselineDir = auto #[path of the baseline dir], i.e.: ./baselines. For isce3/Dolphin: the burst dir, e.g. ../../baselines/t124*/ ##---------interferogram stack: mintpy.load.unwFile = auto #[path pattern of unwrapped interferogram files] mintpy.load.corFile = auto #[path pattern of spatial coherence files] @@ -51,6 +51,7 @@ mintpy.load.azOffStdFile = auto #[path pattern of azimuth offset variance fi mintpy.load.rgOffStdFile = auto #[path pattern of range offset variance file], optional but recommended mintpy.load.offSnrFile = auto #[path pattern of offset signal-to-noise ratio file], optional ##---------geometry: +mintpy.load.geomSrcDir = auto #[path of static-layer h5 files for geometry info], optional, default to the directory of DEM file. For isce3: burst dir(s) with static_layers*.h5, e.g. ../../t124*/20210104/ mintpy.load.demFile = auto #[path of DEM file] mintpy.load.lookupYFile = auto #[path of latitude /row /y coordinate file], not required for geocoded data mintpy.load.lookupXFile = auto #[path of longitude/column/x coordinate file], not required for geocoded data diff --git a/src/mintpy/ifgram_inversion.py b/src/mintpy/ifgram_inversion.py index 16c8e1dc2..daa7ede04 100644 --- a/src/mintpy/ifgram_inversion.py +++ b/src/mintpy/ifgram_inversion.py @@ -521,7 +521,7 @@ def calc_weight_sqrt(stack_obj, box, weight_func='var', dropIfgram=True, chunk_s L = float(stack_obj.metadata['NCORRLOOKS']) else: # use the typical ratio of resolution vs pixel size of Sentinel-1 IW mode - L = int(stack_obj.metadata['ALOOKS']) * int(stack_obj.metadata['RLOOKS']) + L = int(float(stack_obj.metadata['ALOOKS'])) * int(float(stack_obj.metadata['RLOOKS'])) L /= 1.94 # make sure L >= 1 L = max(np.rint(L).astype(int), 1) diff --git a/src/mintpy/load_data.py b/src/mintpy/load_data.py index 1ec38a044..277186e7a 100644 --- a/src/mintpy/load_data.py +++ b/src/mintpy/load_data.py @@ -24,7 +24,7 @@ from mintpy.utils import ptime, readfile, utils as ut ################################################################# -PROCESSOR_LIST = ['isce', 'aria', 'hyp3', 'gmtsar', 'snap', 'gamma', 'roipac', 'cosicorr', 'nisar'] +PROCESSOR_LIST = ['isce', 'aria', 'hyp3', 'gmtsar', 'snap', 'gamma', 'roipac', 'cosicorr', 'nisar', 'isce3'] # primary observation dataset names OBS_DSET_NAMES = ['unwrapPhase', 'rangeOffset', 'azimuthOffset'] @@ -449,7 +449,7 @@ def read_inps_dict2geometry_dict_object(iDict, dset_name2template_key): # for processors with lookup table in geo-coordinates, remove latitude/longitude dset_name2template_key.pop('latitude') dset_name2template_key.pop('longitude') - elif iDict['processor'] in ['aria', 'gmtsar', 'hyp3', 'snap', 'cosicorr']: + elif iDict['processor'] in ['aria', 'gmtsar', 'hyp3', 'snap', 'cosicorr', 'isce3']: # for processors with geocoded products support only, do nothing for now. # check again when adding products support in radar-coordiantes pass @@ -676,6 +676,49 @@ def prepare_metadata(iDict): except: warnings.warn('prep_nisar.py failed. Assuming its result exists and continue...') + elif processor == 'isce3': + + meta_files = sorted(glob.glob(iDict['mintpy.load.metaFile'])) if iDict.get('mintpy.load.metaFile') else [] + meta_file = meta_files[0] if meta_files else 'auto' + + baseline_dir = iDict.get('mintpy.load.baselineDir', None) + + # Geometry source directory (contains burst subdirs with HDF5) + geom_src_dir = iDict.get('mintpy.load.geomSrcDir', None) + if geom_src_dir is None or geom_src_dir.lower() == 'auto': + dem_path = iDict.get('mintpy.load.demFile', '') + if dem_path and dem_path.lower() != 'auto': + geom_src_dir = os.path.dirname(dem_path) + else: + geom_src_dir = os.path.abspath('.') + # NOTE: keep glob patterns (e.g. "../CSLC/t124*/20240611/") as-is, + # prep_isce3.py expands them to support multiple burst directories. + geom_src_dir = os.path.abspath(geom_src_dir) + + # Output directory for merged geometry (same as demFile's directory) + dem_file = iDict.get('mintpy.load.demFile', '') + if dem_file and dem_file.lower() != 'auto': + out_dir = os.path.dirname(os.path.abspath(dem_file)) + else: + first_match = (sorted(glob.glob(geom_src_dir)) or [geom_src_dir])[0] + out_dir = os.path.join(os.path.dirname(first_match), 'merged_geom') + + obs_keys = ['mintpy.load.unwFile', 'mintpy.load.corFile', 'mintpy.load.connCompFile'] + obs_paths = [iDict[key] for key in obs_keys if iDict.get(key, 'auto').lower() != 'auto'] + obs_paths = [x for x in obs_paths if glob.glob(x)] + + iargs = ['-m', meta_file, '-g', geom_src_dir, '--out-dir', out_dir] + if baseline_dir: + iargs += ['-b', baseline_dir] + if obs_paths: + iargs += ['-f'] + obs_paths + if not iDict.get('updateMode', True): + iargs.append('--force') + + ut.print_command_line('prep_isce3.py', iargs) + prep_module = importlib.import_module('mintpy.cli.prep_isce3') + prep_module.main(iargs) + elif processor == 'isce': from mintpy.utils import isce_utils, s1_utils diff --git a/src/mintpy/objects/sensor.py b/src/mintpy/objects/sensor.py index 42e62adcc..3e54ec3de 100644 --- a/src/mintpy/objects/sensor.py +++ b/src/mintpy/objects/sensor.py @@ -476,6 +476,7 @@ def get_unavco_mission_name(meta_dict): 'IW2' : {'range_resolution' : 3.1, 'azimuth_resolution': 22.7}, 'IW3' : {'range_resolution' : 3.5, 'azimuth_resolution': 22.6}, 'noise_equivalent_sigma_zero': -22, # dB + 'incidence_angle' : [20, 47], # degrees for Strip Map mode; 31-46 for IW mode } @@ -651,6 +652,7 @@ def get_unavco_mission_name(meta_dict): 'chirp_bandwidth' : 84e6, # Hz, 84/42/28 'range_resolution' : 3, # m 'noise_equivalent_sigma_zero': -20, # dB, -20/-24/-28 + 'incidence_angle' : [30, 56], # degrees for Strip Map mode } # SAOCOM-1A/B stripmap @@ -695,6 +697,7 @@ def get_unavco_mission_name(meta_dict): 'range_pixel_size' : 1.67, # m 'range_resolution' : 2.50, # m 'noise_equivalent_sigma_zero': -28, # dB + 'incidence_angle' : [20,46], # degree for STRIP1/2 InSAR } # UAVSAR-L @@ -744,6 +747,7 @@ def get_unavco_mission_name(meta_dict): '80MHz' : 1.87, # m }, 'noise_equivalent_sigma_zero': -25, # dB + 'incidence_angle' : [34, 48], # degrees } diff --git a/src/mintpy/objects/stack.py b/src/mintpy/objects/stack.py index c35068590..b67be407d 100644 --- a/src/mintpy/objects/stack.py +++ b/src/mintpy/objects/stack.py @@ -12,6 +12,7 @@ import itertools import os import time +import warnings import h5py import numpy as np @@ -993,6 +994,10 @@ def nonzero_mask(self, datasetName=None, print_msg=True, dropIfgram=True): for i in range(num2read): prog_bar.update(i+1, suffix=f'{i+1}/{num2read}') data = dset[idx2read[i], :, :] + if np.all(data == 0.) or np.all(np.isnan(data)): + if print_msg: + print(f'WARNING: ifgram {idx2read[i]} has all-zero/all-NaN {datasetName}, skipping') + continue mask[data == 0.] = 0 mask[np.isnan(data)] = 0 prog_bar.close() @@ -1048,6 +1053,11 @@ def temporal_average(self, datasetName='coherence', dropIfgram=True, max_memory= # referencing / normalizing for phase if 'unwrapPhase' in datasetName: + # keep MintPy's convention that a phase value of 0.0 indicates + # no-data: convert them to NaN so they are excluded from the + # temporal average (np.nanmean below), consistent with the + # data != 0. check in read_stack_obs(). + data[data == 0.] = np.nan # spatial referencing if ref_val is not None: data -= np.tile(ref_val.reshape(-1, 1, 1), (1, data.shape[1], data.shape[2])) @@ -1056,7 +1066,11 @@ def temporal_average(self, datasetName='coherence', dropIfgram=True, max_memory= data[j,:,:] *= (phase2range / tbase[j]) # use nanmean to better handle NaN values - dmean[r0:r1, :] = np.nanmean(data, axis=0) + # suppress the Mean of empty slice RuntimeWarning for pixels + # with no valid observation (all 0.0 / NaN) + with warnings.catch_warnings(): + warnings.simplefilter('ignore', RuntimeWarning) + dmean[r0:r1, :] = np.nanmean(data, axis=0) prog_bar.close() return dmean diff --git a/src/mintpy/objects/stackDict.py b/src/mintpy/objects/stackDict.py index 54c340f1d..95a324e0e 100644 --- a/src/mintpy/objects/stackDict.py +++ b/src/mintpy/objects/stackDict.py @@ -18,6 +18,7 @@ import h5py import numpy as np +from osgeo import gdal from skimage.transform import resize from mintpy.multilook import multilook_data @@ -391,8 +392,8 @@ def get_size(self, family=IFGRAM_DSET_NAMES[0]): def get_perp_baseline(self, family=IFGRAM_DSET_NAMES[0]): self.file = self.datasetDict[family] metadata = readfile.read_attribute(self.file) - self.bperp_top = float(metadata['P_BASELINE_TOP_HDR']) - self.bperp_bottom = float(metadata['P_BASELINE_BOTTOM_HDR']) + self.bperp_top = float(metadata.get('P_BASELINE_TOP_HDR', 0)) + self.bperp_bottom = float(metadata.get('P_BASELINE_BOTTOM_HDR', 0)) self.bperp = (self.bperp_top + self.bperp_bottom) / 2.0 return self.bperp @@ -608,6 +609,71 @@ def get_metadata(self, family=None): return self.metadata + 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) + """ + import tempfile + + src_file = self.datasetDict[dsName] + 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]) + + src_ds = gdal.Open(src_file) + src_srs = src_ds.GetProjection() + src_ds = None + if not src_srs: + # source file may lack embedded projection (e.g. GeoTIFF + # written before PROJ_DATA was set). Assume EPSG:4326. + print(' (source water mask has no projection, assuming EPSG:4326)') + src_srs = 'EPSG:4326' + + # Use the GDAL Python API (same engine as the gdalwarp CLI), with a + # guaranteed temp-file cleanup via try/finally. + fd, tmp_f = tempfile.mkstemp(suffix='.tif') + os.close(fd) + try: + warp_kwargs = dict( + format='GTiff', + srcSRS=src_srs, + dstSRS=ref_srs, + outputBounds=(xmin, ymin, xmax, ymax), + xRes=dx, + yRes=dy, + resampleAlg='near', + creationOptions=['COMPRESS=LZW'], + ) + tmp_ds = gdal.Warp(tmp_f, src_file, **warp_kwargs) + if tmp_ds is None: + raise RuntimeError(f'gdal.Warp failed for water mask: {src_file}') + result = tmp_ds.GetRasterBand(1).ReadAsArray() + tmp_ds = None + finally: + if os.path.exists(tmp_f): + os.remove(tmp_f) + + # verify the output grid matches the reference geometry grid + # (float rounding in -te/-tr can rarely produce a +/-1 pixel difference) + if result.shape != (target_length, target_width): + print(f' WARNING: warped waterMask size {result.shape} != ' + f'target ({target_length}, {target_width}); cropping to target size') + result = result[:target_length, :target_width] + return result + def write2hdf5(self, outputFile='geometryRadar.h5', access_mode='w', box=None, xstep=1, ystep=1, compression='lzf', extra_metadata=None): """Save/write to HDF5 file with structure defined in: @@ -706,6 +772,12 @@ def write2hdf5(self, outputFile='geometryRadar.h5', access_mode='w', box=None, x 'convert to water mask (False/True for water/land).'.format(fname)) elif dsName == 'waterMask': + # auto-align if source has different resolution/CRS + # (e.g. 1-arcsec global water mask vs 20 m interferogram grid) + if data.shape != (length, width): + print(f' auto-aligning waterMask: {data.shape} -> ({length}, {width})') + data = self._warp_water_mask(dsName, length, width) + # GMTSAR water/land mask: 1 for land, and nan for water / no data if np.sum(np.isnan(data)) > 0: print(' convert NaN value for waterMask to zero.') diff --git a/src/mintpy/prep_isce3.py b/src/mintpy/prep_isce3.py new file mode 100644 index 000000000..17bc75357 --- /dev/null +++ b/src/mintpy/prep_isce3.py @@ -0,0 +1,325 @@ +import glob +import os +import re +from pathlib import Path + +from mintpy.utils import isce3_utils, ptime, readfile, writefile + + +######################################################################### +def add_ifgram_metadata(metadata_in, dates=None, baseline_dict=None): + """Add metadata unique for each interferogram. + + Parameters: metadata_in : dict, input common metadata for the entire dataset + dates : list of str in YYYYMMDD format + baseline_dict : dict, output of baseline_timeseries() + Returns: metadata : dict, updated metadata + """ + dates = dates or [] + baseline_dict = baseline_dict or {} + metadata = metadata_in.copy() + metadata['DATE12'] = f'{dates[0][2:]}-{dates[1][2:]}' + + if baseline_dict: + if dates[0] in baseline_dict and dates[1] in baseline_dict: + bperp_top = baseline_dict[dates[1]][0] - baseline_dict[dates[0]][0] + bperp_bottom = baseline_dict[dates[1]][1] - baseline_dict[dates[0]][1] + metadata['P_BASELINE_TOP_HDR'] = str(bperp_top) + metadata['P_BASELINE_BOTTOM_HDR'] = str(bperp_bottom) + 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): + """Prepare geometry files from ISCE3/Dolphin static_layers HDF5. + + Parameters + ---------- + geom_dir : str + Base geometry directory containing burst subdirectories. + geom_dirs : list of str, optional + Additional directories with static_layers*.h5 files. + geom_files : list of str, optional + Direct paths to static_layers*.h5 files. Defaults to + isce3_utils.GEOMETRY_FILENAMES. + metadata : dict + Metadata dictionary to be updated with LENGTH/WIDTH. + out_dir : str + Output directory for merged geometry files. + processor : str + Processor name (e.g., 'isce3'). + ref_int_file : str, optional + Reference interferogram for output extent. + target_shape : tuple of (length, width), optional + Interferogram dimensions. + update_mode : bool, optional + Update mode (not used here, for consistency). + + """ + print('preparing geometry files from ISCE3/Dolphin static layers') + geom_dir = os.path.abspath(geom_dir) + out_dir = os.path.abspath(out_dir) + os.makedirs(out_dir, exist_ok=True) + + if geom_files is None: + geom_files = isce3_utils.GEOMETRY_FILENAMES + + # Step 0: Read full-resolution pixel size from first static_layers.h5 + burst_full_dx = None + burst_full_dy = None + num_bursts = 0 + 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] + elif h5_files_root: + first_h5_path = h5_files_root[0] + num_bursts = len(burst_subdirs) if burst_subdirs else (1 if h5_files_root else 0) + if geom_dirs and len(geom_dirs) > 1: + for edir in geom_dirs[1:]: + epath = Path(edir) + if epath.is_dir(): + 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: + try: + import h5py + 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]) + 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}') + + # Step 1: Merge and crop geometry to interferogram extent and resolution + extra_dirs = geom_dirs[1:] if (geom_dirs and len(geom_dirs) > 1) else None + geometry_dict = isce3_utils.extract_merge_geometry( + geom_dir=geom_dir, + output_dir=out_dir, + geom_types=geom_files, + ref_int_file=ref_int_file, + metadata=None, + extra_dirs=extra_dirs, + ) + + # Step 2: Compute ALOOKS/RLOOKS from full-res pixel size vs target pixel size + if metadata is not None: + lks_y = 1 + 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})') + + # Write .rsc files + 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' + # Remove any stale .rsc from a previous run when the merged GeoTIFF + # was regenerated, so that read_attribute() reads the fresh tif + # (readfile gives the .rsc sidecar priority over the .tif itself). + if (update_mode and os.path.isfile(rsc_file) + and os.path.getmtime(geom_path) > os.path.getmtime(rsc_file)): + os.remove(rsc_file) + 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) + + # Update metadata with final dimensions + height_file = geometry_dict.get('height.tif') + if height_file and os.path.isfile(height_file): + geom_atr = readfile.read_attribute(str(height_file)) + if metadata is not None: + metadata['LENGTH'] = geom_atr['LENGTH'] + metadata['WIDTH'] = geom_atr['WIDTH'] + print(f"Updated metadata with LENGTH={metadata['LENGTH']}, WIDTH={metadata['WIDTH']}") + + return metadata + + +def prepare_stack_isce3(obs_file, metadata=None, baseline_dict=None, update_mode=True): + print(f'preparing RSC file for: {obs_file}') + tif_files = sorted(glob.glob(obs_file)) + if not tif_files: + raise FileNotFoundError(f'No file found with pattern: {obs_file}') + + meta = metadata.copy() if metadata else {} + num_file = len(tif_files) + + # Normalize ALOOKS/RLOOKS to integer strings if present + for key in ['ALOOKS', 'RLOOKS']: + if key in meta: + try: + meta[key] = str(int(float(meta[key]))) + except (ValueError, TypeError): + pass + + prog_bar = ptime.progressBar(maxValue=num_file, print_msg=(num_file > 5)) + num_valid = 0 + for i, tif_file in enumerate(tif_files): + try: + fbase = os.path.basename(tif_file) + + # Extract date pair(s) from the basename first, then from the parent + # directory name (e.g. Dolphin: ifgrams/YYYYMMDD_YYYYMMDD/fullres.unw.tif) + date_nums = re.findall(r'\d{8}', fbase) + if len(date_nums) < 2: + date_pair = re.findall(r'(\d{8})_(\d{8})', os.path.dirname(tif_file)) + if date_pair: + date_nums = list(date_pair[0]) + else: + date_nums = re.findall(r'\d{8}', os.path.dirname(tif_file)) + if len(date_nums) < 2: + prog_bar.update(i+1, suffix=f'skipped {i+1}/{num_file}') + continue + dates = sorted(set(date_nums))[:2] + + if not all(19000101 <= int(d) <= 20991231 for d in dates): + prog_bar.update(i+1, suffix=f'skipped {i+1}/{num_file}') + continue + + num_valid += 1 + prog_bar.update(i+1, suffix=f'{dates[0]}_{dates[1]} {i+1}/{num_file}') + + rsc_file = tif_file + '.rsc' + # regenerate a stale sidecar written by an older version of this + # script (missing PROCESSOR=isce3, e.g. carrying the buggy 90-deg + # CENTER_INCIDENCE_ANGLE), instead of propagating old values + if os.path.isfile(rsc_file): + old_processor = readfile.read_roipac_rsc(rsc_file).get('PROCESSOR') + if not update_mode or old_processor != 'isce3': + os.remove(rsc_file) + + ifg_meta = meta.copy() + ifg_meta.update(readfile.read_attribute(tif_file)) + # the declared processor from the common metadata is authoritative; + # restore it in case read_attribute detected a generic gdal product + # (e.g. Dolphin "fullres.unw.tif" before the naming-pattern detection) + if meta.get('PROCESSOR'): + ifg_meta['PROCESSOR'] = meta['PROCESSOR'] + ifg_meta = add_ifgram_metadata(ifg_meta, dates, baseline_dict) + + 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 + + prog_bar.close() + if num_valid == 0: + print(f'WARNING: no valid files processed for pattern: {obs_file}') + return + + +######################################################################### +def prep_isce3(inps): + """Prepare ISCE3/Dolphin metadata files.""" + # If no meta file is provided or it is 'auto', generate one from static_layers + if not inps.meta_file or inps.meta_file == 'auto': + print('No meta file provided. Generating burst XML from static_layers.h5...') + geom_path = Path(inps.geom_dir) + # Find all subdirectories containing static_layers*.h5 + burst_dirs = sorted([d for d in geom_path.iterdir() + if d.is_dir() and list(d.glob('static_layers*.h5'))]) + # Also check for HDF5 files directly in geom_dir + 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}') + + # Read common metadata from reference burst XML (now guaranteed to exist) + metadata = {} + if inps.meta_file: + metadata = isce3_utils.extract_isce3_metadata( + inps.meta_file, + update_mode=inps.update_mode + ) + + # Determine target shape from first interferogram if available + target_shape = None + ref_int = None + if inps.obs_files: + int_list = glob.glob(inps.obs_files[0]) + if int_list: + 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}') + + # Prepare geometry (updates metadata with LENGTH/WIDTH) + 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) + ) + + # Read baseline info + baseline_dict = {} + if inps.baseline_dir: + baseline_dict = isce3_utils.read_baseline_timeseries_isce3( + inps.baseline_dir, + processor=inps.processor + ) + + # Prepare metadata for interferogram stack(s) + 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.') + return diff --git a/src/mintpy/unwrap_error_phase_closure.py b/src/mintpy/unwrap_error_phase_closure.py index c4c871d2b..6bbcf547f 100644 --- a/src/mintpy/unwrap_error_phase_closure.py +++ b/src/mintpy/unwrap_error_phase_closure.py @@ -153,9 +153,17 @@ def calc_num_triplet_with_nonzero_integer_ambiguity(ifgram_file, mask_file=None, print_msg=False, ).reshape(num_ifgram, -1) + # keep MintPy's convention that a phase value of 0.0 indicates no-data: + # only count triplets whose 3 interferograms are all valid (non-zero), + # consistent with the data != 0. check in read_stack_obs(). + valid = unw != 0. + num_valid_leg = np.dot(np.abs(np.asarray(C)) > 0, valid.astype(np.float32)) + tri_valid = num_valid_leg == 3 + # calculate based on equation (8-9) and T_int equation inline. closure_pha = np.dot(C, unw) closure_int = np.round((closure_pha - ut.wrap(closure_pha)) / (2.*np.pi)) + closure_int[~tri_valid] = 0. num_nonzero_closure[r0:r1, :] = np.sum(closure_int != 0, axis=0).reshape(-1, width) prog_bar.update(i+1, every=1, suffix=f'line {r0} / {length}') @@ -272,9 +280,16 @@ def get_common_region_int_ambiguity(ifgram_file, cc_mask_file, water_mask_file=N print_msg=False, ).reshape(num_ifgram, -1) + # keep 0.0 == no-data: zero out closure for triplets with invalid legs + valid = unw != 0. + num_valid_leg = np.dot(np.abs(np.asarray(C)) > 0, valid.astype(np.float32)) + tri_valid = num_valid_leg == 3 + # calculate closure_int closure_pha = np.dot(C, unw) - closure_int = matrix(np.round((closure_pha - ut.wrap(closure_pha)) / (2.*np.pi))) + closure_int_arr = np.round((closure_pha - ut.wrap(closure_pha)) / (2.*np.pi)) + closure_int_arr[~tri_valid] = 0. + closure_int = matrix(closure_int_arr) # solve for U U[:,j] = np.round(l1regls( diff --git a/src/mintpy/utils/isce3_utils.py b/src/mintpy/utils/isce3_utils.py new file mode 100644 index 000000000..09456095b --- /dev/null +++ b/src/mintpy/utils/isce3_utils.py @@ -0,0 +1,1421 @@ +import glob +import os +import re +import shutil +import tempfile +import xml.dom.minidom +import xml.etree.ElementTree as ET +from collections import Counter, defaultdict +from datetime import datetime +from pathlib import Path +from typing import Any, Dict, List, Optional, Union + +# Ensure PROJ/GDAL data paths are set when running without conda activation, +# otherwise osr.ImportFromEPSG() fails with "proj_create_from_database" errors. +# This must run BEFORE "from osgeo import gdal, osr" so that the PROJ context +# is initialized with the correct data directory. +def _setup_gdal_proj_data(): + import sys + for env_key, subdir, probe in [('PROJ_DATA', 'share/proj', 'proj.db'), + ('GDAL_DATA', 'share/gdal', None)]: + if env_key in os.environ or (env_key == 'PROJ_DATA' and 'PROJ_LIB' in os.environ): + continue + data_dir = os.path.join(sys.prefix, *subdir.split('/')) + if os.path.isdir(data_dir) and (probe is None or os.path.isfile(os.path.join(data_dir, probe))): + os.environ[env_key] = data_dir + +_setup_gdal_proj_data() + +import h5py # noqa: E402 +import numpy as np # noqa: E402 +from osgeo import gdal, osr # noqa: E402 + +from mintpy.objects import sensor # noqa: E402 +from mintpy.utils import readfile # noqa: E402 + +try: + from scipy.interpolate import CubicHermiteSpline +except ImportError: + CubicHermiteSpline = None + print("Warning: scipy not available. Orbit interpolation will fall back to linear.") + + +# Default geometry layers to extract from the ISCE3/Dolphin static_layers HDF5 +# file, together with their dataset names under /data. Single source of truth +# used by cli/prep_isce3.py and prep_isce3.py as well. +GEOMETRY_FILENAMES = [ + 'height.tif', + 'los_east.tif', + 'los_north.tif', + 'layover_shadow_mask.tif', + 'local_incidence_angle.tif', +] + +GEOMETRY_DSET_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 extract_isce3_metadata(meta_file: str, update_mode: bool = True) -> dict: + """Extract common metadata from an ISCE3/Dolphin burst XML file. + + Parameters + ---------- + meta_file : str + Path to the reference burst XML file. + update_mode : bool + Not used here (kept for consistency). + + Returns + ------- + dict + Common metadata dictionary with keys required by MintPy. + + """ + # Parse XML file (Python 3.8+ disables entity expansion by default) + tree = ET.parse(meta_file) + root = tree.getroot() + burst_elem = root.find('burst_attributes') + if burst_elem is None: + raise ValueError(f'Missing in {meta_file}') + + # Helper to extract text from element + def get_value(tag): + elem = burst_elem.find(tag) + if elem is not None: + return elem.text.strip() + return None + + meta = {} + + # Basic radar parameters + meta['prf'] = get_value('prf') + meta['startUTC'] = get_value('burstStartUTC') + meta['stopUTC'] = get_value('burstStopUTC') + meta['radarWavelength'] = get_value('radarWavelength') + meta['startingRange'] = get_value('startingRange') + meta['passDirection'] = get_value('passDirection') + meta['polarization'] = get_value('polarization') + meta['trackNumber'] = get_value('trackNumber') + meta['orbitNumber'] = get_value('orbitNumber') + + # Platform name (Sentinel-1) + meta['PLATFORM'] = 'sen' + + # Sensing mid time and center line UTC + sensing_mid = get_value('sensingMid') + if sensing_mid: + try: + dt = datetime.strptime(sensing_mid, '%Y-%m-%d %H:%M:%S.%f') + seconds_of_day = dt.hour * 3600.0 + dt.minute * 60.0 + dt.second + dt.microsecond / 1e6 + meta['CENTER_LINE_UTC'] = str(seconds_of_day) + except (ValueError, TypeError): + meta['CENTER_LINE_UTC'] = '0' + else: + meta['CENTER_LINE_UTC'] = '0' + + # Pixel sizes + az_time_interval = get_value('azimuthTimeInterval') + range_pixel_size = get_value('rangePixelSize') + satellite_speed = get_value('satelliteSpeed') + + if az_time_interval and satellite_speed: + try: + az_pixel_size = float(satellite_speed) * float(az_time_interval) + meta['azimuthPixelSize'] = str(az_pixel_size) + except ValueError: + meta['azimuthPixelSize'] = '0.0' + else: + meta['azimuthPixelSize'] = '0.0' + + if range_pixel_size: + meta['rangePixelSize'] = range_pixel_size + else: + meta['rangePixelSize'] = '0.0' + + # Spatial resolution from sensor database (Sentinel-1 IW) + # Determine swath (e.g., 'IW2') + swath_num = get_value('swathNumber') + if swath_num: + iw_str = f'IW{swath_num}' + else: + # fallback: try to infer from filename + base = os.path.basename(meta_file) + if base.startswith('IW'): + iw_str = base.split('.')[0] + else: + iw_str = 'IW2' + try: + meta['azimuthResolution'] = str(sensor.SENSOR_DICT['sen'][iw_str]['azimuth_resolution']) + meta['rangeResolution'] = str(sensor.SENSOR_DICT['sen'][iw_str]['range_resolution']) + except KeyError: + meta['azimuthResolution'] = '0.0' + meta['rangeResolution'] = '0.0' + + # Heading, earth radius, altitude + meta['HEADING'] = get_value('HEADING') + meta['earthRadius'] = get_value('earthRadius') + meta['altitude'] = get_value('altitude') + + # Beam mode and swath + meta['beam_mode'] = 'IW' + meta['swathNumber'] = swath_num if swath_num else '2' + + # Frame numbers (may be 0 for ISCE3 products) + meta['firstFrameNumber'] = get_value('firstFrameNumber') or '0' + meta['lastFrameNumber'] = get_value('lastFrameNumber') or '0' + + # Antenna side (default to -1, right-looking) + meta['ANTENNA_SIDE'] = '-1' + + # Processor + meta['PROCESSOR'] = 'isce3' + + # Height / Earth radius (map from lowercase keys used in the XML) + if meta.get('altitude') and not meta.get('HEIGHT'): + meta['HEIGHT'] = meta['altitude'] + if meta.get('earthRadius') and not meta.get('EARTH_RADIUS'): + meta['EARTH_RADIUS'] = meta['earthRadius'] + + # Compute center incidence angle if possible. + # NOTE: the old formula (look_angle = asin(R/(R+H)) followed by + # asin((R+H)/R * sin(look_angle))) always collapses to exactly 90 degrees. + # Use the standard law-of-cosines solution in the triangle formed by the + # satellite, the target and the earth center, i.e. the same convention as + # ut.incidence_angle(atr, dimension=0) in utils0.py (measured from the + # local vertical; near-/mid-range slant range). + if meta.get('HEIGHT') and meta.get('EARTH_RADIUS') and meta.get('startingRange'): + try: + H = float(meta['HEIGHT']) + R = float(meta['EARTH_RADIUS']) + rho = float(meta['startingRange']) + # use the mid-range slant range when the scene width is known + if meta.get('WIDTH') and meta.get('rangePixelSize'): + rho += float(meta['WIDTH']) / 2.0 * float(meta['rangePixelSize']) + cos_inc = ((R + H) ** 2 - R ** 2 - rho ** 2) / (2.0 * R * rho) + inc_angle = np.rad2deg(np.arccos(np.clip(cos_inc, -1.0, 1.0))) + meta['CENTER_INCIDENCE_ANGLE'] = str(inc_angle) + except (ValueError, TypeError): + pass + + # Looks will be computed later from burst vs interferogram dimensions + meta['ALOOKS'] = '1' + meta['RLOOKS'] = '1' + + # Convert all values to strings (MintPy standard) + for key, value in meta.items(): + if value is not None: + meta[key] = str(value) + else: + meta[key] = '' + + # Standardize keys + meta = readfile.standardize_metadata(meta) + + # Note: LENGTH and WIDTH are not in this XML; they will be added later + # from geometry or interferogram files. + + return meta + + +def read_baseline_timeseries_isce3(baseline_dir: str, processor: str = 'tops') -> Dict: + """Read baseline time series from ISCE3/Dolphin baseline directory. + + Expected structure: baseline_dir/*.txt where each filename is YYYYMMDD_YYYYMMDD.txt + File content example: + Bperp average (m): -19.081015753493432 + Bpar average (m): -6.823039248934357 + + Parameters + ---------- + baseline_dir : str + Path to the baseline directory. + processor : str + Processor name (unused, kept for compatibility). + + Returns + ------- + dict + Dictionary of baseline values keyed by date (YYYYMMDD). + Each value is [bperp_top, bperp_bottom] (identical for both). + + """ + baseline_dict = {} + + # Expand any glob pattern in the baseline dir itself first (e.g. the + # Dolphin layout ".../baselines/t124*/" or ".../baselines/t124*/20210104/"), + # then collect YYYYMMDD_YYYYMMDD.txt files from each matched directory. + txt_files = [] + for bdir in sorted(glob.glob(baseline_dir)) or [baseline_dir]: + if not os.path.isdir(bdir): + continue + txt_files += sorted(glob.glob(os.path.join(bdir, '[0-9]*_[0-9]*.txt'))) + txt_files = sorted(set(txt_files)) + + # Fall back to a recursive search for the nested Dolphin layout, e.g. + # ".../baselines/t124_xxx/YYYYMMDD_YYYYMMDD.txt" when the baseline dir + # is given without any glob (e.g. just "../../baselines"). + if not txt_files: + for bdir in sorted(glob.glob(baseline_dir)) or [baseline_dir]: + if not os.path.isdir(bdir): + continue + txt_files += sorted( + glob.glob(os.path.join(bdir, '**', '[0-9]*_[0-9]*.txt'), + recursive=True)) + txt_files = sorted(set(txt_files)) + + if not txt_files: + print(f'WARNING: no baseline text files found in {os.path.abspath(baseline_dir)}') + return baseline_dict + + # Identify the common reference date from filenames (first part before underscore) + ref_date_candidates = [os.path.basename(f).split('_')[0] for f in txt_files] + ref_date = Counter(ref_date_candidates).most_common(1)[0][0] + + # Filter files that start with the reference date + bFiles = [f for f in txt_files if os.path.basename(f).startswith(ref_date)] + + # Read each file + for bFile in bFiles: + filename = os.path.basename(bFile) + date_pair = filename.replace('.txt', '') + dates = date_pair.split('_') + if len(dates) != 2: + continue + _, date2 = dates + + # Parse file content (robust: handles Bperp average (m), Bperp (m), bperp, etc.) + bperp = 0.0 + with open(bFile) as f: + for line in f: + line = line.strip() + match = re.match(r'[Bb]perp\b.*?:\s*([-+]?\d*\.?\d+(?:[eE][-+]?\d+)?)', line) + if match: + try: + bperp = float(match.group(1)) + except ValueError: + pass + break + + # For tops, top and bottom are assumed equal (average) + baseline_dict[date2] = [bperp, bperp] + + # Set reference date baseline to [0, 0] + baseline_dict[ref_date] = [0.0, 0.0] + + return baseline_dict + + +def extract_h5_geometry( + h5_file: Union[str, Path], + output_dir: Path, + geom_types: List[str], + dataset_mapping: Optional[Dict[str, str]] = None, + x_coords_name: str = "x_coordinates", + y_coords_name: str = "y_coordinates", + projection_name: str = "projection" +) -> Dict[str, Dict]: + """Extract geometry datasets from a static_layers HDF5 file to GeoTIFF. + + Parameters + ---------- + h5_file : Path or str + Path to the HDF5 file. + output_dir : Path + Directory to save extracted GeoTIFFs. + geom_types : list of str + Desired output filenames (e.g., ['height.tif', 'los_east.tif']). + dataset_mapping : dict, optional + Mapping from output filename to HDF5 dataset path. + x_coords_name, y_coords_name, projection_name : str + Names of coordinate and projection datasets. + + Returns + ------- + dict + Keys are geometry types, values are dicts with 'file_list' and 'nodata'. + + """ + if dataset_mapping is None: + dataset_mapping = GEOMETRY_DSET_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'] + + # 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: + wkt = proj_ds.attrs['spatial_ref'].decode('utf-8') + srs = osr.SpatialReference() + srs.ImportFromWkt(wkt) + if srs.GetAuthorityCode(None): + epsg = int(srs.GetAuthorityCode(None)) + elif proj_ds.shape == () and np.issubdtype(proj_ds.dtype, np.integer): + epsg = int(proj_ds[()]) + + # Read coordinates + 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) + + # Extract each geometry type + for geom_type in geom_types: + ds_path = dataset_mapping.get(geom_type) + if ds_path is None: + continue + if ds_path not in data_group: + continue + dataset = data_group[ds_path] + data = dataset[:] + + # Determine nodata value + nodata = None + if '_FillValue' in dataset.attrs: + nodata = dataset.attrs['_FillValue'].item() + elif np.issubdtype(data.dtype, np.floating): + nodata = np.nan + elif np.issubdtype(data.dtype, np.integer): + nodata = np.iinfo(data.dtype).max + + # Write GeoTIFF + h5_stem = h5_file.stem + file_output_dir = output_dir / h5_stem + file_output_dir.mkdir(parents=True, exist_ok=True) + out_file = file_output_dir / geom_type + + driver = gdal.GetDriverByName('GTiff') + ds_out = driver.Create( + str(out_file), + data.shape[1], data.shape[0], + 1, + gdal.GDT_Float32 if data.dtype == np.float32 else gdal.GDT_Float64, + options=['COMPRESS=LZW'] + ) + if geotransform: + ds_out.SetGeoTransform(geotransform) + try: + srs = osr.SpatialReference() + if srs.ImportFromEPSG(epsg) != 0: + raise RuntimeError(f'ImportFromEPSG({epsg}) failed') + ds_out.SetProjection(srs.ExportToWkt()) + except Exception as crs_err: + print(f'WARNING: could not set CRS EPSG:{epsg} for {geom_type} ' + f'({crs_err}). Writing GeoTIFF without projection. ' + f'Check PROJ_DATA environment variable / conda activation.') + band = ds_out.GetRasterBand(1) + band.WriteArray(data) + if nodata is not None: + band.SetNoDataValue(nodata) + band.FlushCache() + ds_out = None + + extracted[geom_type]['file_list'] = out_file + extracted[geom_type]['nodata'] = nodata + + except Exception as e: + print(f"Error processing {h5_file}: {e}") + + return extracted + +_GDAL_DTYPE_MAP = { + 'float32': 'Float32', + 'float64': 'Float64', +} + + +def build_vrt_from_h5( + h5_file: Union[str, Path], + output_dir: Path, + geom_types: List[str], + dataset_mapping: Optional[Dict[str, str]] = None, +) -> Dict[str, Dict]: + """Build tiny VRT files referencing HDF5 subdatasets directly (no raster copy). + + Fast alternative to extract_h5_geometry(): only the 1D coordinate arrays and + dataset attributes are read; the actual raster data stays in the HDF5 file and + is streamed block-wise by gdalwarp later. + + Parameters + ---------- + h5_file : Path or str + Path to the static_layers HDF5 file. + output_dir : Path + Directory to save the VRT files (one subdir per HDF5 stem). + geom_types : list of str + Desired output filenames (e.g., ['height.tif', 'los_east.tif']). + dataset_mapping : dict, optional + Mapping from output filename to HDF5 dataset name under /data. + + Returns + ------- + dict + Keys are geometry types, values are dicts with 'file_list' (VRT path) + and 'nodata'. + + """ + if dataset_mapping is None: + dataset_mapping = GEOMETRY_DSET_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 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]) + 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) + + srs = osr.SpatialReference() + if srs.ImportFromEPSG(epsg) != 0: + raise RuntimeError(f'ImportFromEPSG({epsg}) failed (check PROJ_DATA)') + srs_wkt = srs.ExportToWkt() + + vrt_dir = Path(output_dir) / h5_file.stem + vrt_dir.mkdir(parents=True, exist_ok=True) + + for geom_type in geom_types: + ds_name = dataset_mapping.get(geom_type) + if ds_name is None or ds_name not in data_group: + continue + dataset = data_group[ds_name] + length, width = dataset.shape + + # Match legacy extract_h5_geometry dtype behavior: float32 or float64 + gdal_dtype = _GDAL_DTYPE_MAP.get(dataset.dtype.name, 'Float32') + + # Determine nodata (same rules as extract_h5_geometry) + if '_FillValue' in dataset.attrs: + nodata = dataset.attrs['_FillValue'].item() + elif np.issubdtype(dataset.dtype, np.floating): + nodata = np.nan + elif np.issubdtype(dataset.dtype, np.integer): + nodata = np.iinfo(dataset.dtype).max + else: + nodata = None + + src_fname = f'HDF5:"{h5_file}"://data/{ds_name}' + nodata_xml = '' + if nodata is not None: + nodata_xml = f' {nodata}\n' + vrt_xml = ( + f'\n' + f' {srs_wkt}\n' + f' {", ".join(repr(float(v)) for v in geotransform)}\n' + f' \n' + f'{nodata_xml}' + f' \n' + f' {src_fname}\n' + f' 1\n' + f' \n' + f' \n' + f' \n' + f' \n' + f'\n' + ) + vrt_file = vrt_dir / (geom_type + '.vrt') + with open(vrt_file, 'w') as f: + f.write(vrt_xml) + + # Sanity check: VRT (and the HDF5 subdataset behind it) must be readable + ds_check = gdal.Open(str(vrt_file)) + if ds_check is None: + raise RuntimeError(f'GDAL cannot open generated VRT: {vrt_file}') + if ds_check.RasterXSize != width or ds_check.RasterYSize != length: + raise RuntimeError(f'VRT size mismatch for {vrt_file}') + ds_check = None + + extracted[geom_type]['file_list'] = vrt_file + extracted[geom_type]['nodata'] = nodata + + return extracted + + +def merge_geometry_files( + burst_ids: List[str], + geom_dir: Path, + output_dir: Path, + geom_types: List[str], + ref_int_file: Optional[str] = None, + keep_temp: bool = False +) -> Dict[str, Path]: + """Merge geometry files across multiple burst IDs and compute incidence angle. + + Parameters + ---------- + burst_ids : list of str + List of burst IDs (subdirectory names). + geom_dir : Path + Base directory containing burst subdirectories. + output_dir : Path + Output directory for merged files. + geom_types : list of str + List of geometry file basenames to merge. + ref_int_file : str, optional + Path to a reference interferogram GeoTIFF to set output grid. + keep_temp : bool + If True, keep temporary extracted files. + + Returns + ------- + dict + Mapping from geometry type to merged file path. + + """ + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + # Temporary directory for per-burst VRT (fast path) or GeoTIFF (legacy path) + # files. Created next to the output dir to avoid exhausting the system /tmp. + if keep_temp: + temp_dir = output_dir / "temp_extracted" + temp_dir.mkdir(exist_ok=True) + else: + temp_dir = Path(tempfile.mkdtemp(prefix='geom_merge_', dir=str(output_dir))) + + def _collect(use_vrt): + """Collect per-burst sources (VRT or extracted GeoTIFF) per geometry type.""" + 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")): + if use_vrt: + extracted = build_vrt_from_h5(h5_file, temp_dir, geom_types) + else: + extracted = 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 + + merged = {} + try: + # Fast path: VRTs referencing HDF5 subdatasets directly, no full-res raster + # copy. Fall back to the legacy GeoTIFF extraction on any failure. + try: + geometry_files, nodata_dict = _collect(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(use_vrt=False) + + if not any(geometry_files.values()): + print('#' * 60) + print('WARNING: no geometry layers extracted from any static_layers*.h5!') + print(' Merged geometry files will NOT be (re)generated;') + print(' stale files in the output directory may be reused downstream.') + print(' Check the "Error processing ..." messages above.') + print('#' * 60) + + # Determine output bounds and resolution from reference interferogram if provided + bounds, xres, yres, expected_size = None, None, None, None + if ref_int_file and os.path.exists(ref_int_file): + ds = gdal.Open(str(ref_int_file)) + gt = ds.GetGeoTransform() + xmin = gt[0] + ymax = gt[3] + xmax = xmin + gt[1] * ds.RasterXSize + ymin = ymax + gt[5] * ds.RasterYSize + bounds = (xmin, ymin, xmax, ymax) + xres, yres = abs(gt[1]), abs(gt[5]) + expected_size = (ds.RasterXSize, ds.RasterYSize) + ds = None + + # Merge each geometry type using gdal.Warp (python API, streams block-wise) + for gtype, file_list in geometry_files.items(): + if not file_list: + continue + out_file = output_dir / gtype + if out_file.exists(): + out_file.unlink() + + print(f"Merging {gtype} from {len(file_list)} bursts using gdal.Warp...") + warp_kwargs = dict( + format='GTiff', + creationOptions=['COMPRESS=LZW'], + resampleAlg='near', + multithread=True, + warpOptions=['NUM_THREADS=ALL_CPUS'], + ) + if bounds is not None: + warp_kwargs.update(outputBounds=bounds, xRes=xres, yRes=yres) + nodata = nodata_dict.get(gtype) + if nodata is not None: + warp_kwargs.update(dstNodata=nodata) + + ds_out = gdal.Warp(str(out_file), [str(f) for f in file_list], **warp_kwargs) + if ds_out is None: + raise RuntimeError(f'gdal.Warp failed for {gtype} ' + f'(inputs: {[str(f) for f in file_list]})') + + # Validate output grid and report coverage + out_size = (ds_out.RasterXSize, ds_out.RasterYSize) + if expected_size is not None and out_size != expected_size: + raise RuntimeError(f'merged {gtype} size {out_size} does not match ' + f'reference interferogram size {expected_size}') + arr = ds_out.GetRasterBand(1).ReadAsArray() + if nodata is not None and not (isinstance(nodata, float) and np.isnan(nodata)): + valid_ratio = float(np.mean(arr != nodata)) + else: + valid_ratio = float(np.mean(np.isfinite(arr))) + print(f' size: {out_size[0]} x {out_size[1]}, valid pixels: {valid_ratio:.1%}') + ds_out = None + merged[gtype] = out_file + + # Compute incidence angle and azimuth angle if los_east and los_north are present + los_east = merged.get('los_east.tif') + los_north = merged.get('los_north.tif') + if los_east and los_north: + inc_file = output_dir / 'incidenceAngle.tif' + compute_incidence_angle( + los_east, los_north, inc_file, + nodata=nodata_dict.get('los_east.tif') + ) + merged['incidenceAngle.tif'] = inc_file + + az_file = output_dir / 'azimuthAngle.tif' + compute_azimuth_angle( + los_east, los_north, az_file, + nodata=nodata_dict.get('los_east.tif') + ) + merged['azimuthAngle.tif'] = az_file + + return merged + + finally: + # Always clean up the temporary directory, even when the merge failed. + if not keep_temp: + shutil.rmtree(temp_dir, ignore_errors=True) + +def compute_azimuth_angle(los_east_file: Path, los_north_file: Path, output_file: Path, nodata: float = None): + """Compute azimuth angle from LOS east and north components. + + The azimuth angle is defined as the angle from the North, measured + anti‑clockwise as positive (standard mathematical convention). + + Parameters + ---------- + los_east_file, los_north_file : Path + Paths to GeoTIFF files of LOS vector components. + output_file : Path + Output path for azimuthAngle.tif. + nodata : float, optional + No-data value to use. + + """ + 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() + + # 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 + if nodata is not None: + mask = np.isnan(east) | np.isnan(north) + az_angle[mask] = nodata + + driver = gdal.GetDriverByName('GTiff') + ds_out = driver.Create(str(output_file), ds_east.RasterXSize, ds_east.RasterYSize, + 1, gdal.GDT_Float32, options=['COMPRESS=LZW']) + ds_out.SetGeoTransform(ds_east.GetGeoTransform()) + ds_out.SetProjection(ds_east.GetProjection()) + band = ds_out.GetRasterBand(1) + band.WriteArray(az_angle) + if nodata is not None: + band.SetNoDataValue(nodata) + band.FlushCache() + ds_out = None + ds_east = None + ds_north = None + +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. + + """ + 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() + + # incidence = arccos(up), where up = sqrt(1 - east^2 - north^2) + 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)) + + # Mask nodata + if nodata is not None: + mask = np.isnan(east) | np.isnan(north) + inc_angle[mask] = nodata + + driver = gdal.GetDriverByName('GTiff') + ds_out = driver.Create(str(output_file), ds_east.RasterXSize, ds_east.RasterYSize, + 1, gdal.GDT_Float32, options=['COMPRESS=LZW']) + ds_out.SetGeoTransform(ds_east.GetGeoTransform()) + ds_out.SetProjection(ds_east.GetProjection()) + band = ds_out.GetRasterBand(1) + band.WriteArray(inc_angle) + if nodata is not None: + band.SetNoDataValue(nodata) + band.FlushCache() + ds_out = None + ds_east = None + ds_north = None + + +def extract_merge_geometry( + geom_dir: str, + output_dir: str, + geom_types: List[str], + ref_int_file: Optional[str] = None, + metadata: Optional[Dict] = None, + extra_dirs: Optional[List[str]] = None +) -> Dict[str, Path]: + """High-level function to extract, merge, and prepare geometry. + + Parameters + ---------- + geom_dir : str + Base geometry directory containing burst subdirs. + output_dir : str + Output directory. + geom_types : list + Desired geometry file names. + ref_int_file : str, optional + Reference interferogram GeoTIFF for extent. + metadata : dict, optional + Metadata dictionary to be updated with LENGTH and WIDTH. + extra_dirs : list of str, optional + Additional directories containing static_layers*.h5 files + (e.g. individual burst dirs from glob expansion). + + Returns + ------- + dict + Merged geometry file paths. + + """ + geom_path = Path(geom_dir) + # Find burst IDs (subdirectories containing static_layers*.h5) + burst_ids = [] + for subdir in geom_path.iterdir(): + if subdir.is_dir() and list(subdir.glob("static_layers*.h5")): + burst_ids.append(subdir.name) + if not burst_ids: + if list(geom_path.glob("static_layers*.h5")): + burst_ids = ['.'] + + # When extra dirs are provided (multi-burst from glob expansion), + # build a temporary directory tree with symlinks so that + # merge_geometry_files can process all bursts together. + base_dir = geom_path + temp_base_dir = None + if extra_dirs: + temp_base_dir = Path(tempfile.mkdtemp()) + new_burst_ids = [] + for bid in burst_ids: + src = (geom_path / bid).resolve() + dst_name = f"_m_{bid}".replace('.', 'r') if bid == '.' else f"_m_{bid}" + os.symlink(src, temp_base_dir / dst_name, target_is_directory=True) + new_burst_ids.append(dst_name) + for i, edir in enumerate(extra_dirs): + edir_path = Path(edir) + if edir_path.is_dir() and list(edir_path.glob("static_layers*.h5")): + dst_name = f"_x_{i}" + os.symlink(edir_path.resolve(), temp_base_dir / dst_name, target_is_directory=True) + new_burst_ids.append(dst_name) + burst_ids = new_burst_ids + base_dir = temp_base_dir + + if not burst_ids: + raise FileNotFoundError(f"No static_layers HDF5 found in {geom_dir}") + + output_path = Path(output_dir) + merged = merge_geometry_files( + burst_ids=burst_ids, + geom_dir=base_dir, + output_dir=output_path, + geom_types=geom_types, + ref_int_file=ref_int_file, + keep_temp=False + ) + + if temp_base_dir: + shutil.rmtree(temp_base_dir, ignore_errors=True) + + # Update metadata with dimensions from the merged height file + if metadata is not None: + height_file = merged.get('height.tif') + if height_file and height_file.exists(): + 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']}, WIDTH={metadata['WIDTH']}") + + return merged + +############################################################################### +# XML generation for ISCE3 burst metadata (optional) +############################################################################### +def _to_seconds(t_str: str, ref_epoch_str: str) -> float: + """Convert datetime string to seconds relative to reference epoch.""" + t = datetime.strptime(t_str, '%Y-%m-%d %H:%M:%S.%f') + ref = datetime.strptime(ref_epoch_str, '%Y-%m-%d %H:%M:%S.%f') + return (t - ref).total_seconds() + + +def _compute_heading(state_vector): + """Calculate ENU heading angle from satellite state vectors. + + Parameters + ---------- + state_vector : tuple + (position, velocity) in ECEF coordinates. + position: [x, y, z] in meters. + velocity: [vx, vy, vz] in m/s. + + Returns + ------- + heading : float + ENU heading angle in degrees (0-360), clockwise from North. + + """ + # Unpack state vector + state_pos, state_vel = state_vector + + # Convert ECEF position to LLH + # Using WGS84 ellipsoid parameters + a = 6378137.0 # semi-major axis + f = 1.0 / 298.257223563 # flattening + b = a * (1 - f) # semi-minor axis + e2 = 1 - (b**2 / a**2) # eccentricity squared + + x, y, z = state_pos + + # Longitude + lon = np.arctan2(y, x) + + # Latitude using iterative method + p = np.sqrt(x**2 + y**2) + lat = np.arctan2(z, p * (1 - e2)) + + # Iterate to improve latitude accuracy + for _ in range(10): + N = a / np.sqrt(1 - e2 * np.sin(lat)**2) + h = p / np.cos(lat) - N + lat_new = np.arctan2(z, p * (1 - e2 * N / (N + h))) + if np.abs(lat_new - lat) < 1e-12: + break + lat = lat_new + + # Calculate ENU basis vectors at satellite position + sin_lat = np.sin(lat) + cos_lat = np.cos(lat) + sin_lon = np.sin(lon) + cos_lon = np.cos(lon) + + # ECEF to ENU rotation matrix + R = np.array([ + [-sin_lon, cos_lon, 0.0], + [-sin_lat * cos_lon, -sin_lat * sin_lon, cos_lat], + [cos_lat * cos_lon, cos_lat * sin_lon, sin_lat] + ]) + + # Rotate velocity from ECEF to ENU + vel_enu = np.dot(R, state_vel) + + # Extract East and North components + v_east = vel_enu[0] + v_north = vel_enu[1] + + # Calculate heading angle + heading_rad = np.arctan2(v_east, v_north) + heading_deg = np.degrees(heading_rad) + + # Normalize to 0-360 + if heading_deg < 0: + heading_deg += 360.0 + + return heading_deg + + +def _orbit_interp_hermite(metadata, time): + """Interpolate orbit state vectors at given time using Hermite interpolation. + + Parameters + ---------- + metadata : dict + Dictionary containing orbit metadata with keys: + - 'time': array of times in seconds relative to reference_epoch + - 'position_x', 'position_y', 'position_z': arrays of position components + - 'velocity_x', 'velocity_y', 'velocity_z': arrays of velocity components + - 'reference_epoch': reference epoch datetime string + time : float or array_like + Time(s) in seconds relative to reference_epoch for interpolation + + Returns + ------- + dict + Dictionary containing interpolated position and velocity at given time(s) + + """ + # Extract time array + t_array = np.array(metadata['time']) + + if CubicHermiteSpline is None: + # scipy is not available: fall back to linear interpolation of the + # position with velocity estimated by central differences. + def _pos(vals, t): + return np.interp(t, t_array, vals) + + def _vel(vals, t): + dt = np.abs(t_array[1] - t_array[0]) * 0.5 + t_lo = np.clip(np.asarray(t, dtype=float) - dt, t_array[0], t_array[-1]) + t_hi = np.clip(np.asarray(t, dtype=float) + dt, t_array[0], t_array[-1]) + return (np.interp(t_hi, t_array, vals) - np.interp(t_lo, t_array, vals)) / (2.0 * dt) + + if np.isscalar(time): + position = np.array([_pos(metadata[f'position_{c}'], time) for c in 'xyz']) + velocity = np.array([_vel(metadata[f'velocity_{c}'], time) for c in 'xyz']) + else: + time_array = np.asarray(time) + position = np.column_stack([_pos(metadata[f'position_{c}'], time_array) for c in 'xyz']) + velocity = np.column_stack([_vel(metadata[f'velocity_{c}'], time_array) for c in 'xyz']) + return position, velocity + + # Create Hermite interpolators for each component + # Position interpolation with velocity as derivatives + pos_x_interp = CubicHermiteSpline(t_array, metadata['position_x'], metadata['velocity_x']) + pos_y_interp = CubicHermiteSpline(t_array, metadata['position_y'], metadata['velocity_y']) + pos_z_interp = CubicHermiteSpline(t_array, metadata['position_z'], metadata['velocity_z']) + + # For velocity, we can either use the derivative of position interpolators + # or create separate velocity interpolators. Using derivative ensures consistency. + + # Calculate interpolated values + if np.isscalar(time): + position = np.array([ + float(pos_x_interp(time)), + float(pos_y_interp(time)), + float(pos_z_interp(time)) + ]) + # Get velocity from derivative of position interpolators + velocity = np.array([ + float(pos_x_interp.derivative()(time)), + float(pos_y_interp.derivative()(time)), + float(pos_z_interp.derivative()(time)) + ]) + else: + time_array = np.asarray(time) + position = np.column_stack([ + pos_x_interp(time_array), + pos_y_interp(time_array), + pos_z_interp(time_array) + ]) + # Get velocity from derivative of position interpolators + velocity = np.column_stack([ + pos_x_interp.derivative()(time_array), + pos_y_interp.derivative()(time_array), + pos_z_interp.derivative()(time_array) + ]) + + return position, velocity + + +def read_burst_metadata_h5( + h5_file: Path, + layer_names: List[str] = None, + group_path: str = "/metadata/processing_information/input_burst_metadata/" +) -> Dict[str, Any]: + """Read metadata for burst attributes from an HDF5 file. + + Parameters + ---------- + h5_file : Path + Path to the HDF5 file + layer_names : List[str], optional + List of burst metadata attribute names to read. + If None, all datasets directly under the group will be read. + Defaults to None. + group_path : str, optional + Path to the burst metadata group in the HDF5 file. + Defaults to '/metadata/processing_information/input_burst_metadata/' + + Returns + ------- + Dict[str, Any] + Metadata dictionary for the burst attributes + + """ + metadata = {} + + try: + with h5py.File(h5_file, 'r') as f: + # Check if group exists + if group_path not in f: + print(f"Group {group_path} not found in {h5_file}") + return metadata + + # Get the group object + group = f[group_path] + + # If layer_names is not provided, get all direct datasets in the group + if layer_names is None: + # Get all items in the group, filter only datasets (not subgroups) + layer_names = [] + for name, item in group.items(): + if isinstance(item, h5py.Dataset): + layer_names.append(name) + + # Read each requested layer + for layer_name in layer_names: + dataset_path = f"{group_path.rstrip('/')}/{layer_name}" + + if dataset_path in f: + # Get the dataset + dataset = f[dataset_path] + + # Read the value + value = dataset[()] + + # Handle different data types + if isinstance(value, np.ndarray): + # For string arrays, decode bytes to string + if value.dtype.kind == 'S' or value.dtype.kind == 'O': + if value.size == 1: + metadata[layer_name] = value.item().decode('utf-8') if isinstance(value.item(), bytes) else value.item() + else: + metadata[layer_name] = [v.decode('utf-8') if isinstance(v, bytes) else v for v in value.tolist()] + else: + # For shape array, extract length and width + if layer_name == 'shape' and value.shape == (2,): + metadata['length'] = int(value[0]) + metadata['width'] = int(value[1]) + metadata['shape'] = value.tolist() + else: + metadata[layer_name] = value.tolist() + elif isinstance(value, bytes): + metadata[layer_name] = value.decode('utf-8') + else: + # Convert scalar numpy types to Python types + metadata[layer_name] = value.item() if hasattr(value, 'item') else value + + # Optionally add dataset attributes if they exist + if dataset.attrs: + metadata[f"{layer_name}_attrs"] = { + key: (val.tolist() if isinstance(val, np.ndarray) else val) + for key, val in dataset.attrs.items() + } + else: + print(f"Dataset {dataset_path} not found in {h5_file}") + + except Exception as e: + print(f"Failed to read burst metadata from {h5_file}: {e}") + return metadata + + # Calculate sensing_mid if we have sensing_start and sensing_stop + if 'sensing_start' in metadata and 'sensing_stop' in metadata: + try: + # Parse datetime strings + start_dt = datetime.strptime(metadata['sensing_start'], '%Y-%m-%d %H:%M:%S.%f') + stop_dt = datetime.strptime(metadata['sensing_stop'], '%Y-%m-%d %H:%M:%S.%f') + + # Calculate time difference and mid time + time_diff = stop_dt - start_dt + mid_dt = start_dt + (time_diff / 2) + + # Format back to string + metadata['sensing_mid'] = mid_dt.strftime('%Y-%m-%d %H:%M:%S.%f') + + except Exception as e: + print(f"Failed to calculate sensing_mid: {e}") + # If calculation fails, try to read it directly from file + try: + dataset_path = f"{group_path.rstrip('/')}/sensing_mid" + with h5py.File(h5_file, 'r') as f: + if dataset_path in f: + value = f[dataset_path][()] + if isinstance(value, bytes): + metadata['sensing_mid'] = value.decode('utf-8') + elif hasattr(value, 'item'): + metadata['sensing_mid'] = value.item() + else: + metadata['sensing_mid'] = value + except Exception as read_e: + print(f"Failed to read sensing_mid from file: {read_e}") + + # Calculate mid_range if we have starting_range, width, and range_pixel_spacing + if all(key in metadata for key in ['starting_range', 'width', 'range_pixel_spacing']): + try: + mid_range = (metadata['starting_range'] + + (metadata['width'] / 2) * + metadata['range_pixel_spacing']) + metadata['mid_range'] = mid_range + except Exception as e: + print(f"Failed to calculate mid_range: {e}") + + return metadata + + +def prepare_mintpy_metadata(metafile: Path) -> Dict[str, Any]: + """Prepare metadata dictionary from MintPy HDF5 file for processing. + + This function extracts specific metadata groups from a MintPy HDF5 file, + combines them into a single dictionary, and extracts additional derived + metadata fields such as swath number. + + Parameters + ---------- + metafile : Path + Path to the MintPy HDF5 metadata file. + + Returns + ------- + Dict[str, Any] + Combined metadata dictionary containing: + - All metadata from specified HDF5 groups + - Derived fields like 'swathNumber' + + Notes + ----- + The function reads metadata from three predefined HDF5 group paths: + 1. Processing information and input burst metadata + 2. Orbit information + 3. Identification information + + The swath number is extracted from the 'burst_id' field using regex + pattern matching for IW1, IW2, or IW3 swaths. + + """ + # Define HDF5 group paths to extract metadata from + group_path_list = [ + '/metadata/processing_information/input_burst_metadata/', + '/metadata/orbit', + '/identification' + ] + + # Initialize empty metadata dictionary + metadata = {} + + # Iterate through each group path and read metadata + for group_path in group_path_list: + # Update metadata dictionary with contents from current group + metadata.update(read_burst_metadata_h5(metafile, group_path=group_path)) + + # Extract swath number from burst_id using regex pattern matching + # Pattern matches iw0, iw1, or iw2 (case-insensitive) and extracts the digit + if re.search(r"iw([012])", metadata.get('burst_id', ''), re.IGNORECASE): + metadata['swathNumber'] = int( + re.search(r"iw([012])", metadata['burst_id'], re.IGNORECASE).group(1) + ) + else: + metadata['swathNumber'] = None + + return metadata + + +def extract_required_attributes(metadata): + """Extract only the burst attributes needed by extract_tops_metadata.""" + isce3_available = True + try: + import isce3 + except ImportError: + isce3_available = False + print("WARNING: isce3 not available. Some metadata fields will use defaults.") + + meta = {} + + # Direct burst attributes used in original function (with fallbacks) + 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)) + # Convert slant_range_time (seconds) to meters if needed + if isinstance(meta['startingRange'], (int, float)): + try: + if meta['startingRange'] < 1.0: + SPEED_OF_LIGHT = 299792458.0 + meta['startingRange'] = meta['startingRange'] * SPEED_OF_LIGHT / 2.0 + except Exception: + meta['startingRange'] = meta.get('starting_range', 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) + + # Additional attributes needed for calculations + 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 + + # Calculate satellite speed (Vs) from orbit data (if available) + has_orbit = all(k in metadata for k in ['time', 'position_x', 'velocity_x']) + if has_orbit and isce3_available: + try: + sensingMid = _to_seconds(metadata['sensing_mid'], metadata['reference_epoch']) + sv = _orbit_interp_hermite(metadata, sensingMid) + velocity = np.linalg.norm(sv[1]) + position = sv[0] + ellipsoid = isce3.core.Ellipsoid() + # xyz_to_lon_lat() returns (lon, lat, h) in radians/meters + llh = ellipsoid.xyz_to_lon_lat(position) + heading = _compute_heading(sv) + meta['satelliteSpeed'] = velocity + meta['position'] = position + meta['HEADING'] = heading + meta['earthRadius'] = ellipsoid.r_dir(llh[1], llh[0]) + meta['altitude'] = llh[2] + except Exception as e: + print(f"WARNING: orbit interpolation failed: {e}. Using default values.") + isce3_available = False + + if not has_orbit or not isce3_available: + meta['satelliteSpeed'] = 7545.0 + meta['position'] = [0.0, 0.0, 0.0] + meta['HEADING'] = 0.0 + meta['earthRadius'] = 6371000.0 + meta['altitude'] = 693000.0 + + # Calculate frame numbers + if meta['ascendingNodeTime'] is not None and meta['burstStartUTC']: + try: + start_dt = datetime.strptime(meta['burstStartUTC'], '%Y-%m-%d %H:%M:%S.%f') + if isinstance(meta['ascendingNodeTime'], str): + node_dt = datetime.strptime(meta['ascendingNodeTime'], '%Y-%m-%d %H:%M:%S.%f') + else: + node_dt = meta['ascendingNodeTime'] + time_diff_start = (start_dt - node_dt).total_seconds() + time_diff_stop = (datetime.strptime(meta['burstStopUTC'], '%Y-%m-%d %H:%M:%S.%f') - node_dt).total_seconds() + meta['firstFrameNumber'] = int(0.2 * time_diff_start) + meta['lastFrameNumber'] = int(0.2 * time_diff_stop) + except Exception: + meta['firstFrameNumber'] = 0 + meta['lastFrameNumber'] = 0 + else: + meta['firstFrameNumber'] = 0 + meta['lastFrameNumber'] = 0 + + return meta + + +def save_burst_attributes_to_xml(burst_attrs: Dict[str, Any], output_path: Union[str, Path]) -> bool: + """Save burst attributes to XML file in MintPy-compatible format.""" + root = ET.Element('burst_metadata') + info = ET.SubElement(root, 'info') + ET.SubElement(info, 'generation_time').text = datetime.now().isoformat() + 'Z' + ET.SubElement(info, 'source').text = 'static_layers extraction' + + attrs_elem = ET.SubElement(root, 'burst_attributes') + attributes_to_save = [ + 'prf', 'burstStartUTC', 'burstStopUTC', 'radarWavelength', + 'startingRange', 'passDirection', 'polarization', 'trackNumber', + 'orbitNumber', 'sensingMid', 'azimuthTimeInterval', 'rangePixelSize', + 'swathNumber', 'ascendingNodeTime', 'satelliteSpeed', 'position', + 'HEADING', 'earthRadius', 'altitude', 'firstFrameNumber', 'lastFrameNumber' + ] + for key in attributes_to_save: + if key in burst_attrs: + elem = ET.SubElement(attrs_elem, key) + value = burst_attrs[key] + if value is None: + elem.text = 'None' + elem.set('type', 'NoneType') + elif isinstance(value, (int, np.integer)): + elem.text = str(int(value)) + elem.set('type', 'int') + elif isinstance(value, (float, np.floating)): + elem.text = f"{value:.12f}" + elem.set('type', 'float') + elif isinstance(value, str): + elem.text = value + elem.set('type', 'str') + elif isinstance(value, np.ndarray): + elem.text = str(value) + elem.set('type', 'np.ndarray') + else: + elem.text = str(value) + elem.set('type', 'other') + + xml_str = ET.tostring(root, encoding='utf-8') + dom = xml.dom.minidom.parseString(xml_str) + pretty_xml = dom.toprettyxml(indent=' ') + lines = [line for line in pretty_xml.split('\n') if line.strip()] + pretty_xml = '\n'.join(lines) + with open(output_path, 'w') as f: + f.write(pretty_xml) + return True + + +def generate_burst_xml_from_static(static_h5_file: Union[str, Path], output_xml: Union[str, Path]) -> str: + """Generate burst XML file from a single static_layers HDF5.""" + print(f'Generating burst XML from {static_h5_file}') + meta = prepare_mintpy_metadata(static_h5_file) + attrs = extract_required_attributes(meta) + save_burst_attributes_to_xml(attrs, output_xml) + print(f'Saved burst XML to {output_xml}') + return str(output_xml) diff --git a/src/mintpy/utils/readfile.py b/src/mintpy/utils/readfile.py index a41718166..c35fe704f 100644 --- a/src/mintpy/utils/readfile.py +++ b/src/mintpy/utils/readfile.py @@ -719,8 +719,8 @@ def read_binary_file(fname, datasetName=None, box=None, xstep=1, ystep=1): if 'byte order' in atr.keys() and atr['byte order'] == '0': byte_order = 'little-endian' - # GDAL / GMTSAR / ASF HyP3 - elif processor in ['gdal', 'gmtsar', 'hyp3', 'cosicorr', 'uavsar']: + # GDAL / GMTSAR / ASF HyP3 / ISCE3 geo + elif processor in ['gdal', 'gmtsar', 'hyp3', 'cosicorr', 'uavsar', 'isce3']: # try to recognize custom dataset names if specified and recognized. if datasetName: slice_list = get_slice_list(fname) @@ -738,7 +738,7 @@ def read_binary_file(fname, datasetName=None, box=None, xstep=1, ystep=1): xstep=xstep, ystep=ystep, ) - if processor in ['gdal', 'gmtsar', 'hyp3', 'cosicorr']: + if processor in ['gdal', 'gmtsar', 'hyp3', 'cosicorr', 'isce3']: data = read_gdal(fname, **kwargs) else: @@ -1275,6 +1275,22 @@ def get_hdf5_dataset(name, obj): elif fext in GDAL_FILE_EXTS: atr['PROCESSOR'] = 'gdal' + # Recognize ISCE3/Dolphin geocoded products by naming pattern. + # Dolphin may name the files "fullres.unw.tif" with the date pair + # in the parent directory name (e.g. ifgrams/20220105_20220117/). + if (fext in ['.tif', '.tiff'] + and (fbase.endswith('.unw') or fbase.endswith('.cor') + or fbase.endswith('.int') or fbase.endswith('.unw.conncomp'))): + has_date_pair = bool(re.search(r'\d{8}_\d{8}', fbase)) + if not has_date_pair: + has_date_pair = bool( + re.search(r'\d{8}_\d{8}', os.path.basename(os.path.dirname(fname)))) + if has_date_pair: + atr['PROCESSOR'] = 'isce3' + # supplement with ISCE3-specific attributes (DATA_TYPE, + # EPSG, DATE12, ...); any .rsc sidecar (read below) still + # takes priority over these values. + atr.update(read_isce3_geotiff(fname)) if 'PROCESSOR' not in atr.keys(): atr['PROCESSOR'] = 'mintpy' @@ -1322,7 +1338,9 @@ def get_hdf5_dataset(name, obj): elif meta_ext in ['.vrt'] + GDAL_FILE_EXTS: atr.update(read_gdal_vrt(metafile)) - atr['FILE_TYPE'] = fext + # keep the ISCE3-specific FILE_TYPE (e.g. '.unw' for *.unw.tif) + # set by read_isce3_geotiff() above + atr.setdefault('FILE_TYPE', fext) # DATA_TYPE for ISCE/ROI_PAC products data_type_dict = { @@ -1375,9 +1393,95 @@ def get_hdf5_dataset(name, obj): atr = standardize_metadata(atr) + # Fill in missing HEIGHT/EARTH_RADIUS for isce3 geocoded products + if atr.get('PROCESSOR', '').startswith('isce'): + # try altitude -> HEIGHT if standardization didn't catch it + if 'HEIGHT' not in atr or not atr['HEIGHT']: + for alt_key in ['altitude', 'SC_height']: + if alt_key in atr and atr[alt_key]: + atr['HEIGHT'] = str(atr[alt_key]) + break + if ('EARTH_RADIUS' not in atr or not atr['EARTH_RADIUS']) and 'earthRadius' in atr: + if atr['earthRadius']: + atr['EARTH_RADIUS'] = str(atr['earthRadius']) + return atr +def read_isce3_geotiff(fname): + """Read attributes from ISCE3/Dolphin GeoTIFF file (e.g., .int.tif / .unw.tif). + + NOTE: called from read_attribute() for ISCE3/Dolphin geocoded products to + supplement the common attributes (DATA_TYPE, EPSG, DATE12, ...). Any .rsc + sidecar file is still given priority over these values. + """ + from osgeo import gdal, osr + + ds = gdal.Open(fname, gdal.GA_ReadOnly) + if ds is None: + raise ValueError(f"Cannot open {fname} with GDAL") + + meta = {} + # Image dimensions + meta['LENGTH'] = str(ds.RasterYSize) + meta['WIDTH'] = str(ds.RasterXSize) + + # Geotransform + gt = ds.GetGeoTransform() + meta['X_FIRST'] = str(gt[0]) + meta['Y_FIRST'] = str(gt[3]) + meta['X_STEP'] = str(abs(gt[1])) + meta['Y_STEP'] = str(abs(gt[5])) + meta['X_UNIT'] = 'meters' + meta['Y_UNIT'] = 'meters' + + # Projection + proj = ds.GetProjection() + if proj: + srs = osr.SpatialReference() + srs.ImportFromWkt(proj) + if srs.IsGeographic(): + meta['X_UNIT'] = 'degrees' + meta['Y_UNIT'] = 'degrees' + epsg = srs.GetAuthorityCode(None) + if epsg: + meta['EPSG'] = epsg + + # Data type (exact numpy-style mapping, e.g. float32/float64/int16/uint8) + band = ds.GetRasterBand(1) + data_type = DATA_TYPE_GDAL2NUMPY.get(band.DataType) + if data_type: + meta['DATA_TYPE'] = data_type.replace('>', '').replace('<', '') + + # No-data value (skip when the band has none, instead of writing 'None') + ndv = band.GetNoDataValue() + if ndv is not None: + meta['NO_DATA_VALUE'] = 'nan' if np.isnan(ndv) else str(float(ndv)) + + # File type and processor + fbase = os.path.basename(fname).lower() + if fbase.endswith('.int.tif'): + meta['FILE_TYPE'] = '.int' + elif fbase.endswith('.cor.tif'): + meta['FILE_TYPE'] = '.cor' + elif 'unw' in fbase: + meta['FILE_TYPE'] = '.unw' + meta['PROCESSOR'] = 'isce3' + + # Extract DATE12 from the file name (or its parent directory, e.g. Dolphin + # "ifgrams/20220105_20220117/fullres.unw.tif") + # Pattern: YYYYMMDD_YYYYMMDD.int.tif + match = re.search(r'(\d{8})_(\d{8})', fname) + if not match: + match = re.search(r'(\d{8})_(\d{8})', os.path.basename(os.path.dirname(fname))) + if match: + d1, d2 = match.groups() + meta['DATE12'] = f"{d1[2:]}-{d2[2:]}" + + ds = None + return meta + + def auto_no_data_value(meta): """Get default no-data-value for the given file's metadata. diff --git a/src/mintpy/version.py b/src/mintpy/version.py index 0d13d295a..00d1edd46 100644 --- a/src/mintpy/version.py +++ b/src/mintpy/version.py @@ -14,6 +14,7 @@ ########################################################################### Tag = collections.namedtuple('Tag', 'version date') release_history = ( + Tag('1.6.4', '2026-07-25'), Tag('1.6.3', '2025-11-24'), Tag('1.6.2', '2025-07-07'), Tag('1.6.1', '2024-07-31'), diff --git a/src/mintpy/view.py b/src/mintpy/view.py index ef10f871d..f85ea320d 100644 --- a/src/mintpy/view.py +++ b/src/mintpy/view.py @@ -1649,19 +1649,22 @@ def plot(self): self.msk = np.ones(data.shape, dtype=np.bool_) # masking - NO_DATA_VALUE - no_data_val = readfile.get_no_data_value(self.file) - if self.no_data_value is not None: - vprint(f'masking pixels with NO_DATA_VALUE of {self.no_data_value}') - # convert integer to floating to enable masking with nan - if np.issubdtype(data.dtype, np.integer): - data = np.array(data, np.float32) - data[data == self.no_data_value] = np.nan - elif no_data_val is not None and not np.isnan(no_data_val): - vprint(f'masking pixels with NO_DATA_VALUE of {no_data_val}') - # convert integer to floating to enable masking with nan - if np.issubdtype(data.dtype, np.integer): - data = np.array(data, np.float32) - data[data == no_data_val] = np.nan + # skip for mask files: value 0 (False) is a valid display value, + # not no-data (maskConnComp/maskTempCoh/waterMask/shadowMask, etc.) + if self.key != 'mask': + no_data_val = readfile.get_no_data_value(self.file) + if self.no_data_value is not None: + vprint(f'masking pixels with NO_DATA_VALUE of {self.no_data_value}') + # convert integer to floating to enable masking with nan + if np.issubdtype(data.dtype, np.integer): + data = np.array(data, np.float32) + data[data == self.no_data_value] = np.nan + elif no_data_val is not None and not np.isnan(no_data_val): + vprint(f'masking pixels with NO_DATA_VALUE of {no_data_val}') + # convert integer to floating to enable masking with nan + if np.issubdtype(data.dtype, np.integer): + data = np.array(data, np.float32) + data[data == no_data_val] = np.nan # update/save mask info if np.ma.is_masked(data): diff --git a/tests/test_prep_isce3.py b/tests/test_prep_isce3.py new file mode 100644 index 000000000..1e9c7f068 --- /dev/null +++ b/tests/test_prep_isce3.py @@ -0,0 +1,228 @@ +"""Unit tests for the ISCE3/Dolphin interface (prep_isce3, isce3_utils, readfile).""" + +import os + +import numpy as np +import pytest +from osgeo import gdal, osr + +from mintpy.objects.stackDict import geometryDict +from mintpy.prep_isce3 import add_ifgram_metadata +from mintpy.utils import isce3_utils, readfile + + +######################################################################### +# baseline time series +######################################################################### +def _write_baseline_file(dir_path, date_pair, bperp, fmt='Bperp (m)'): + fname = os.path.join(str(dir_path), f'{date_pair}.txt') + with open(fname, 'w') as f: + f.write(f'{fmt}: {bperp}\n') + f.write('Bpar average (m): 12.345\n') + return fname + + +def test_read_baseline_timeseries_isce3(tmp_path): + # plain directory layout + _write_baseline_file(tmp_path, '20210104_20210110', -85.65046462265644) + _write_baseline_file(tmp_path, '20210104_20210116', -19.081015753493432, + fmt='Bperp average (m)') + # a pair NOT starting with the reference date must be ignored + _write_baseline_file(tmp_path, '20210110_20210116', 99.9) + + bdict = isce3_utils.read_baseline_timeseries_isce3(str(tmp_path)) + assert bdict['20210104'] == [0.0, 0.0] + assert bdict['20210110'] == pytest.approx([-85.65046462265644] * 2) + assert bdict['20210116'] == pytest.approx([-19.081015753493432] * 2) + assert len(bdict) == 3 + + +def test_read_baseline_timeseries_isce3_nested_glob(tmp_path): + # Dolphin layout: baselines/t124_xxx/YYYYMMDD_YYYYMMDD.txt + burst_dir = tmp_path / 'baselines' / 't124_264305_iw2' + burst_dir.mkdir(parents=True) + _write_baseline_file(burst_dir, '20210104_20210110', -85.65) + _write_baseline_file(burst_dir, '20210104_20210116', -19.08) + + # glob in the baseline dir itself must be expanded + bdict = isce3_utils.read_baseline_timeseries_isce3( + str(tmp_path / 'baselines' / 't124_*/')) + assert bdict['20210104'] == [0.0, 0.0] + assert bdict['20210110'] == pytest.approx([-85.65] * 2) + + # plain dir without glob must fall back to a recursive search + bdict2 = isce3_utils.read_baseline_timeseries_isce3(str(tmp_path / 'baselines')) + assert bdict2['20210110'] == pytest.approx([-85.65] * 2) + + # a non-existent dir must NOT crash, just return an empty dict + assert isce3_utils.read_baseline_timeseries_isce3( + str(tmp_path / 'nope')) == {} + + +def test_add_ifgram_metadata(): + meta = add_ifgram_metadata( + {'PROCESSOR': 'isce3'}, + dates=['20210104', '20210110'], + baseline_dict={'20210104': [0.0, 0.0], '20210110': [-85.65, -85.65]}, + ) + assert meta['DATE12'] == '210104-210110' + assert meta['P_BASELINE_TOP_HDR'] == '-85.65' + assert meta['P_BASELINE_BOTTOM_HDR'] == '-85.65' + + +######################################################################### +# burst XML metadata extraction (CENTER_INCIDENCE_ANGLE regression) +######################################################################### +def _write_burst_xml(path, attrs): + isce3_utils.save_burst_attributes_to_xml(attrs, str(path)) + return path + + +def test_extract_isce3_metadata_center_incidence_angle(tmp_path): + # realistic S1 IW2 values (Hawaii T124, from an actual dolphin run) + attrs = { + 'prf': 1451.627112193990, + 'burstStartUTC': '2021-01-04 04:30:18.346686', + 'burstStopUTC': '2021-01-04 04:30:21.430020', + 'radarWavelength': 0.05546576, + 'startingRange': 845516.028030508198, + 'passDirection': 'Ascending', + 'polarization': 'VV', + 'trackNumber': 124, + 'orbitNumber': 25000, + 'sensingMid': '2021-01-04 04:30:19.888353', + 'azimuthTimeInterval': 0.0020555563, + 'rangePixelSize': 2.329562114715, + 'swathNumber': 2, + 'satelliteSpeed': 7597.883725532828, + 'HEADING': 347.728905505086, + 'earthRadius': 6346500.407216606, + 'altitude': 697692.716961, + } + xml_file = _write_burst_xml(tmp_path / 'IW2.burst.xml', attrs) + meta = isce3_utils.extract_isce3_metadata(str(xml_file)) + + # the old formula returned exactly ~90 deg; the fixed one must be the + # realistic near-range incidence angle (~36 deg for these values) + inc = float(meta['CENTER_INCIDENCE_ANGLE']) + assert 30.0 < inc < 45.0 + assert abs(inc - 90.0) > 45.0 + + # sanity of a few other keys + assert meta['PROCESSOR'] == 'isce3' + assert meta['PLATFORM'] == 'sen' + assert float(meta['CENTER_LINE_UTC']) == pytest.approx(4 * 3600 + 30 * 60 + 19.888353) + assert float(meta['HEIGHT']) == pytest.approx(697692.716961) + + +def test_read_baseline_missing_dir_warns(tmp_path, capsys): + isce3_utils.read_baseline_timeseries_isce3(str(tmp_path / 'empty')) + assert 'no baseline text files found' in capsys.readouterr().out + + +######################################################################### +# read_attribute on ISCE3/Dolphin GeoTIFF products +######################################################################### +def _write_synthetic_geotiff(fname, shape=(8, 10), gt=(500000.0, 10.0, 0.0, 4000000.0, 0.0, -10.0), + epsg=32605, nodata=0.0): + driver = gdal.GetDriverByName('GTiff') + ds = driver.Create(str(fname), shape[1], shape[0], 1, gdal.GDT_Float32) + ds.SetGeoTransform(gt) + srs = osr.SpatialReference() + srs.ImportFromEPSG(epsg) + ds.SetProjection(srs.ExportToWkt()) + band = ds.GetRasterBand(1) + data = np.ones(shape, dtype=np.float32) * 1.5 + data[0, :] = nodata + band.WriteArray(data) + if nodata is not None: + band.SetNoDataValue(nodata) + ds = None + return fname + + +def test_read_attribute_isce3_unw_tif(tmp_path): + # Dolphin-style name with the date pair in the filename + fname = str(tmp_path / '20220105_20220117.unw.tif') + _write_synthetic_geotiff(fname) + atr = readfile.read_attribute(fname) + assert atr['PROCESSOR'] == 'isce3' + assert atr['DATE12'] == '220105-220117' + assert atr['FILE_TYPE'] == '.unw' + assert atr['LENGTH'] == 8 + assert atr['WIDTH'] == 10 + assert atr['EPSG'] == '32605' + assert atr['NO_DATA_VALUE'] == '0.0' + + +def test_read_attribute_isce3_fullres_parent_dir(tmp_path): + # Dolphin "fullres.unw.tif" layout: date pair only in the parent dir name + ifg_dir = tmp_path / '20220105_20220117' + ifg_dir.mkdir() + fname = str(ifg_dir / 'fullres.unw.tif') + _write_synthetic_geotiff(fname) + atr = readfile.read_attribute(fname) + assert atr['PROCESSOR'] == 'isce3' + assert atr['DATE12'] == '220105-220117' + assert atr['FILE_TYPE'] == '.unw' + + +def test_read_attribute_isce3_rsc_priority(tmp_path): + # a .rsc sidecar written by prep_isce3 must take priority over the + # attributes derived from the GeoTIFF itself + fname = str(tmp_path / '20220105_20220117.unw.tif') + _write_synthetic_geotiff(fname) + with open(fname + '.rsc', 'w') as f: + f.write('PROCESSOR isce3\n') + f.write('DATE12 220105-220117\n') + f.write('P_BASELINE_TOP_HDR -2.5048807779724456\n') + f.write('P_BASELINE_BOTTOM_HDR -2.5048807779724456\n') + + atr = readfile.read_attribute(fname) + assert atr['PROCESSOR'] == 'isce3' + assert atr['P_BASELINE_TOP_HDR'] == '-2.5048807779724456' + assert 'FILE_PATH' in atr + + +def test_read_isce3_geotiff_no_nodata(tmp_path): + # a band without nodata must not produce the string 'None' + fname = str(tmp_path / '20220105_20220117.int.tif') + _write_synthetic_geotiff(fname, nodata=None) + atr = readfile.read_attribute(fname) + assert atr['PROCESSOR'] == 'isce3' + assert 'NO_DATA_VALUE' not in atr or atr['NO_DATA_VALUE'] != 'None' + + +######################################################################### +# water mask auto-align (geometryDict._warp_water_mask) +######################################################################### +def test_warp_water_mask(tmp_path): + # reference grid: 10 x 8 pixels at 10 m + ref_fname = str(tmp_path / 'ref.tif') + _write_synthetic_geotiff(ref_fname, shape=(8, 10)) + + # water mask on a coarser grid: 5 x 4 pixels at 20 m, same CRS/extent + src_fname = str(tmp_path / 'water_mask.tif') + driver = gdal.GetDriverByName('GTiff') + ds = driver.Create(src_fname, 5, 4, 1, gdal.GDT_Float32) + ds.SetGeoTransform((500000.0, 20.0, 0.0, 4000000.0, 0.0, -20.0)) + srs = osr.SpatialReference() + srs.ImportFromEPSG(32605) + ds.SetProjection(srs.ExportToWkt()) + band = ds.GetRasterBand(1) + mask = np.ones((4, 5), dtype=np.float32) # 1 = land + mask[0, :] = 0.0 # 0 = water + band.WriteArray(mask) + ds = None + + geom_obj = geometryDict( + processor='isce3', + datasetDict={'height': ref_fname, 'waterMask': src_fname}, + extraMetadata={'dummy': 'x'}, + ) + result = geom_obj._warp_water_mask('waterMask', 8, 10) + assert result.shape == (8, 10) + assert set(np.unique(result)).issubset({0.0, 1.0}) + # the top row of the source (water) stays water after the warp + assert np.all(result[0, :] == 0.0) + assert np.all(result[-1, :] == 1.0)