From 22822a282cf11756f41742c37bbe662de7f99a5f Mon Sep 17 00:00:00 2001 From: william Wu Date: Wed, 26 Aug 2026 17:17:58 +0800 Subject: [PATCH 1/2] Added Stack processing for Nisar L1 RSLC products, handles pair wise dense offsets and inverstions in the isce2 style using isce3 cuda optimized modules. Setup instruction, run and configure manuals are included as well. --- bin/geo2rdr.py | 1 - bin/rdr2geo.py | 1 - .../packages/isce3/atmosphere/tec_product.py | 1 - .../packages/isce3/cal/point_target_info.py | 2 - .../isce3/core/rdr_geo_block_generator.py | 1 - python/packages/isce3/focus/valid_regions.py | 1 - python/packages/isce3/io/background.py | 3 - python/packages/nisar/antenna/beamformer.py | 1 - python/packages/nisar/antenna/pattern.py | 1 - .../nisar/cal/corner_reflector_slc_func.py | 1 - .../nisar/cal/faraday_rotation_angle_slc.py | 1 - .../nisar/cal/pol_channel_imbalance_slc.py | 1 - .../nisar/products/insar/granule_id.py | 1 - .../products/readers/GenericProduct/GCOV.py | 1 - .../readers/instrument/instrument_parser.py | 1 - .../nisar/products/readers/rslc_cal.py | 2 +- python/packages/nisar/products/writers/SLC.py | 1 - python/packages/nisar/workflows/focus.py | 1 - python/packages/nisar/workflows/geogrid.py | 1 - share/stack/.gitignore | 1 + share/stack/README.md | 133 + share/stack/ScanList.py | 167 ++ share/stack/StackProcessor.py | 243 ++ share/stack/coarse_alignment.py | 195 ++ share/stack/configuration-guide.md | 78 + share/stack/dense-offset-and-inversion.md | 204 ++ share/stack/insar.yaml | 2503 +++++++++++++++++ share/stack/invert_offsets.py | 884 ++++++ share/stack/make_stack.py | 122 + share/stack/offset_util.py | 265 ++ share/stack/stack_utils.py | 161 ++ .../extensions/pybind/geometry/rdr2geo.py | 1 - .../python/extensions/pybind/unwrap/phass.py | 1 - tests/python/packages/estimate_dop.py | 1 - tests/python/packages/isce3/core/gpu_check.py | 1 - .../packages/isce3/io/gdal/gdal_raster.py | 1 - tools/imagesets/imgset.py | 8 +- tools/run.py | 2 +- 38 files changed, 4965 insertions(+), 29 deletions(-) create mode 100644 share/stack/.gitignore create mode 100644 share/stack/README.md create mode 100644 share/stack/ScanList.py create mode 100644 share/stack/StackProcessor.py create mode 100644 share/stack/coarse_alignment.py create mode 100644 share/stack/configuration-guide.md create mode 100644 share/stack/dense-offset-and-inversion.md create mode 100755 share/stack/insar.yaml create mode 100644 share/stack/invert_offsets.py create mode 100644 share/stack/make_stack.py create mode 100644 share/stack/offset_util.py create mode 100644 share/stack/stack_utils.py diff --git a/bin/geo2rdr.py b/bin/geo2rdr.py index f10d7589be..d4e6e963f5 100755 --- a/bin/geo2rdr.py +++ b/bin/geo2rdr.py @@ -3,7 +3,6 @@ import numpy as np import argparse from osgeo import gdal -import sys import os import isce3 diff --git a/bin/rdr2geo.py b/bin/rdr2geo.py index 802a6043cc..771bafdc6a 100755 --- a/bin/rdr2geo.py +++ b/bin/rdr2geo.py @@ -3,7 +3,6 @@ import numpy as np import argparse from osgeo import gdal -import sys import os import isce3 from nisar.products.readers import SLC diff --git a/python/packages/isce3/atmosphere/tec_product.py b/python/packages/isce3/atmosphere/tec_product.py index 9c239e2bcd..77bdb5258a 100644 --- a/python/packages/isce3/atmosphere/tec_product.py +++ b/python/packages/isce3/atmosphere/tec_product.py @@ -1,7 +1,6 @@ ''' Package to compute TEC LUT from JSON file ''' -from datetime import datetime, timedelta import json import os diff --git a/python/packages/isce3/cal/point_target_info.py b/python/packages/isce3/cal/point_target_info.py index 9f13dfd90b..d73838564d 100755 --- a/python/packages/isce3/cal/point_target_info.py +++ b/python/packages/isce3/cal/point_target_info.py @@ -21,12 +21,10 @@ class MissingNull(Exception): """Raised when mainlobe null(s) cannot be determined""" - pass class UnsupportedWindow(Exception): """Raised if window_type input is not supported.""" - pass def get_chip(x: DatasetReader, i: float, j: float, nchip: int = 64) -> np.ndarray: """ diff --git a/python/packages/isce3/core/rdr_geo_block_generator.py b/python/packages/isce3/core/rdr_geo_block_generator.py index fd6bf220ec..1bb0abf97c 100644 --- a/python/packages/isce3/core/rdr_geo_block_generator.py +++ b/python/packages/isce3/core/rdr_geo_block_generator.py @@ -1,6 +1,5 @@ import journal import numpy as np -import os import isce3 diff --git a/python/packages/isce3/focus/valid_regions.py b/python/packages/isce3/focus/valid_regions.py index 269a1634a4..301c6cac9b 100644 --- a/python/packages/isce3/focus/valid_regions.py +++ b/python/packages/isce3/focus/valid_regions.py @@ -2,7 +2,6 @@ from .sar_duration import get_sar_duration from copy import deepcopy from dataclasses import dataclass -from functools import reduce import isce3 import numpy as np import shapely diff --git a/python/packages/isce3/io/background.py b/python/packages/isce3/io/background.py index 9d2e72cf95..d4c75e36a4 100644 --- a/python/packages/isce3/io/background.py +++ b/python/packages/isce3/io/background.py @@ -72,7 +72,6 @@ def process(self, *args, **kw): """ User-defined task to operate in background thread. """ - pass def queue_work(self, *args, **kw): """ @@ -146,7 +145,6 @@ def write(self, *args, **kw): """ User-defined method for writing data. """ - pass class BackgroundReader(BackgroundWorker): @@ -199,4 +197,3 @@ def read(self, *args, **kw): """ User-defined method for reading a chunk of data. """ - pass diff --git a/python/packages/nisar/antenna/beamformer.py b/python/packages/nisar/antenna/beamformer.py index 4acc7cf9bd..e7f04a7cd4 100644 --- a/python/packages/nisar/antenna/beamformer.py +++ b/python/packages/nisar/antenna/beamformer.py @@ -18,7 +18,6 @@ class BadHPACalWarning(Warning): """Warning for Bad HPA CAL values""" - pass class ElevationBeamformer(ABC): diff --git a/python/packages/nisar/antenna/pattern.py b/python/packages/nisar/antenna/pattern.py index e338da5706..193370e06b 100644 --- a/python/packages/nisar/antenna/pattern.py +++ b/python/packages/nisar/antenna/pattern.py @@ -1,4 +1,3 @@ -from warnings import warn from collections import defaultdict from isce3.core import Orbit, Attitude, Linspace from isce3.geometry import DEMInterpolator diff --git a/python/packages/nisar/cal/corner_reflector_slc_func.py b/python/packages/nisar/cal/corner_reflector_slc_func.py index 9a07fcea56..bae2c6c764 100644 --- a/python/packages/nisar/cal/corner_reflector_slc_func.py +++ b/python/packages/nisar/cal/corner_reflector_slc_func.py @@ -46,7 +46,6 @@ class CRInfoSlc: class OutOfSlcBoundWarning(UserWarning): """Warning when CR LLH is out of SLC azimuth-range bound""" - pass # functions diff --git a/python/packages/nisar/cal/faraday_rotation_angle_slc.py b/python/packages/nisar/cal/faraday_rotation_angle_slc.py index aa8c187f56..647dce76a6 100644 --- a/python/packages/nisar/cal/faraday_rotation_angle_slc.py +++ b/python/packages/nisar/cal/faraday_rotation_angle_slc.py @@ -448,7 +448,6 @@ def estimate(self, azt_blk_size=5.0, sr_blk_size=3000.0, Too small block sizes in either range or azimuth. """ - pass class FaradayRotEstBickelBates(FaradayRotAngleSlc): diff --git a/python/packages/nisar/cal/pol_channel_imbalance_slc.py b/python/packages/nisar/cal/pol_channel_imbalance_slc.py index 7446c899dc..7f23d0de4b 100644 --- a/python/packages/nisar/cal/pol_channel_imbalance_slc.py +++ b/python/packages/nisar/cal/pol_channel_imbalance_slc.py @@ -18,7 +18,6 @@ class OutOfSlcBoundError(Exception): """Azimuth time or slant range values are out of SLC data boundary""" - pass @dataclass(frozen=True) diff --git a/python/packages/nisar/products/insar/granule_id.py b/python/packages/nisar/products/insar/granule_id.py index 87e50ba890..33305d736a 100644 --- a/python/packages/nisar/products/insar/granule_id.py +++ b/python/packages/nisar/products/insar/granule_id.py @@ -1,4 +1,3 @@ -from datetime import datetime import h5py import journal diff --git a/python/packages/nisar/products/readers/GenericProduct/GCOV.py b/python/packages/nisar/products/readers/GenericProduct/GCOV.py index 42cc777bb8..bea02ec4d8 100644 --- a/python/packages/nisar/products/readers/GenericProduct/GCOV.py +++ b/python/packages/nisar/products/readers/GenericProduct/GCOV.py @@ -1,7 +1,6 @@ # -*- coding: utf-8 -*- from __future__ import annotations -import os import h5py import pyre diff --git a/python/packages/nisar/products/readers/instrument/instrument_parser.py b/python/packages/nisar/products/readers/instrument/instrument_parser.py index a26eaa6e43..d7bf07dc14 100644 --- a/python/packages/nisar/products/readers/instrument/instrument_parser.py +++ b/python/packages/nisar/products/readers/instrument/instrument_parser.py @@ -15,7 +15,6 @@ class MissingInstrumentFieldWarning(UserWarning): """Used for missing field in instrument table """ - pass class InstrumentParser: diff --git a/python/packages/nisar/products/readers/rslc_cal.py b/python/packages/nisar/products/readers/rslc_cal.py index 34598cce53..e2b92d678d 100644 --- a/python/packages/nisar/products/readers/rslc_cal.py +++ b/python/packages/nisar/products/readers/rslc_cal.py @@ -6,7 +6,7 @@ import logging from numpy import angle, deg2rad, rad2deg, exp, pi from ruamel.yaml import YAML -from typing import List, Optional, Tuple, Union +from typing import Optional, Tuple, Union log = logging.getLogger("rslc_cal") diff --git a/python/packages/nisar/products/writers/SLC.py b/python/packages/nisar/products/writers/SLC.py index 257cceeae2..576eee6440 100644 --- a/python/packages/nisar/products/writers/SLC.py +++ b/python/packages/nisar/products/writers/SLC.py @@ -4,7 +4,6 @@ import numpy as np from numpy.linalg import norm from numpy.testing import assert_allclose -import os from shapely import Polygon from typing import Optional import isce3 diff --git a/python/packages/nisar/workflows/focus.py b/python/packages/nisar/workflows/focus.py index 11fc1009c0..087075fcee 100644 --- a/python/packages/nisar/workflows/focus.py +++ b/python/packages/nisar/workflows/focus.py @@ -5,7 +5,6 @@ from functools import reduce import h5py from itertools import chain -import json import logging import math import os diff --git a/python/packages/nisar/workflows/geogrid.py b/python/packages/nisar/workflows/geogrid.py index 2c5057b89a..133273724b 100644 --- a/python/packages/nisar/workflows/geogrid.py +++ b/python/packages/nisar/workflows/geogrid.py @@ -1,7 +1,6 @@ ''' collection of functions for determining and setting geogrid ''' -import copy import numpy as np from osgeo import osr diff --git a/share/stack/.gitignore b/share/stack/.gitignore new file mode 100644 index 0000000000..7b2511e445 --- /dev/null +++ b/share/stack/.gitignore @@ -0,0 +1 @@ +pyproject.toml \ No newline at end of file diff --git a/share/stack/README.md b/share/stack/README.md new file mode 100644 index 0000000000..af09fc595e --- /dev/null +++ b/share/stack/README.md @@ -0,0 +1,133 @@ +# ISCE3 Stack Processing + +A working first draft to align a stack of **NISAR scans** to a reference. Currently works for single GPU setups. + +## Preperation + +1. Regular ISCE DEM files are needed with the following postfix: + +``` +.dem.wgs84 + +.dem.wgs84.vrt + +.dem.wgs84.xml +``` + + A convenient way to produce them is the [dem-stitcher tool for isce workflow](https://github.com/ACCESS-Cloud-Based-InSAR/dem-stitcher/blob/dev/notebooks/Staging_a_DEM_for_ISCE2.ipynb). + +2. Download the L1 RSLC h5 product for the scans you wish to process and put them in the same folder. The h5 files should be in the original naming convention:#### + +**Example:** + +```bash +NISAR_L1_PR_RSLC_004_122_D_067_4005_DHDH_A_20251106T160541_20251106T160622_X05009_N_F_J_001.h5 +``` + +**Example aligned to format:** + +```bash +NISAR_L1_PR_RSLC_________ _____.h5 +NISAR_L1_PR_RSLC_ 004 _ 122 _ D _ 067 _ 4005 _ DHDH _ A _20251106T160541_20251106T160622_X05009_ N _ F _ J _ 001 .h5 +``` + + + +Stack alignment will be performed on all scans for a specific **Track,Frame and Direction*** combination. + +If multiple combinations are present, they will be processed as separate stacks. + +## Aligning the stacks + +Run the following command to start the alignment workflow: + +```python +python run_stack.py \ + --slc-folder \ + --out-folder \ + --dem-file +``` + +You can also run the workflow from a custom python script as shown below: + +```python +from StackProcessor import StackProcessor +processor = self = StackProcessor( + slc_folder='', + out_folder='', + dem_file='', + ) +processor.coarse_register_scans() +processor.dense_offset_pairs() +processor.invert_dense_offset_pairs() +processor.rubbersheet() +processor.resample_slcs() +``` + +To do this you would need to add the stack folder to python path by doing: + +```bash +export ISCE_STACK={full_path_to_your_share/stack} +export PYTHONPATH=${PYTHONPATH}:${ISCE_STACK} +``` + +For detailed config of the cmd and script way of runing the pipleine, see [configuration guide](configuration-guide.md) + +## Alignment Steps + +the scans in the stack is aligned using the following steps, again a stack is all the scans found in the input folder with the same track, frame and direction. By default the ** **earliest scan** is used as the reference + +- coarse_register_scans + + - geometry of the reference is extracted using ```isce3.cuda.geometry.Rdr2Geo``` + + - coarse offsets for each secondary scans are calculated with ```isce3.cuda.geometry.Geo2Rdr``` + + - secondary h5s are than resampled to coarse aligned slcs with ```nisar.workflows.resample_slc_v2.resample_secondary_rslc_onto_reference``` + +- dense_offset_pairs + + - dense offsets are calculated between date pairs using ```isce3.cuda.matchtemplate.PyCuAmpcor``` with configs provided to ```StackProcessor.dense_offset_pairs```. See more details about which pairs are choosen [here](./dense-offset-and-inversion.md#about-dense-offset-pairs) + +- invert_dense_offset_pairs + + - Offsets are then inverted to obtain offset to the reference scan. + This process prevents files with large offsets from skewing offsets fir the rest of the stack. You can find more detail [here](./dense-offset-and-inversion.md#dense-offset-inversion) + +- rubbersheet + + - The inverted offsets are then upscaled to full resoltuion through interpolation + +- resample_slcs + + - the secondary h5s are then sampled to aligned slcs that can be used for downstream processing with Mintpy or dolphin + +## Output Folders + +This sections descrbes all the output folders produced by the pipeline. + +Each stack is stored in a directory identified by its track, frame, and orbit direction, for example: + +```textile +output_folder/ +├── processor.pkl # a pickle of the processor object +├── t_f_/ +├── tXX_fXX_A/ +├── tXX_fXX_D/ +``` + +A processed stack contains the following directories: + +```textile +tXX_fXX_A/ +├── ref_geom/ [coarse] rdr2geo reference radar geometry +├── coarse_offsets/ [coarse] geo2rdr range/azimuth offsets +├── stack/ [coarse] coarsely registered SLC stack +├── dense_pairwise/ [dense] pairwise Ampcor dense offsets +├── inverted_offsets/ [invert] inverted offsets to reference SLC +├── h5_stack/ [rubbersheet] by product of rubbersheet +├── rubber/ [rubbersheet] final rubbersheeted offsets +├── merged/ [resample] FINAL aligned/resampled SLC stack +``` + + diff --git a/share/stack/ScanList.py b/share/stack/ScanList.py new file mode 100644 index 0000000000..fdd70b4b38 --- /dev/null +++ b/share/stack/ScanList.py @@ -0,0 +1,167 @@ + + +import numpy as np +import yaml +from pathlib import Path +import numpy as np +import yaml +import os +import re +from datetime import datetime +import isce3 +from nisar.products.readers import RSLC + + +def baseline(ref_file, sec_file, freq="A"): + ref = RSLC(hdf5file=ref_file) + sec = RSLC(hdf5file=sec_file) + + rg1 = ref.getRadarGrid(freq) + rg2 = sec.getRadarGrid(freq) + + orb1 = ref.getOrbit().copy() + orb2 = sec.getOrbit().copy() + + # Make radar-grid times and orbit times use identical epochs + orb1.update_reference_epoch(rg1.ref_epoch) + orb2.update_reference_epoch(rg2.ref_epoch) + + # Center of reference radar grid + t1 = rg1.sensing_mid + r1 = rg1.mid_range + + ellipsoid = isce3.core.Ellipsoid() + + # Ground point observed at reference scene center + llh = isce3.geometry.rdr2geo( + t1, + r1, + orb1, + rg1.lookside, + doppler=0.0, + wavelength=rg1.wavelength, + ellipsoid=ellipsoid, + ) + + xyz = np.asarray( + ellipsoid.lon_lat_to_xyz(llh), + dtype=float + ) + + # Find where same ground point occurs in secondary acquisition + zero_doppler = isce3.core.LUT2d() + + t2, r2 = isce3.geometry.geo2rdr( + llh, + ellipsoid, + orb2, + zero_doppler, + rg2.wavelength, + rg2.lookside, + ) + + p1, v1 = orb1.interpolate(t1) + p2, v2 = orb2.interpolate(t2) + + p1 = np.asarray(p1) + p2 = np.asarray(p2) + v1 = np.asarray(v1) + + B = p2 - p1 + + # Unit LOS: reference satellite -> ground + los = xyz - p1 + los /= np.linalg.norm(los) + + # Along-track direction + along = v1 / np.linalg.norm(v1) + + # Cross-track direction perpendicular to LOS + perp = np.cross(along, los) + perp /= np.linalg.norm(perp) + + Bpar = np.dot(B, los) + Bperp = np.dot(B, perp) + Balong = np.dot(B, along) + + return { + "B": np.linalg.norm(B), + "Bperp": Bperp, + "Bpar": Bpar, + "Balong": Balong, + "t_ref": t1, + "t_sec": t2, + "r_ref": r1, + "r_sec": r2, + } + + +def grep_date6(fn): + m = re.search(r'20\d{6}', os.path.basename(fn)) + return m.group() if m else os.path.basename(fn) + + +def grep_date8(s): + m = re.search(r'\d{8}', s) + if m: + date = m.group() + return date + + +class Scan: + def __init__(self, path): + self.path = Path(path) + self.path_str = str(path) + self.stem = self.path.stem + self.date_str = grep_date8(self.stem) + self.date = datetime.strptime(self.date_str, "%Y%m%d") + # self.scan_path = [f't{}_' self.stem.split('_')[4:6]] + parts = Path(self.stem).stem.split("_") + self.cycle, self.track, self.direction, self.frame, \ + self.mode, self.polarization = int(parts[4]), int(parts[5]), \ + parts[6], int(parts[7]), parts[8], parts[9] + self.scan_path_id = f't{self.track}_f{self.frame}_{self.direction}' + + def __str__(self): + return f"Scan {self.stem}" + + def __repr__(self): + return f"Scan {self.stem}" + + +def print_baseline_choices(reference, secondaries): + ref = reference.path_str + print("REFERENCE:", grep_date6(ref)) + for i in secondaries: + sec = i.path_str + print(i.scan_path_id) + b = baseline(ref, sec) + + print( + f"{grep_date6(sec):10s} " + f"{b['B']:12.2f} " + f"{b['Bperp']:12.2f} " + f"{b['Bpar']:14.2f} " + f"{b['Balong']:12.2f}" + ) + + +class ScanList: + def __init__(self, files): + self.files = files + self.scans = [Scan(i) for i in files] + self.scan_path_ids = np.unique([i.scan_path_id for i in self.scans]) + self.scan_by_id = {} + self.ref_scans = {} + for idi in self.scan_path_ids: + scan_id_list = [i for i in self.scans if i.scan_path_id == idi] + sort_id = np.argsort([i.date for i in scan_id_list]) + self.scan_by_id[idi] = [scan_id_list[i] for i in sort_id] + self.ref_scans[idi] = scan_id_list[sort_id[0]] + + def iter_path_ref_sec(self, path): + val = self.scan_by_id[path] + for i in range(len(val) - 1): + reference = val[i] + secondary = val[i + 1] + yield reference, secondary diff --git a/share/stack/StackProcessor.py b/share/stack/StackProcessor.py new file mode 100644 index 0000000000..9d8d7837b2 --- /dev/null +++ b/share/stack/StackProcessor.py @@ -0,0 +1,243 @@ +from nisar.workflows.rubbersheet import run_rubbersheet_with_interpolation +import os +from pathlib import Path +import shutil +import pickle + +from coarse_alignment import align_secondary, extract_ref_topo +from invert_offsets import invert_offsets_isce3 +from nisar.workflows.helpers import copy_raster +from offset_util import (create_minimal_rifg_for_rubbersheet, + dense_offset_coregistered) +from ScanList import ScanList +from stack_utils import resample, relative_symlink_contents, set_nested, load_config + +pol_freq = {'A': ['HH']} +runcfg = load_config() +cfg = runcfg.cfg["runconfig"]["groups"] +dense_cfg = cfg["processing"]["dense_offsets"] +rdr2geo_cfg = cfg['processing']['rdr2geo'] +rubbersheet_cfg = cfg['processing']['rubbersheet'] + + +class StackProcessor(ScanList): + def __init__( + self, + slc_folder, + out_folder, + dem_file, + pair_level=2, + freq='A', + pol='HH', + overwrite=False, + debug=False): + self.slc_folder = Path(slc_folder) + self.out_folder = Path(out_folder) + self.config_save = self.out_folder / 'processor.pkl' + self.dem_file = Path(dem_file) + self.overwrite = overwrite + self.freq = freq + self.pol = pol + self.debug = debug + h5_files = [str(i) for i in self.slc_folder.glob('*1.h5')] + super().__init__(h5_files) + self.pol_freq = pol_freq + self.pair_level = pair_level + self.init_pairs() + if not os.path.exists(self.config_save) or self.overwrite: + pickle.dump(self, open(self.config_save, 'wb')) + + def init_pairs(self): + self.dense_pairs = {} + for path, scans in self.scan_by_id.items(): + nscans = len(scans) + path_pairs = [] + for i in range(self.pair_level): + if nscans - i > 1: + for j in range(nscans - i - 1): + path_pairs.append((scans[j], scans[j + i + 1])) + self.dense_pairs[path] = path_pairs + + def iter_path_dense_pairs(self, path): + yield from self.dense_pairs[path] + + def coarse_register_scans(self, rdr2geo_cfg=rdr2geo_cfg): + for path in self.scan_by_id: + main_ref = self.ref_scans[path] + path_folder = self.out_folder / Path(path) + stack_folder = path_folder / Path('stack') + ref_geom_folder = path_folder / Path('ref_geom') + coarse_offset_folder = path_folder / Path('coarse_offsets') + os.makedirs(path_folder, exist_ok=True) + os.makedirs(stack_folder, exist_ok=True) + os.makedirs(coarse_offset_folder, exist_ok=True) + reference_slc = stack_folder / Path(f'{main_ref.date_str}.slc') + if not os.path.exists(ref_geom_folder) or self.overwrite: + os.makedirs(ref_geom_folder, exist_ok=True) + extract_ref_topo( + main_ref.path_str, + self.dem_file, + rdr2geo_cfg, + ref_geom_folder) + if not os.path.exists(reference_slc) or self.overwrite: + copy_raster( + main_ref.path_str, + self.freq, # frequency + self.pol, # polarization + 1024, # lines per block + reference_slc, + file_type="ENVI", + ) + pairs = {} + for ref, sec in self.iter_path_ref_sec(path): + coarse_offset_folderi = coarse_offset_folder / \ + Path(sec.date_str) + val = {} + val['azimuth'] = coarse_offset_folderi / \ + Path(f'geo2rdr/freq{self.freq}/azimuth.off') + val['range'] = coarse_offset_folderi / \ + Path(f'geo2rdr/freq{self.freq}/range.off') + pairs[(ref.date_str, sec.date_str)] = val + if not os.path.exists(coarse_offset_folderi) or self.overwrite: + os.makedirs(coarse_offset_folderi, exist_ok=True) + align_secondary( + sec.path_str, + self.dem_file, + ref_geom_folder, + coarse_offset_folderi) + + sec_slc = stack_folder / Path(f'{sec.date_str}.slc') + if not os.path.exists(sec_slc) or self.overwrite: + print(f'calculating coarse reg for {path} {sec.date_str}') + offset_dir = coarse_offset_folderi / "geo2rdr" / "freqA" + resample( + stack_folder, + offset_dir, + main_ref.path_str, + sec.path_str, + sec.date_str) + + def dense_offset_pairs(self, dense_cfg): + for path in self.scan_by_id: + for ref, sec in self.iter_path_dense_pairs(path): + path_folder = self.out_folder / Path(path) + dense_folder = path_folder / Path('dense_pairwise') + stack_folder = path_folder / Path('stack') + os.makedirs(dense_folder, exist_ok=True) + starting_pair = f'{ref.date_str}_{sec.date_str}' + out_dir = dense_folder / Path(starting_pair) + if not os.path.exists(out_dir) or self.overwrite: + print(f'calculating dense offset {path} {starting_pair}') + os.makedirs(out_dir, exist_ok=True) + reference_slc = stack_folder / Path(f'{ref.date_str}.slc') + secondary_slc = stack_folder / Path(f"{sec.date_str}.slc") + dense_offset_coregistered( + reference_slc, + secondary_slc, + out_dir, + dense_cfg=dense_cfg, + gpu_id=0, + ) + + def invert_dense_offset_pairs(self): + for path in self.scan_by_id: + pairs = {} + path_folder = self.out_folder / Path(path) + dense_folder = path_folder / Path('dense_pairwise') + inverted_offsets_folder = path_folder / Path('inverted_offsets') + if not os.path.exists(inverted_offsets_folder) or self.debug: + os.makedirs(inverted_offsets_folder, exist_ok=True) + for ref, sec in self.iter_path_dense_pairs(path): + starting_pair = f'{ref.date_str}_{sec.date_str}' + out_dir = dense_folder / Path(starting_pair) + pairs[(ref.date_str, sec.date_str) + ] = out_dir / 'dense_offsets' + if len(pairs) > 0: + outputs = invert_offsets_isce3( + pair_offsets=pairs, + output_dir=str(inverted_offsets_folder), + reference_date=self.ref_scans[path].date_str, + ) + + def rubbersheet(self, rubbersheet_cfg=rubbersheet_cfg): + for path in self.scan_by_id: + path_folder = self.out_folder / Path(path) + coarse_offset_folder = path_folder / Path('coarse_offsets') + inverted_offsets_folder = path_folder / Path('inverted_offsets') + dense_folder = path_folder / Path('dense_pairwise') + rubber_folder = path_folder / Path('rubber') + h5_stack_folder = path_folder / Path('h5_stack') + os.makedirs(h5_stack_folder, exist_ok=True) + os.makedirs(rubber_folder, exist_ok=True) + main_ref = self.ref_scans[path] + for ref, sec in self.iter_path_ref_sec(path): + coarse_offset_folderi = coarse_offset_folder / \ + Path(sec.date_str) + out_dir = h5_stack_folder / sec.date_str + rubberi = rubber_folder / sec.date_str + if not os.path.exists(rubberi) or self.overwrite: + print(f'calculating rubbersheet {path} {sec.date_str}') + cfg['processing']['rubbersheet'] = rubbersheet_cfg + inverted_offset_pathi = inverted_offsets_folder / sec.date_str + starting_pair = f'{ref.date_str}_{sec.date_str}' + densei = dense_folder / Path(starting_pair) + link_to = path_folder / Path('dense_offsets/freqA/HH') + relative_symlink_contents(densei, link_to) + relative_symlink_contents(inverted_offset_pathi, link_to) + cfg['dynamic_ancillary_file_group']['dem_file'] = str( + self.dem_file) + set_nested( + cfg, [ + 'input_file_group', 'reference_rslc_file'], main_ref.path_str) + set_nested( + cfg, [ + 'product_path_group', 'sas_output_file'], out_dir) + set_nested( + cfg, [ + 'product_path_group', 'scratch_path'], path_folder) + cfg['processing']['rubbersheet']['geo2rdr_offsets_path'] = coarse_offset_folderi + cfg['processing']['rubbersheet']['dense_offsets_path'] = path_folder + cfg['processing']['input_subset']['list_of_frequencies'] = { + self.freq: [self.pol]} + create_minimal_rifg_for_rubbersheet( + out_dir, cfg) + run_rubbersheet_with_interpolation( + cfg, out_dir) + shutil.move( + path_folder / + Path('rubbersheet_offsets'), + rubberi) + + def resample_slcs(self): + for path in self.scan_by_id: + path_folder = self.out_folder / Path(path) + merged_folder = path_folder / Path('merged') + rubber_folder = path_folder / Path('rubber') + os.makedirs(merged_folder, exist_ok=True) + main_ref = self.ref_scans[path] + ref_slc_folder = merged_folder / main_ref.date_str + reference_slc = ref_slc_folder / Path(f'{main_ref.date_str}.slc') + if not os.path.exists(reference_slc) or self.overwrite: + print(f'copying reference slc for {path}') + os.makedirs(ref_slc_folder, exist_ok=True) + copy_raster( + main_ref.path_str, + self.freq, # frequency + self.pol, # polarization + 1024, # lines per block + reference_slc, + file_type="ENVI", + ) + for _, sec in self.iter_path_ref_sec(path): + print(f'calculating final offset for {path} {sec.date_str}') + rubberi = rubber_folder / sec.date_str / \ + f'freq{self.freq}' / self.pol + pairs_out = merged_folder / sec.date_str + if not os.path.exists(pairs_out) or self.overwrite: + print(f'resampling {path} {sec.date_str}') + resample( + pairs_out, + rubberi, + main_ref.path_str, + sec.path_str, + out_tag=sec.date_str) diff --git a/share/stack/coarse_alignment.py b/share/stack/coarse_alignment.py new file mode 100644 index 0000000000..040a705232 --- /dev/null +++ b/share/stack/coarse_alignment.py @@ -0,0 +1,195 @@ +#!/usr/bin/env python3 + +''' +wrapper for rdr2geo +''' + +import pathlib +import time + +from osgeo import gdal + +import journal +import isce3 +from isce3.core import crop_external_orbit +from nisar.products.readers import SLC +from nisar.products.readers.orbit import load_orbit_from_xml +from nisar.workflows.rdr2geo_runconfig import Rdr2geoRunConfig +from nisar.workflows.yaml_argparse import YamlArgparse + +def get_raster_obj(out_path: str, radargrid: isce3.product.RadarGridParameters, + write2disk: bool, dtype: int) -> None: + '''Function that returns io.Raster or None based on write2disk bool + + dtype has to be a GDAL datatype + ''' + if not write2disk: + return None + + return isce3.io.Raster(out_path, radargrid.width, radargrid.length, 1, + dtype, 'ENVI') + +def align_secondary(sec_hdf5,dem_file,topo_path,out_path,freqs=['A'],threshold = 1e-8, + numiter = 25, lines_per_block = 1000,sec_orbit = None,gpu_enabled = True, + gpu_id = 0): + ''' + run rdr2geo + ''' + out_path = pathlib.Path(out_path) + slc = SLC(hdf5file=sec_hdf5) + + # Get orbit + orbit = slc.getOrbit() + if sec_orbit is not None: + # SLC will get first radar grid whose frequency is available. + # Reference epoch and orbit have no frequency dependency. + external_orbit = load_orbit_from_xml(sec_orbit, slc.getRadarGrid().ref_epoch) + orbit = crop_external_orbit(external_orbit, orbit) + + # Set ellipsoid based on DEM epsg + dem_raster = isce3.io.Raster(str(dem_file)) + epsg = dem_raster.get_epsg() + proj = isce3.core.make_projection(epsg) + ellipsoid = proj.ellipsoid + + # NISAR RSLC products are always zero doppler + doppler_grid = isce3.core.LUT2d() + + info_channel = journal.info('geo2rdr.run') + info_channel.log("starting geo2rdr") + + # check if gpu use if required + use_gpu = isce3.core.gpu_check.use_gpu(gpu_enabled,gpu_id) + + if use_gpu: + # set CUDA device + device = isce3.cuda.core.Device(gpu_id) + isce3.cuda.core.set_device(device) + + t_all = time.time() + + for freq in freqs: + + # Get parameters specific for that frequency + radar_grid = slc.getRadarGrid(frequency=freq) + + # Create geo2rdr directory + geo2rdr_scratch_path = out_path / 'geo2rdr' / f'freq{freq}' + geo2rdr_scratch_path.mkdir(parents=True, exist_ok=True) + + # Initialize CPU or GPU geo2rdr object accordingly + if use_gpu: + Geo2Rdr = isce3.cuda.geometry.Geo2Rdr + else: + Geo2Rdr = isce3.geometry.Geo2Rdr + + geo2rdr_obj = Geo2Rdr(radar_grid, orbit, ellipsoid, doppler_grid, + threshold, numiter, lines_per_block) + + # Open Topo Raster + topo_path = pathlib.Path(topo_path) + rdr2geo_topo_path = topo_path / 'rdr2geo' / f'freq{freq}' / 'topo.vrt' + topo_raster = isce3.io.Raster(str(rdr2geo_topo_path)) + + # Run geo2rdr + geo2rdr_obj.geo2rdr(topo_raster, str(geo2rdr_scratch_path)) + + t_all_elapsed = time.time() - t_all + info_channel.log(f"Successfully ran geo2rdr in {t_all_elapsed:.3f} seconds") + + +def extract_ref_topo(input_hdf5,dem_file,rdr2geo_cfg,out_path,freqs=['A'],threshold = 1e-8, + numiter = 25, extraiter = 10 ,lines_per_block = 1000,ref_orbit = None, + gpu_enabled = True,gpu_id = 0): + + + out_path = pathlib.Path(out_path) + slc = SLC(hdf5file=input_hdf5) + + # Get orbit + orbit = slc.getOrbit() + if ref_orbit is not None: + # SLC will get first radar grid whose frequency is available. + # Reference epoch and orbit have no frequency dependency. + external_orbit = load_orbit_from_xml(ref_orbit, slc.getRadarGrid().ref_epoch) + orbit = crop_external_orbit(external_orbit, orbit) + + # set defaults shared by both frequencies + dem_raster = isce3.io.Raster(str(dem_file)) + epsg = dem_raster.get_epsg() + proj = isce3.core.make_projection(epsg) + ellipsoid = proj.ellipsoid + + # NISAR RSLC products are always zero doppler + grid_doppler = isce3.core.LUT2d() + + info_channel = journal.info("rdr2geo.run") + info_channel.log("starting rdr2geo") + + # check if gpu ok to use + use_gpu = isce3.core.gpu_check.use_gpu(gpu_enabled,gpu_id) + if use_gpu: + # Set the current CUDA device. + device = isce3.cuda.core.Device(gpu_id) + isce3.cuda.core.set_device(device) + + t_all = time.time() + for freq in freqs: + # get frequency specific parameters + radargrid = slc.getRadarGrid(freq) + + # create separate directory within scratch dir for rdr2geo run + rdr2geo_scratch_path = out_path / 'rdr2geo' / f'freq{freq}' + rdr2geo_scratch_path.mkdir(parents=True, exist_ok=True) + + # init CPU or CUDA object accordingly + if use_gpu: + Rdr2Geo = isce3.cuda.geometry.Rdr2Geo + else: + Rdr2Geo = isce3.geometry.Rdr2Geo + + rdr2geo_obj = Rdr2Geo(radargrid, orbit, ellipsoid, grid_doppler, + threshold=threshold, numiter=numiter, + extraiter=extraiter, + lines_per_block=lines_per_block) + + # dict of layer names keys to tuples of their output name and GDAL types + layers = {'x':('x', gdal.GDT_Float64), 'y':('y', gdal.GDT_Float64), + 'z':('z', gdal.GDT_Float64), + 'incidence':('incidence', gdal.GDT_Float32), + 'heading':('heading', gdal.GDT_Float32), + 'local_incidence':('localIncidence', gdal.GDT_Float32), + 'local_psi':('localPsi', gdal.GDT_Float32), + 'simulated_amplitude':('simamp', gdal.GDT_Float32), + 'layover_shadow':('layoverShadowMask', gdal.GDT_Byte)} + + # rdr2geo_cfg = cfg['processing']['rdr2geo'] + + # list comprehend rasters to be written from layers dict + raster_list = [ + get_raster_obj(f'{str(rdr2geo_scratch_path)}/{fname}.rdr', + radargrid, rdr2geo_cfg[f'write_{key_name}'], + dtype) + for key_name, (fname, dtype) in layers.items()] + + # extract individual elements from dict as args for topo + x_raster, y_raster, height_raster, incidence_raster,\ + heading_raster, local_incidence_raster, local_psi_raster,\ + simulated_amplitude_raster, shadow_raster = raster_list + + # run topo - with east and north unit vector components of ground to + # satellite layers permanently disabled. + rdr2geo_obj.topo(dem_raster, x_raster, y_raster, height_raster, + incidence_raster, heading_raster, local_incidence_raster, + local_psi_raster, simulated_amplitude_raster, + shadow_raster, None, None) + + # remove undesired/None rasters from raster list + raster_list = [raster for raster in raster_list if raster is not None] + + # save non-None rasters to vrt + output_vrt = isce3.io.Raster(f'{str(rdr2geo_scratch_path)}/topo.vrt', raster_list) + output_vrt.set_epsg(rdr2geo_obj.epsg_out) + + t_all_elapsed = time.time() - t_all + info_channel.log(f"successfully ran rdr2geo in {t_all_elapsed:.3f} seconds") \ No newline at end of file diff --git a/share/stack/configuration-guide.md b/share/stack/configuration-guide.md new file mode 100644 index 0000000000..5f95417881 --- /dev/null +++ b/share/stack/configuration-guide.md @@ -0,0 +1,78 @@ +# Configuration Guide + +This is a guide for all the options for configuration of the cmd line and script version of running the pipeline + +## Options for the commandline tool + +```bash +usage: run_stack.py --slc-folder SLC_FOLDER --out-folder OUT_FOLDER + --dem-file DEM_FILE [--pair-level PAIR_LEVEL] + [--freq FREQ] [--pol POL] [--overwrite] [--debug] + [--steps {coarse,dense,invert,rubbersheet,resample} + [{coarse,dense,invert,rubbersheet,resample} ...]] + +Run NISAR StackProcessor pipeline. + +options: + + --slc-folder SLC_FOLDER + Folder containing NISAR RSLC HDF5 files. + + --out-folder OUT_FOLDER + Output stack folder. + + --dem-file DEM_FILE DEM file used by the stack processing workflow. + + --pair-level PAIR_LEVEL + Dense-offset pair level. Default: 2. + + --freq FREQ Frequency to process. Default: A. + + --pol POL Polarization to process. Default: HH. + + --overwrite Overwrite existing outputs. + + --steps {coarse,dense,invert,rubbersheet,resample} ... + Pipeline stages to run. If omitted, all stages run in + this order: coarse dense invert rubbersheet resample. +``` + + + +## Options for the script pipeline + +The optional inputs above can also be supplied to the script version of the pipeline: + +```python +processor = self = StackProcessor( + slc_folder='', + out_folder='', + dem_file='', + + # Optional arguments + pair_level=2, # dense-offset pairing level + # 1 = adjacent acquisitions only + # 2 = first- and second-neighbor pairs + freq='A', # NISAR frequency + pol='HH', # polarization + overwrite=False, # overwrite existing outputs +) +``` + +the ```rdr2geo``` ,```dense_offsets``` and```rubbersheet``` steps of the workflow can be configured with isce3 format. The default is as defined in nisar's example [insar.yaml](https://github.com/isce-framework/isce3/blob/develop/share/nisar/defaults/insar.yaml). + +example + +```python +from stack_utils import load_config + +runcfg = load_config('') +cfg = runcfg.cfg["runconfig"]["groups"] +dense_cfg = cfg["processing"]["dense_offsets"] +rdr2geo_cfg = cfg['processing']['rdr2geo'] +rubbersheet_cfg = cfg['processing']['rubbersheet'] + +processor.coarse_register_scans(rdr2geo_cfg=rdr2geo_cfg) +processor.dense_offset_pairs(dense_cfg=dense_cfg) +processor.rubbersheet(rubbersheet_cfg=rubbersheet_cfg) +``` diff --git a/share/stack/dense-offset-and-inversion.md b/share/stack/dense-offset-and-inversion.md new file mode 100644 index 0000000000..c33506c2ed --- /dev/null +++ b/share/stack/dense-offset-and-inversion.md @@ -0,0 +1,204 @@ +## About Dense Offset Pairs + +Dense offsets are calculated between selected pairs of scans in the stack rather than between every possible combination of acquisitions. + +This is done as a systematic way to form pair networks for equally spaces nisar scans + +For a chronologically ordered stack: + +$$ +S_0,\;S_1,\;S_2,\;S_3,\;S_4 +$$ + +the parameter `pair_level` controls how far apart in the scan sequence two acquisitions may be when forming a dense-offset pair. + +### `pair_level = 1` + +Only adjacent scans are paired: + +$$ +(S_0,S_1),\; +(S_1,S_2),\; +(S_2,S_3),\; +(S_3,S_4) +$$ + +This gives the minimum set of pairwise measurements required to connect the stack. In this case, each scan offset is propagated directly along the date chain, and inversion provides no additional correction. + +### `pair_level = 2` + +Both adjacent scans and second-neighbor scans are paired: + +$$ +(S_0,S_1),\; +(S_1,S_2),\; +(S_2,S_3),\; +(S_3,S_4) +$$ + +and: + +$$ +(S_0,S_2),\; +(S_1,S_3),\; +(S_2,S_4) +$$ + +The additional pairs introduce redundant measurements into the scan network. + +For example, the offset between $S_0$ and $S_2$ can now be constrained in two ways: + +$$ +S_0 \rightarrow S_1 \rightarrow S_2 +$$ + +and directly by: + +$$ +S_0 \rightarrow S_2 +$$ + +Ideally, + +$$ +d_{01}+d_{12}=d_{02} +$$ + +but in practice the dense-offset measurements contain correlation noise, so these values will not be exactly equal. + +The redundant pair measurements reduce the influence of Ampcor fluctuations and measurement errors by jointly solving all available pairwise offsets to estimate a consistent offset for each scan relative to the stack reference. + +## Dense Offset Inversion + +The inversion process essentially derives, for each secondary scan, a dense offset to the reference that provides the **least-squares fit** to all calculated pairwise offset. + +Suppose the stack contains four scans: + +$$ +S_0,\;S_1,\;S_2,\;S_3 +$$ + +and $S_0$ is chosen as the reference scan. + +Dense-offset matching does not initially give the offset of every scan relative to $S_0$. Instead, it gives offsets for selected scan pairs, for example: + +$$ +d_{01},\quad d_{12},\quad d_{23},\quad d_{02},\quad d_{13} +$$ + +where $d_{ij}$ is the measured range or azimuth offset between scans $S_i$ and $S_j$. + +Define the unknown offsets relative to the reference as: + +$$ +x_0=0,\quad x_1,\quad x_2,\quad x_3 +$$ + +Each pairwise measurement gives one linear constraint: + +$$ +x_j-x_i=d_{ij} +$$ + +For the example above: + +$$ +\begin{aligned} +x_1-x_0 &= d_{01} \\ +x_2-x_1 &= d_{12} \\ +x_3-x_2 &= d_{23} \\ +x_2-x_0 &= d_{02} \\ +x_3-x_1 &= d_{13} +\end{aligned} +$$ + +Since $x_0=0$, this becomes: + +$$ +\begin{bmatrix} +1 & 0 & 0 \\ +-1 & 1 & 0 \\ +0 & -1 & 1 \\ +0 & 1 & 0 \\ +-1 & 0 & 1 +\end{bmatrix} +\begin{bmatrix} +x_1 \\ +x_2 \\ +x_3 +\end{bmatrix} += +\begin{bmatrix} +d_{01} \\ +d_{12} \\ +d_{23} \\ +d_{02} \\ +d_{13} +\end{bmatrix} +$$ + +In compact form: + +$$ +Ax=d +$$ + +where: + +- $A$ describes the scan-pair network. +- $d$ contains the measured pairwise dense offsets. +- $x$ contains the unknown offset of each scan relative to the reference scan. + +Because several scan pairs constrain the same unknown offsets, the system is usually overdetermined. The solution is therefore obtained in a least-squares sense: + +$$ +\hat{x} += +\arg\min_x \|Ax-d\|_2^2 +$$ + +This is useful because the pairwise measurements are generally not perfectly consistent due to correlation noise. + +For example: + +$$ +d_{01}+d_{12} +$$ + +may not exactly equal: + +$$ +d_{02} +$$ + +Even though both describe the displacement between $S_0$ and $S_2$. + +The inversion finds the set of per-scan offsets that best satisfies all available pairwise measurements simultaneously. + +For a standard least-squares problem, the solution can be written as: + +$$ +\hat{x} += +(A^T A)^{-1}A^T d +$$ + +when $A^T A$ is invertible. + +The same inversion is performed independently at every pixel in the dense-offset grid and separately for the two offset components: + +- range offset +- azimuth offset + +The output therefore gives, for each scan $S_i$, a range and azimuth offset field relative to the reference scan $S_0$: + +$$ +S_i +\rightarrow +\begin{cases} +\text{range offset relative to } S_0 \\ +\text{azimuth offset relative to } S_0 +\end{cases} +$$ + +These inverted offset fields are then combined with the coarse alignment during rubbersheeting and are used to produce the final resampled SLC stack. diff --git a/share/stack/insar.yaml b/share/stack/insar.yaml new file mode 100755 index 0000000000..c019819973 --- /dev/null +++ b/share/stack/insar.yaml @@ -0,0 +1,2503 @@ +runconfig: + name: gunw_workflow_default + + groups: + + logging: + path: insar.log + write_mode: w + + pge_name_group: + pge_name: INSAR_L_PGE + + input_file_group: + # Required - File path to reference RSLC HDF5 file + reference_rslc_file: ... + # Required - File path to secondary RSLC HDF5 file + secondary_rslc_file: ... + + # REQUIRED for QA. NOT REQUIRED if only running Product SAS. + # These field(s) should contain the filename of the + # corresponding NISAR product(s) for QA to process. + # If Product SAS and QA SAS are run back-to-back, + # the(se) field(s) should be identical to filename convention + # defined in `sas_output_file`. + qa_rifg_input_file: RIFG_product.h5 + qa_runw_input_file: RUNW_product.h5 + qa_gunw_input_file: GUNW_product.h5 + qa_roff_input_file: ROFF_product.h5 + qa_goff_input_file: GOFF_product.h5 + + dynamic_ancillary_file_group: + # REQUIRED - Digital elevation model + dem_file: ... + + # Digital elevation model file description, optional + dem_file_description: + + water_mask_file: + + # Water mask file description, optional + water_mask_file_description: + + orbit_files: + reference_orbit_file: + secondary_orbit_file: + tec_file: + + # Troposphere weather model files + troposphere_weather_model_files: + reference_troposphere_file: + secondary_troposphere_file: + + product_path_group: + # Directory where PGE will place results. Irrelevant to SAS. + product_path: product_path + + # Directory where SAS can write temporary data + scratch_path: . + + # SAS writes output product to the following file. PGE may rename. + # After the SAS completes, the PGE wrapper renames the product file + # according to proper file naming conventions. + # In INSAR mode, GUNW, GOFF, RIFG, ROFF, and RUNW will be prepended + # to sas_output_file + sas_output_file: _product.h5 + + # REQUIRED for QA. NOT REQUIRED if only running Product SAS. + # Output directory to store all QA output files. + # Because multiple InSAR products can be generated by this + # runconfig, QA will make new product-specific directories + # here to store the corresponding output files in. + # Ex: if the output dir is set to './qa' and QA is requested for + # RIFG and RUNW products, then QA will create './qa/rifg' and + # './qa/runw' and save the respective QA outputs in these. + qa_output_dir: ./qa + + primary_executable: + product_type: RIFG + product_version: + rifg_version: + roff_version: + runw_version: + goff_version: + gunw_version: + product_doi: + rifg_doi: + roff_doi: + runw_doi: + goff_doi: + gunw_doi: + static_layers_data_access: + debug_level_group: + debug_switch: false + + # TODO OPTIONAL - To setup type of worker + worker: + # OPTIONAL - To remove the intermediate files. Default False + intermediate_files_removal_enabled: False + + # OPTIONAL - To prevent downloading DEM / other data automatically. Default True + internet_access: False + + # OPTIONAL - To explicitly use GPU capability if available. Default False + gpu_enabled: True + + # Index of the GPU to use for processing, optional. Defaults to the + # first available CUDA device. Ignored if *gpu_enabled* is False. + gpu_id: 0 + + # ADT section - isce3 + pyre workflow + processing: + # Process single co-pol in offset and rubbersheeting + process_single_co_pol_offset: True + + input_subset: + # Frequencies and polarizations to be processed + list_of_frequencies: + # Keys for frequency A and B followed by a list of polarizations to process: + # {frequency}: [polarizations] + # + # - Only frequencies listed (i.e, 'A' and/or 'B') will be processed. + # - If the [polarizations] list is empty, all available polarizations + # for given frequency will be processed + A: + - HH + B: + - HH + + # TODO OPTIONAL - Only checked when internet access is available + dem_download: + # OPTIONAL - s3 bucket / curl URL / local file + source: + top_left: + x: + y: + bottom_right: + x: + y: + + # OPTIONAL - Mechanism to specify output posting and DEM + geocode: + # Optional - path to RUNW product to geocode + runw_path: + # Optional - path to ROFF product to geocode + roff_path: + # Optional - path to RIFG product to geocode + rifg_path: + + # OPTIONAL - To control output grid in same units as output EPSG + x_snap: + + # OPTIONAL - To control output grid in same units as output EPSG + y_snap: + + # OPTIONAL - Same as input DEM if not provided. + output_epsg: + + # OPTIONAL - Spacing between pixels, in same units as output EPSG. + # If not provided, spacing values will be taken from provided DEM. + # All postings/spacings must be > 0. + # ISCE3 output rasters always have North-up West-left orientation + output_posting: + A: + x_posting: + y_posting: + B: + x_posting: + y_posting: + + # OPTIONAL - Spacing between pixels for wrapped interferogram, + # in same units as output EPSG. + # If not provided, spacing values will be taken from provided DEM. + # All postings/spacings must be > 0. + # ISCE3 output rasters always have North-up West-left orientation + wrapped_interferogram: + # Default interpolation method is SINC method + interp_method: SINC + # OPTIONAL - To control output grid in same units as output EPSG + x_snap: + # OPTIONAL - To control output grid in same units as output EPSG + y_snap: + + output_posting: + A: + x_posting: + y_posting: + B: + x_posting: + y_posting: + + # OPTIONAL - Can control with absolute values or with snap values + top_left: + # OPTIONAL - Set top-left y in same units as output EPSG + y_abs: + # OPTIONAL - Set top-left x in same units as output EPSG + x_abs: + + # OPTIONAL - Can control with absolute values or with snap values + bottom_right: + y_abs: + x_abs: + + # OPTIONAL - Set lines to be processed per block + lines_per_block: 1000 + + # OPTIONAL - Set default data interpolation method to bilinear + interp_method: BILINEAR + + radar_grid_cubes: + + # List of heights in meters, where the parallelBaseline, perpendicularBaseline, and slantRangeSolidEarthTidesPhase + # datacubes contain ONLY TWO height layers (the bottom and top heights), i.e., -500 and 9000 from the default heights below + heights: [-500, 0, 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500, 7000, 7500, 8000, 8500, 9000] + + # OPTIONAL - Same as the geocode group output_epsg if not provided + output_epsg: + + # OPTIONAL - Spacing between pixels, in same units as output EPSG. + # If not provided, spacing values are 0.005 degrees for EPSG:4326 + # and spacing values of other EPSG are 500 meters. + # All postings/spacings must be > 0. + # ISCE3 output rasters always have North-up West-left orientation + output_posting: + x_posting: + y_posting: + + # OPTIONAL - To control output grid in same units as output EPSG + y_snap: + + # OPTIONAL - To control output grid in same units as output EPSG + x_snap: + + # OPTIONAL - Can control with absolute values or with snap values + top_left: + # OPTIONAL - Set top-left y in same units as output EPSG + y_abs: + # OPTIONAL - Set top-left x in same units as output EPSG + x_abs: + + # OPTIONAL - Can control with absolute values or with snap values + bottom_right: + y_abs: + x_abs: + + bandpass: + # Number of lines per blocks + lines_per_block: 1000 + # Set window function for range spectrum bandpass filtering + window_function: tukey + # Shape parameter. For Kaiser, 0 <= shape < Inf. + # For Cosine and Tukey, 0 <= shape <= 1. + window_shape: 0.25 + # FFT size in range. + range_fft_size: + + ionosphere_phase_correction: + enabled: False + # Number of lines per blocks + lines_per_block: 1000 + # ionosphere estimation methods + #(split_main_band, main_side_band, main_diff_ms_band, + # main_diff_low_high_subband) + spectral_diversity: main_diff_low_high_subband + # Polarization to be used for ionosphere estimation. + # Split_main_band method does not require frequency B. + list_of_frequencies: + A: + - HH + B: + - HH + # Options for split_main_band method + split_range_spectrum: + lines_per_block: 2048 + # Bandwidth for lower sub-band SLC in Hz + # If unspecified, it is automatically set by the workflow + low_band_bandwidth: + # Bandwidth for upper sub-band SLC in Hz + # If unspecified, it is automatically set by the workflow + high_band_bandwidth: + window_function: tukey + # Shape parameter. For Kaiser, 0 <= shape < Inf. + # For Cosine and Tukey, 0 <= shape <= 1. + window_shape: 0.25 + # Options for dispersive filter + dispersive_filter: + # Apply filtering to dispersive layer + enabled: True + # Filter type for dispersive raster + # (median_filter, coherence, connected_components) + # When multiple validation sources are enabled, + # take the intersection of their invalid masks + # to form a combined mask. + filter_mask_type: ['coherence'] + # Median filter size for mask generation + median_filter_size: 15 + median_filter_threshold: 0.6 + # Coherence thresholds 0<= thresholds <=1 + filter_coherence_threshold: 0.5 + # Filling method for invalid regions + # (nearest, smoothed) + filling_method: nearest + # Method used to build the smooth guide image for outlier detection + # and gap filling. {'median_gaussian', 'gaussian', 'none'}, optional + filling_guide_filter_method: median_gaussian + # Window size of the NaN-aware median filter used when + # `guide_filter_method='median_gaussian'` + filling_guide_median_size: 3 + # Threshold multiplier for robust outlier detection. Pixels are marked + # as outliers when + # abs(residual - median(residual)) > + # outlier_threshold * robust_scale + # where `robust_scale = max(1.4826 * MAD, outlier_min_scale)`. + filling_outlier_threshold: 3.5 + # Lower bound for the robust scale used in outlier detection. + # This prevents the threshold from becoming too small in very + # clean regions, + filling_outlier_min_scale: 0.0 + # Scale factor that converts MAD (median absolute deviation) + # into a sigma-like robust spread estimate for outlier detection. + # For a standard normal variable Z ~ N(0, 1), + # median(|Z|) = Phi^{-1}(0.75) ≈ 0.67449, + # so std ≈ MAD / 0.67449 ≈ 1.4826 * MAD. + filling_outlier_mad_scale_factor: 1.4826 + + # Kernel size in pixels + kernel_range: 100 + kernel_azimuth: 100 + # Sigma for gaussian filtering in pixel + sigma_range: 33 + sigma_azimuth: 33 + # Number of iteration for low pass filtering + filter_iterations: 1 + # Unwrapping correction bool. If False, only filtering + # is applied + unwrap_correction: True + # Before running dispersive_filter, remove every cluster + # whose size is strictly smaller than `min_cluster_pixels`. + # When `min_cluster_pixels` is 0, this pre-removal step is + # skipped entirely. + min_cluster_pixels: 2 + + bridge_algorithm_enabled: True + # Maximum radius used when bridging disconnected + # regions. + bridge_radius: 500 + # Size of the structuring element used for erosion + # during bridge labeling. + bridge_erosion_size: 2 + # Minimum number of pixel of connected components + # to retain during processing. + bridge_minimum_samples: 14 + # Ramp type to correct before bridge algoritm + bridge_ramp_type: + # During deramping, uniform sampling is applied + # when number of samples is larger than maximum pixel + bridge_ramp_maximum_pixel: 1e6 + + troposphere_delay: + # Boolean flag to activate (True) troposphere delay computation + enabled: False + # pyAPS package is used for default + package: pyaps + # Weather model type + weather_model_type: ERA5 + # Conventional line of sight delay + delay_direction: line_of_sight_mapping + # Delay product, where comb product is the summation of wet and hydrostatic + enable_wet_product: True + enable_hydrostatic_product: True + enable_comb_product: False + + geo2rdr: + # No topo_path provided. Default to scratch_path. + threshold: 1.0e-8 + maxiter: 25 + lines_per_block: 1000 + + rdr2geo: + # No topo_path provided. Default to scratch_path. + threshold: 1.0e-7 + numiter: 25 + extraiter: 10 + lines_per_block: 1000 + + # default to save only layers needed for geo2rdr + InSAR + # x, y, and height required for geo2rdr and InSAR + # layover_shadow required for InSAR + write_x: True + write_y: True + write_z: True + write_incidence: False + write_heading: False + write_local_incidence: False + write_local_psi: False + write_simulated_amplitude: False + write_layover_shadow: False + + coarse_resample: + # Path to coarse offsets directory (geo2rdr offsets) + offsets_dir: + # Number of lines and columns per tile to resample in batch. + # If columns_per_tile is 0, the tiles will span across all columns of + # the raster. + lines_per_tile: 1000 + columns_per_tile: 1000 + + dense_offsets: + # Flag to enable/disable dense offsets estimation + enabled: True + # Lines per block to process in batch + lines_per_block: 1000 + # Path to HDF5 file or directory with geometry-coregistered SLCs. + # Not required as we allow the use of intermediate outputs from + # the previous InSAR module, which is not user-specified + coregistered_slc_path: + # Number of columns for chip/template window in reference image + window_range: 64 + # Number of lines for chip/template window in reference image + window_azimuth: 64 + # Number of columns for search chip/template in secondary image + half_search_range: 20 + # Number of lines for search chip/template in secondary image + half_search_azimuth: 20 + # Number of columns to skip in reference image for next offset estimate + skip_range: 32 + # Number of lines to skip in reference image for next offset estimate + skip_azimuth: 32 + # Margin around image edges to avoid for offset computation + margin: 0 + # Number of columns of output offset (covariance, snr) file + offset_width: + # Number of lines of output offset (covariance, snr) file + offset_length: + # Start pixel of the reference image along slant range + start_pixel_range: + # Start pixel of the reference image along azimuth + start_pixel_azimuth: + # Domain (frequency or spatial) where to compute cross-correlation + cross_correlation_domain: 'frequency' + # Anti-aliasing oversampling factor to apply to reference/secondary SLC + # prior to cross-correlation computation + slc_oversampling_factor: 2 + # Deramp data prior to FFT: 'magnitude', 'complex' (linear phase ramp), or 'complex_no_deramp' + deramping_method: 'complex' + # Deramp data axis prior to FFT: range, azimuth, or both + deramping_axis: 'azimuth' + # Flag to use constant range/azimuth offsets in dense offsets estimation + use_gross_offsets: True + # Constant offset along slant range to guide dense offset estimation + gross_offset_range: 0 + # Constant gross offset along az to guide dense offset estimation + gross_offset_azimuth: 0 + # File path to pixel by pixel gross offset (e.g. from reference velocity map) + gross_offset_filepath: + # Flag to add gross offset to the output of dense offset file + merge_gross_offset: + # Square zoom window size for cross-correlation statistics (SNR, covariance) + correlation_statistics_zoom: 21 + # Square zoom window size for sub-pixel offset refinement + correlation_surface_zoom: 8 + # Oversampling factor of the correlation surface for sub-pixel offset refinement + correlation_surface_oversampling_factor: 64 + # Correlation surface oversampling algorithm + correlation_surface_oversampling_method: 'fft' + # Number of cuda streams + cuda_streams: + # Number of offset estimates to process in batch along slant range + windows_batch_range: 10 + # Number of offset estimates to process in batch along azimuth + windows_batch_azimuth: 1 + + offsets_product: + enabled: False + # Path to HDF5 file or directory with geometry-coregistered SLCs. + # Not required as we allow the use of intermediate outputs from + # the previous InSAR module, which is not user-specified + coregistered_slc_path: + # Lines per block + lines_per_block: 1000 + # Margin round ref RSLC edges to exclude during cross-correlation + margin: 0 + # Number of lines of offset raster + offset_length: + # Number of samples of offset raster + offset_width: + # Skip window along range + skip_range: 15 + # Skip window along azimuth + skip_azimuth: 15 + # Start pixel in reference RSLC along range + start_pixel_range: + # Start pixel in reference RSLC along azimuth + start_pixel_azimuth: + # Cross-correlation domain + cross_correlation_domain: 'frequency' + # Deramp data prior to FFT: 'magnitude', 'complex' (linear phase ramp), or 'complex_no_deramp' + deramping_method: 'complex' + # Deramp data axis prior to FFT: range, azimuth, or both + deramping_axis: 'azimuth' + # Anti-aliasing oversampling factor to apply to reference/secondary SLC + # prior to cross-correlation computation + slc_oversampling_factor: 2 + # File path to pixel by pixel gross offset (e.g. from reference velocity map) + gross_offset_filepath: + # Constant offset between ref and secondary RSLCs in slant range + gross_offset_range: 0 + # Constant offset between ref and secondary RSLCs in azimuth + gross_offset_azimuth: 0 + # Merge gross offset file + merge_gross_offset: False + # Correlation surface zoom for statistics + correlation_statistics_zoom: 21 + # Correlation surface zoom + correlation_surface_zoom: 8 + # Correlation surface oversampling factor + correlation_surface_oversampling_factor: 64 + # Correlation surface oversampling method + correlation_surface_oversampling_method: 'sinc' + # Number of cuda streams + cuda_streams: 2 + # Number of windows in range to process in batch + windows_batch_range: 10 + # Number of windows in azimuth to process in batch + windows_batch_azimuth: 5 + # First layer of offsets (same as dense_offsets params) + layer1: + window_range: 64 + window_azimuth: 64 + half_search_range: 16 + half_search_azimuth: 16 + + rubbersheet: + # Path to dense offsets outputs (offsets, snr, covariance). + # Not required as InSAR workflow allows using intermediate outputs + # from previous steps (not user-specified) + dense_offsets_path: + # Path to offsets product outputs (offsets, snr, covariance). + # Not required as InSAR workflow allows using intermediate outputs + # from previous steps (not user-specified) + offsets_product_path: + # Path to geo2rdr offsets + geo2rdr_offsets_path: + # Flag to enable/disable rubbersheet step + enabled: True + # Metric used to identify offsets outliers + culling_metric: median_filter + # Size of median and moving average filter in range. Median filter is used + # for offset outlier identification. Moving average is used to smooth the + # culled offset field prior to interferogram formation + median_filter_size_range: 9 + # Size of median and moving average filter in azimuth (see above) + median_filter_size_azimuth: 9 + # Threshold for outlier identification. Depends on "culling_metric". + # SNR: offsets with SNR below threshold are identified as outliers + # Covariance: offsets with range/azimuth covariance above threshold are outliers + # Median filter: offsets with median absolute deviation above thresholds are outliers + threshold: 0.75 + # Flag to enable mask refinement The intent of this refinement is to remove + # noisy spots on 1-2 pixels after the first outlier identification. + # The refinement is performed by thresholding the median absolute deviation of + # the offsets where outlier locations have been filled with zeros + mask_refine_enabled: True + # Minimum number of valid (non-NaN) neighbors required for a pixel to remain + # unchanged during the rubbersheet refinement step + mask_refine_min_neighbors: 5 + # Size of median filter in range and azimuth for outlier mask refinement + mask_refine_filter_size: 5 + # Method to fill data holes left by outliers culling: + # fill_smoothed: iteratively fill holes with mean value in a neighborhood. + # nearest_neighbor: Fill holes with nearest neighbor interpolation + outlier_filling_method: fill_smoothed + fill_smoothed: + # Number of fill smoothed iterations + iterations: 1 + # Size of the neighborhood to look for smoothed values + kernel_size: 3 + # Interpolation method. Interpolation is used to fill residual outlier + # holes if present + interpolation_method: linear + + # Additional azimuth offset filtering applied after outlier detection + # to suppress offset jumps caused by ionospheric disturbances, + # decorrelation, or processing artifacts + azimuth_offset_filter: + # Spatial filter type: 'none' (disabled), 'median', or 'mean' + offsets_filter: none + # Kernel size in pixels for spatial filtering + kernel_size: 31 + + # Flag to enable the use of subswath mask for the outlier detection + subswath_mask_apply_enabled: True + + # IDW interpolation method parameters + idw_interpolation: + # Power parameter controlling the distance weighting + power: 2 + # Number of nearest neighbors used for interpolation + number: 100 + # Maximum search radius in pixels + radius: 200 + + # It is good practice to filter the offsets prior to interferogram + # formation to reduce noise. We expose: no_filter: do not filter the offsets/ + # degrade offsets resolution. median_filter, boxcar_filter (moving average), + # gaussian_filter + offsets_filter: boxcar + boxcar: + # Filter size (median/boxcar) in slant range + filter_size_range: 5 + filter_size_azimuth: 5 + # Polyfit the dense offsets for the fine resampling + polyfitting: + # Flag to enable/disable polyfitting + enabled: false + # The degrees for the polyfitting + degree: 2 + # Critical value for the w-test to remove the outliers + critical_value: 0.1 + # Number of samples along the range in pixels for polyfitting + samples_along_range: 200 + # Number of samples along the azimuth in pixels for polyfitting + samples_along_azimuth: 200 + + fine_resample: + # Flag to enable/disable fine resampling + enabled: True + # Path to fine offsets (rubbersheet dense offsets) + offsets_dir: + # Number of lines and columns per tile to resample in batch. + # If columns_per_tile is 0, the tiles will span across all columns of + # the raster. + lines_per_tile: 100 + columns_per_tile: 0 + + crossmul: + # Number of range looks to generate wrapped interferogram in RIFG + range_looks: 1 + # Number of azimuth looks to generate wrapped interferogram in RIFG + azimuth_looks: 1 + flatten: True + oversample: 2 + lines_per_block: 1024 + # If True, activates common band range filtering + common_band_range_filter: False + # If True, activates common band azimuth filtering + common_band_azimuth_filter: False + + filter_interferogram: + interferogram_path: + mask: + general: + lines_per_block: 100 + filter_type: no_filter + boxcar: + filter_size_range: 9 + filter_size_azimuth: 9 + + baseline: + # 3D_full: compute baselines for all heights in metadata cube + # top_bottom: compute baselines for top and bottom heights + mode: 'top_bottom' + + phase_unwrap: + # Path to HDF5 file or directory containing the input interferogram + # and coherence (normalized magnitude of complex correlation) rasters + # Not required as we allow the use of intermediate outputs from + # the previous InSAR module, which is not user-specified. + # Required for running stand-alone phase unwrapping + crossmul_path: + # Number of looks in slant range and azimuth directions to generate the + # wrapped interferogram that will be unwrapped in RUNW. + # If range_looks = 1 and azimuth_looks = 1 no further multilooking occurs + # In this case the wrapped interferogram in RIFG will be unwrapped. + range_looks: 1 + azimuth_looks: 1 + # Mask and fill wrapped InSAR phase prior to unwrapping + preprocess_wrapped_phase: + # Enable/disable preprocessing of wrapped igram + enabled: False + # Mask for outlier identification. User-provided or data-driven + mask: + # Path to user-provided mask + mask_path: + # Type of data-driven mask to compute. Data-driven masking + # is ignored if mask_path is provided. mask_type can be: + # 1) median_filter: compute mask with median absolute deviation thresholding + # 2) coherence: compute mask with normalized InSAR coherence thresholding + # 3) water: extract water areas from water mask file in dynamic ancillary data + mask_type: coherence + # buffer from the land. Unit is [km]. Buffer should be integer. + ocean_water_buffer: 1 + inland_water_buffer: 1 + # Threshold for invalid pixels. If mask_type is coherence (median_filter) + # values below (above) threshold are classified as outliers + outlier_threshold: 0.5 + # Size of median filter kernel (same shape in range/azimuth). Only used if + # "mask_type" is "median_filter" + median_filter_size: 25 + # Enable/disable filling invalid values. If disabled, invalid values + # are replaced with zero. + filling_enabled: False + # Algorithm to fill invalid wrapped igram pixels + filling_method: distance_interpolator + # Distance interpolator options based on Chen et al. 2015 + # https://doi.org/10.1002/2015GL065031 + distance_interpolator: + # Distance/radius of the sampling window (pixels) + distance: 2 + # Unwrapping algorithm to use + algorithm: snaphu + # Default parameters for ICU + icu: + # Seed point to initialize ICU algorithm + seed: 0 + # Number of lines in buffer to process simultaneously + buffer_lines: 3700 + # Number of overlapping lines for each data block + overlap_lines: 200 + # Flag to activate neutrons based on phase gradient metric + use_phase_gradient_neutron: False + # Flag to activate neutron based on reference RSLC intensity + use_intensity_neutron: False + # Square window size for phase gradient computation + phase_gradient_window_size: 5 + # Threshold value for phase gradient neutrons (radian) + neutron_phase_gradient_threshold: 3.0 + # Threshold value for intensity neutron (sigma above mean) + neutron_intensity_threshold: 8.0 + # Threshold on the maximum level of intensity for intensity neutron + max_intensity_correlation_threshold: 0.8 + # Maximum number of tree realizations + trees_number: 7 + # Maximum extent for each tree branch + max_branch_length: 64 + # Ratio of pixel spacing along slant range and azimuth + pixel_spacing_ratio: 1.0 + # Threshold for the initial correlation value + initial_correlation_threshold: 0.1 + # Threshold for the maximum correlation value + max_correlation_threshold: 0.9 + # Increments to correlation threshold + correlation_threshold_increments: 0.1 + # Minimum size of an area to unwrap + min_tile_area: 0.003125 + # Number of lines used in the bootstrap phase + bootstrap_lines: 16 + # Minimum bootstrap overlapping area + min_overlap_area: 16 + # Bootstrap phase variance threshold (radian) + phase_variance_threshold: 8 + phass: + # Increments to correlation threshold + correlation_threshold_increments: 0.2 + # Normalized interferometric correlation value of the pixel to unwrap + # Pixels with a good_correlation lower than the assigned value are + # masking out during the phass unwrapping + good_correlation: 0.7 + # Minimum size (in pixels) of a region to unwrap + min_unwrap_area: 200 + # Filepath to power raster (reference RSLC) + power: + snaphu: + # Effective number of looks used to form the sample coherence. + # It provides an estimate of the number of statistically independent + # samples averaged in the multilooked data, taking into account spatial + # correlation due to filtering/oversampling. If not specified, nlooks (n_e) is + # estimated as: + # $$ + # n_e = k_r k_a \frac{d_r d_a}{\rho_r \rho_a} + # $$ + # where $k_r$ and $k_a$ are the number of looks in range and azimuth, + # $d_r$ and $d_a$ are the (single-look) sample spacing in range and + # azimuth, and $\rho_r$ and $\rho_a$ are the range and azimuth resolution. + nlooks: + # Statistical cost mode. One of {smooth, defo} + cost_mode: smooth + # File path to a binary mask of valid pixels. Zeros in this raster indicate + # interferogram pixels that should be masked out. If provided, it must have + # the same dimensions as the input interferogram and boolean or 8-bit integer datatype + mask: + # Algorithm used for initialize unwrapped phase gradients. + # Supported algorithms: MST (Minimum Spanning Tree); MCF + # (Minimum Cost Flow) + initialization_method: mcf + # Minimum size of a single connected component, as a fraction of the total number + # of pixels in the tile. + min_conncomp_frac: 0.01 + # The dimensions, in pixels, of the sliding window used for averaging wrapped + # phase gradients to get the mean non-layover slope, in directions parallel and + # perpendicular to the examined phase difference. This parameter corresponds to + # SNAPHU's `KPARDPSI` and `KPERPDPSI` configuration options. + phase_grad_window: [7, 7] + # Number of tiles along the row/column directions. If `ntiles` is (1, 1), then the + # interferogram will be unwrapped as a single tile. Increasing the number of tiles + # may improve runtime and reduce peak memory utilization, but may also introduce + # tile boundary artifacts in the unwrapped result. + ntiles: [1, 1] + # Overlap, in pixels, between neighboring tiles. Increasing overlap may help to + # avoid phase discontinuities between tiles. If `tile_overlap` is a scalar + # integer, the number of overlapping rows and columns will be the same. + tile_overlap: [4, 4] + # Maximum number of child processes to spawn for parallel tile unwrapping. If + # `nproc` is less than 1, use all available processors. Defaults to 1. + nproc: 20 + # Cost threshold to use for determining boundaries of reliable regions + # (dimensionless; scaled according to other cost constants). Larger cost threshold + # implies smaller regions -- safer, but more expensive computationally + tile_cost_thresh: 500 + # Minimum size, in pixels, of a reliable region in tile mode. + min_region_size: 300 + # If True, after unwrapping with multiple tiles, an additional post-processing + # unwrapping step is performed to re-optimize the unwrapped phase using a single + # tile. This option is disregarded when `ntiles` is (1, 1). It supersedes the + # `regrow_conncomps` option -- if both are enabled, only the single-tile + # re-optimization step will be performed in order to avoid redundant computation. + single_tile_reoptimize: True + # If True, the connected component labels will be re-computed using a single tile + # after first unwrapping with multiple tiles. This option is disregarded when + # `ntiles` is (1, 1). It has no effect if `single_tile_reoptimize` was also + # enabled (since regrowing connected components would be redundant in that case). + regrow_conncomps: True + bridge: + enabled: True + # Maximum radius used when bridging disconnected regions. + bridge_radius: 500 + # Size of the structuring element used for erosion + # during bridge labeling. + bridge_erosion_size: 2 + # Minimum number of pixel of connected components + # to retain during processing. + bridge_minimum_samples: 14 + # Ramp type to correct before bridge algoritm + bridge_ramp_type: + # During deramping, uniform sampling is applied + # when number of samples is larger than maximum pixel + bridge_ramp_maximum_pixel: 1e6 + + correction_luts: + # Boolean flag to activate/deactivate model-based solid earth tide corrections + solid_earth_tides_enabled: True + + # Set of options to enable compression, choose compression level, + # chunk size and shuffle filter + output: + # Enable/Disable the paging + page_enabled: True + # Enable/Disable gzip compression + compression_enabled: True + # Compression level: 1: least compression; 9: best compression + compression_level: 1 + # Chunk size of InSAR raster + chunk_size: [512, 512] + # Enable/disable shuffle filter + shuffle: True + + qa: + rifg: + + workflows: + + # Flag to run `validate` workflow to validate the + # input file against its product spec. + # Default: True + validate: true + + # Flag to run `qa_reports` workflow to generate a + # PDF report, geolocated browse image, compute statistics + # on the input file, etc.. + # Default: True + qa_reports: true + + software_config: + + # True to cache selected dataset(s) into intermediate + # memory-mapped flat file(s), which speeds up repeat access. + # False to always read data directly from the input file. + # Generally, enabling caching should reduce runtime. + # Default: True + use_cache: true + + # True to delete the nested QA scratch directory in + # `scratch_path` and its contents when QA SAS is finished. + # Default: True + delete_qa_scratch_files: true + + validation: + + # True to raise an exception if one or more metadata LUTs + # contain all non-finite (e.g. Nan, +/- Inf) values, or if one or more + # z-dimension height layers in a 3D LUT ("metadata cube") has + # all non-finite values. False to quiet the exception (although + # it will still be logged). + # Default: True + metadata_luts_fail_if_all_nan: true + + qa_reports: + + # browse: + + # # The maximum number of pixels allowed for the longest + # # side of the final 2D browse image. + # # Default: 2048 + # longest_side_max: 2048 + + # # True to generate a version of the native browse PNG+KML that is + # # projected to EPSG 4326 (lon/lat) coordinate system for accurate + # # geolocation in GIS software. + # # False to only generate native browse outputs. + # # The EPSG 4326 PNG+KML version will be in addition to the native + # # browse PNG+KML outputs (where the PNG reflects the input product's + # # native coordinate system). The EPSG 4326 outputs will be suffixed with + # # '_LATLON' (e.g., 'BROWSE_LATLON.png' and 'BROWSE_LATLON.kml'). + # # Default: False + # output_browse_latlon: false + + # # Resampling method for ISCE3 geocoding. Options: 'sinc', + # # 'bilinear', 'bicubic', 'nearest', 'biquintic'. + # # Ignored if `output_browse_latlon` is False. + # # When output browse image data contains phase discontinuities + # # (e.g. RIFG's wrapped interferogram), strongly recommend 'nearest'; + # # other methods can create artifacts when interpolating across the jump. + # # Default: nearest + # resample: nearest + + wrapped_igram: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 95.0 + nan_threshold: 95.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 95.0 + inf_threshold: 95.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 95.0 + fill_threshold: 95.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 95.0 + total_invalid_threshold: 95.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + coh_mag: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 95.0 + nan_threshold: 95.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 95.0 + inf_threshold: 95.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 95.0 + fill_threshold: 95.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: 95.0 + near_zero_threshold: 95.0 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 95.0 + total_invalid_threshold: 95.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: True + zero_is_invalid: true + + az_and_range_offsets: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 95.0 + nan_threshold: 95.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 95.0 + inf_threshold: 95.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 95.0 + fill_threshold: 95.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 95.0 + total_invalid_threshold: 95.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + correlation_surface_peak: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 95.0 + nan_threshold: 95.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 95.0 + inf_threshold: 95.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 95.0 + fill_threshold: 95.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: 95.0 + near_zero_threshold: 95.0 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 95.0 + total_invalid_threshold: 95.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: True + zero_is_invalid: true + + runw: + + workflows: + + # Flag to run `validate` workflow to validate the + # input file against its product spec. + # Default: True + validate: true + + # Flag to run `qa_reports` workflow to generate a + # PDF report, geolocated browse image, compute statistics + # on the input file, etc.. + # Default: True + qa_reports: true + + software_config: + + # True to cache selected dataset(s) into intermediate + # memory-mapped flat file(s), which speeds up repeat access. + # False to always read data directly from the input file. + # Generally, enabling caching should reduce runtime. + # Default: True + use_cache: true + + # True to delete the nested QA scratch directory in + # `scratch_path` and its contents when QA SAS is finished. + # Default: True + delete_qa_scratch_files: true + + validation: + + # True to raise an exception if one or more metadata LUTs + # contain all non-finite (e.g. Nan, +/- Inf) values, or if one or more + # z-dimension height layers in a 3D LUT ("metadata cube") has + # all non-finite values. False to quiet the exception (although + # it will still be logged). + # Default: True + metadata_luts_fail_if_all_nan: true + + qa_reports: + + # browse: + + # # The maximum number of pixels allowed for the longest + # # side of the final 2D browse image. + # # Default: 2048 + # longest_side_max: 2048 + + # # The multiple of pi to rewrap the unwrapped phase image + # # when generating the browse PNG. If None, no rewrapping will occur. + # # Ex: If 3 is provided, the image is rewrapped to the interval [0, 3pi). + # # Default: 7 + # rewrap: 7 + + # # True to generate a version of the native browse PNG+KML that is + # # projected to EPSG 4326 (lon/lat) coordinate system for accurate + # # geolocation in GIS software. + # # False to only generate native browse outputs. + # # The EPSG 4326 PNG+KML version will be in addition to the native + # # browse PNG+KML outputs (where the PNG reflects the input product's + # # native coordinate system). The EPSG 4326 outputs will be suffixed with + # # '_LATLON' (e.g., 'BROWSE_LATLON.png' and 'BROWSE_LATLON.kml'). + # # Default: False + # output_browse_latlon: false + + # # Resampling method for ISCE3 geocoding. Options: 'sinc', + # # 'bilinear', 'bicubic', 'nearest', 'biquintic'. + # # Ignored if `output_browse_latlon` is False. + # # When output browse image data contains phase discontinuities + # # (e.g. RIFG's wrapped interferogram), strongly recommend 'nearest'; + # # other methods can create artifacts when interpolating across the jump. + # # Default: bilinear + # resample: bilinear + + phase_img: + + # The multiple of pi to rewrap the unwrapped phase image in the report + # PDF. If None, no rewrapping will occur. + # Ex: If 3 is provided, the image is rewrapped to the interval [0, 3pi). + # Default: 7 + rewrap: 7 + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 95.0 + nan_threshold: 95.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 95.0 + inf_threshold: 95.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 95.0 + fill_threshold: 95.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: 95.0 + near_zero_threshold: 95.0 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 95.0 + total_invalid_threshold: 95.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: True + zero_is_invalid: true + + coh_mag: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 95.0 + nan_threshold: 95.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 95.0 + inf_threshold: 95.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 95.0 + fill_threshold: 95.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: 95.0 + near_zero_threshold: 95.0 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 95.0 + total_invalid_threshold: 95.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: True + zero_is_invalid: true + + connected_components: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 0 + nan_threshold: 0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 0 + inf_threshold: 0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 95.0 + fill_threshold: 95.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: 99.9 + near_zero_threshold: 99.9 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 99.9 + total_invalid_threshold: 99.9 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: True + zero_is_invalid: true + + # Maximum number of valid connected components allowed. + # If the number of valid connected components (not including + # zero nor the fill value) is greater than this value, + # it will be logged and an exception will be raised. + # If None, this error check will be skipped. + # Default: None + max_num_cc: + + ionosphere_phase_screen: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: -1 + nan_threshold: -1 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 95.0 + inf_threshold: 95.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: -1 + fill_threshold: -1 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: -1 + total_invalid_threshold: -1 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + ionosphere_phase_screen_uncertainty: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 95.0 + nan_threshold: 95.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 95.0 + inf_threshold: 95.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 95.0 + fill_threshold: 95.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 95.0 + total_invalid_threshold: 95.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + az_and_range_offsets: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 95.0 + nan_threshold: 95.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 95.0 + inf_threshold: 95.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 95.0 + fill_threshold: 95.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 95.0 + total_invalid_threshold: 95.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + correlation_surface_peak: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 95.0 + nan_threshold: 95.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 95.0 + inf_threshold: 95.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 95.0 + fill_threshold: 95.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: 95.0 + near_zero_threshold: 95.0 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 95.0 + total_invalid_threshold: 95.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: True + zero_is_invalid: true + + gunw: + + workflows: + + # Flag to run `validate` workflow to validate the + # input file against its product spec. + # Default: True + validate: true + + # Flag to run `qa_reports` workflow to generate a + # PDF report, geolocated browse image, compute statistics + # on the input file, etc.. + # Default: True + qa_reports: true + + software_config: + + # True to cache selected dataset(s) into intermediate + # memory-mapped flat file(s), which speeds up repeat access. + # False to always read data directly from the input file. + # Generally, enabling caching should reduce runtime. + # Default: True + use_cache: true + + # True to delete the nested QA scratch directory in + # `scratch_path` and its contents when QA SAS is finished. + # Default: True + delete_qa_scratch_files: true + + validation: + + # True to raise an exception if one or more metadata LUTs + # contain all non-finite (e.g. Nan, +/- Inf) values, or if one or more + # z-dimension height layers in a 3D LUT ("metadata cube") has + # all non-finite values. False to quiet the exception (although + # it will still be logged). + # Default: True + metadata_luts_fail_if_all_nan: true + + qa_reports: + + # browse: + + # # The maximum number of pixels allowed for the longest + # # side of the final 2D browse image. + # # Default: 2048 + # longest_side_max: 2048 + + # # The multiple of pi to rewrap the unwrapped phase image + # # when generating the browse PNG. If None, no rewrapping will occur. + # # Ex: If 3 is provided, the image is rewrapped to the interval [0, 3pi). + # # Default: 7 + # rewrap: 7 + + # # True to generate a version of the native browse PNG+KML that is + # # projected to EPSG 4326 (lon/lat) coordinate system for accurate + # # geolocation in GIS software. + # # False to only generate native browse outputs. + # # The EPSG 4326 PNG+KML version will be in addition to the native + # # browse PNG+KML outputs (where the PNG reflects the input product's + # # native coordinate system). The EPSG 4326 outputs will be suffixed with + # # '_LATLON' (e.g., 'BROWSE_LATLON.png' and 'BROWSE_LATLON.kml'). + # # Default: False + # output_browse_latlon: false + + # # Resampling algorithm for GDAL reprojection. Options: 'near', + # # 'bilinear', 'cubic', 'cubicspline', 'lanczos', 'average', 'mode'. + # # Ignored if `output_browse_latlon` is False. + # # Default: average + # resample: average + + phase_img: + + # The multiple of pi to rewrap the unwrapped phase image in the report + # PDF. If None, no rewrapping will occur. + # Ex: If 3 is provided, the image is rewrapped to the interval [0, 3pi). + # Default: 7 + rewrap: 7 + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 99.0 + nan_threshold: 99.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 99.0 + inf_threshold: 99.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 99.0 + fill_threshold: 99.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: 99.0 + near_zero_threshold: 99.0 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 99.0 + total_invalid_threshold: 99.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: True + zero_is_invalid: true + + wrapped_igram: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 99.0 + nan_threshold: 99.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 99.0 + inf_threshold: 99.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 99.0 + fill_threshold: 99.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 99.0 + total_invalid_threshold: 99.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + coh_mag: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 99.0 + nan_threshold: 99.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 99.0 + inf_threshold: 99.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 99.0 + fill_threshold: 99.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: 99.0 + near_zero_threshold: 99.0 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 99.0 + total_invalid_threshold: 99.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: True + zero_is_invalid: true + + connected_components: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 0 + nan_threshold: 0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 0 + inf_threshold: 0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 99.0 + fill_threshold: 99.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: 99.9 + near_zero_threshold: 99.9 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 99.9 + total_invalid_threshold: 99.9 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: True + zero_is_invalid: true + + # Maximum number of valid connected components allowed. + # If the number of valid connected components (not including + # zero nor the fill value) is greater than this value, + # it will be logged and an exception will be raised. + # If None, this error check will be skipped. + # Default: None + max_num_cc: + + ionosphere_phase_screen: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: -1 + nan_threshold: -1 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 99.0 + inf_threshold: 99.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: -1 + fill_threshold: -1 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: -1 + total_invalid_threshold: -1 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + ionosphere_phase_screen_uncertainty: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 99.0 + nan_threshold: 99.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 99.0 + inf_threshold: 99.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 99.0 + fill_threshold: 99.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 99.0 + total_invalid_threshold: 99.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + az_and_range_offsets: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 99.0 + nan_threshold: 99.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 99.0 + inf_threshold: 99.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 99.0 + fill_threshold: 99.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 99.0 + total_invalid_threshold: 99.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + correlation_surface_peak: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 99.0 + nan_threshold: 99.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 99.0 + inf_threshold: 99.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 99.0 + fill_threshold: 99.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: 99.0 + near_zero_threshold: 99.0 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 99.0 + total_invalid_threshold: 99.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: True + zero_is_invalid: true + + roff: + + workflows: + + # Flag to run `validate` workflow to validate the + # input file against its product spec. + # Default: True + validate: true + + # Flag to run `qa_reports` workflow to generate a + # PDF report, geolocated browse image, compute statistics + # on the input file, etc.. + # Default: True + qa_reports: true + + software_config: + + # True to cache selected dataset(s) into intermediate + # memory-mapped flat file(s), which speeds up repeat access. + # False to always read data directly from the input file. + # Generally, enabling caching should reduce runtime. + # Default: True + use_cache: true + + # True to delete the nested QA scratch directory in + # `scratch_path` and its contents when QA SAS is finished. + # Default: True + delete_qa_scratch_files: true + + validation: + + # True to raise an exception if one or more metadata LUTs + # contain all non-finite (e.g. Nan, +/- Inf) values, or if one or more + # z-dimension height layers in a 3D LUT ("metadata cube") has + # all non-finite values. False to quiet the exception (although + # it will still be logged). + # Default: True + metadata_luts_fail_if_all_nan: true + + qa_reports: + + # browse: + + # # The maximum number of pixels allowed for the longest + # # side of the final 2D browse image PNG. + # # Default: 2048 + # longest_side_max: 2048 + + # # Stride along each axis of the Frequency A + # # image arrays for decimating the raster for the browse PNG. + # # This takes precedence over `longest_side_max`. + # # Format: [, ] + # # Example: [6,7] + # # If None, QA-SAS will compute the decimation strides based on `longeset_side_max` + # # and decimating the raster such that the pixels are roughly square. + # # Default: None + # browse_decimation_freqa: + + # # Stride along each axis of the Frequency B + # # image arrays for decimating the raster for the browse PNG. + # # This takes precedence over `longest_side_max`. + # # Format: [, ] + # # Example: [6,7] + # # If None, QA-SAS will compute the decimation strides based on `longeset_side_max` + # # and decimating the raster such that the pixels are roughly square. + # # Default: None + # browse_decimation_freqb: + + # # True to generate a version of the native browse PNG+KML that is + # # projected to EPSG 4326 (lon/lat) coordinate system for accurate + # # geolocation in GIS software. + # # False to only generate native browse outputs. + # # The EPSG 4326 PNG+KML version will be in addition to the native + # # browse PNG+KML outputs (where the PNG reflects the input product's + # # native coordinate system). The EPSG 4326 outputs will be suffixed with + # # '_LATLON' (e.g., 'BROWSE_LATLON.png' and 'BROWSE_LATLON.kml'). + # # Default: False + # output_browse_latlon: false + + # # Resampling method for ISCE3 geocoding. Options: 'sinc', + # # 'bilinear', 'bicubic', 'nearest', 'biquintic'. + # # Ignored if `output_browse_latlon` is False. + # # When output browse image data contains phase discontinuities + # # (e.g. RIFG's wrapped interferogram), strongly recommend 'nearest'; + # # other methods can create artifacts when interpolating across the jump. + # # Default: bilinear + # resample: bilinear + + az_and_range_offsets: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 95.0 + nan_threshold: 95.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 95.0 + inf_threshold: 95.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 95.0 + fill_threshold: 95.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 95.0 + total_invalid_threshold: 95.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + quiver_plots: + + # The vmin and vmax values to generate the quiver plot. + # The magnitude of the offsets is clipped to this interval, + # which (in turn) is used for the interval of the colorbar. + # If None, the interval is computed using the min and max + # magnitudes of along track and slant range offset. + # Default: None + colorbar_min_max: + + # Number of arrows (vectors) to plot per the + # longest edge of the raster. + # Default: 20 + arrow_density: 20 + + # Scales the length of the arrow inversely. + # Number of data units per arrow length unit, e.g., m/s per plot width; + # a smaller scale parameter makes the arrow longer. + # If None, a simple autoscaling algorithm is used, based on the average + # vector length and the number of vectors. + # See: The `scaling` parameter for `matplotlib.axes.Axes.quiver()`. + # Default: None + arrow_scaling: + + az_and_rg_variance: + + # The vmin and vmax values to generate the plots for the + # az and slant range offset variance layers for ROFF and GOFF. + # The square root of these layers (i.e. the standard deviation + # of the offsets) is computed, clipped to this interval, and + # then plotted using this interval for the colorbar. + # If None, the interval is computed by: + # If the variance layers' units are 'meters^2': [0.0, 10.0] + # Else-if the variance layers' units are 'pixels^2': [0.0, 0.1] + # Otherwise: [0.0, max(sqrt(), sqrt())] + # Default: None + colorbar_min_max: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 95.0 + nan_threshold: 95.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 95.0 + inf_threshold: 95.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 95.0 + fill_threshold: 95.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 95.0 + total_invalid_threshold: 95.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + cross_offset_variance: + + # The vmin and vmax values to generate the plots + # for the cross offset variance layer for ROFF and GOFF. + # If None, then the colorbar range will be computed based + # on `percentile_for_clipping`. + # Default: None + colorbar_min_max: + + # Percentile range that the cross offset variance raster + # will be clipped to, which determines the colormap interval. + # Must be in range [0.0, 100.0]. + # Superseded by `cbar_min_max` parameter. + # Default: (1.0, 99.0) + percentile_for_clipping: [1.0, 99.0] + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 95.0 + nan_threshold: 95.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 95.0 + inf_threshold: 95.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 95.0 + fill_threshold: 95.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 95.0 + total_invalid_threshold: 95.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + correlation_surface_peak: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 95.0 + nan_threshold: 95.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 95.0 + inf_threshold: 95.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 95.0 + fill_threshold: 95.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: 95.0 + near_zero_threshold: 95.0 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 95.0 + total_invalid_threshold: 95.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: True + zero_is_invalid: true + + goff: + + workflows: + + # Flag to run `validate` workflow to validate the + # input file against its product spec. + # Default: True + validate: true + + # Flag to run `qa_reports` workflow to generate a + # PDF report, geolocated browse image, compute statistics + # on the input file, etc.. + # Default: True + qa_reports: true + + software_config: + + # True to cache selected dataset(s) into intermediate + # memory-mapped flat file(s), which speeds up repeat access. + # False to always read data directly from the input file. + # Generally, enabling caching should reduce runtime. + # Default: True + use_cache: true + + # True to delete the nested QA scratch directory in + # `scratch_path` and its contents when QA SAS is finished. + # Default: True + delete_qa_scratch_files: true + + validation: + + # True to raise an exception if one or more metadata LUTs + # contain all non-finite (e.g. Nan, +/- Inf) values, or if one or more + # z-dimension height layers in a 3D LUT ("metadata cube") has + # all non-finite values. False to quiet the exception (although + # it will still be logged). + # Default: True + metadata_luts_fail_if_all_nan: true + + qa_reports: + + # browse: + + # # The maximum number of pixels allowed for the longest + # # side of the final 2D browse image PNG. + # # Default: 2048 + # longest_side_max: 2048 + + # # Stride along each axis of the Frequency A + # # image arrays for decimating the raster for the browse PNG. + # # This takes precedence over `longest_side_max`. + # # Format: [, ] + # # Example: [6,7] + # # If None, QA-SAS will compute the decimation strides based on `longeset_side_max` + # # and decimating the raster such that the pixels are roughly square. + # # Default: None + # browse_decimation_freqa: + + # # Stride along each axis of the Frequency B + # # image arrays for decimating the raster for the browse PNG. + # # This takes precedence over `longest_side_max`. + # # Format: [, ] + # # Example: [6,7] + # # If None, QA-SAS will compute the decimation strides based on `longeset_side_max` + # # and decimating the raster such that the pixels are roughly square. + # # Default: None + # browse_decimation_freqb: + + # # True to generate a version of the native browse PNG+KML that is + # # projected to EPSG 4326 (lon/lat) coordinate system for accurate + # # geolocation in GIS software. + # # False to only generate native browse outputs. + # # The EPSG 4326 PNG+KML version will be in addition to the native + # # browse PNG+KML outputs (where the PNG reflects the input product's + # # native coordinate system). The EPSG 4326 outputs will be suffixed with + # # '_LATLON' (e.g., 'BROWSE_LATLON.png' and 'BROWSE_LATLON.kml'). + # # Default: False + # output_browse_latlon: false + + # # Resampling algorithm for GDAL reprojection. Options: 'near', + # # 'bilinear', 'cubic', 'cubicspline', 'lanczos', 'average', 'mode'. + # # Ignored if `output_browse_latlon` is False. + # # Default: average + # resample: average + + az_and_range_offsets: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 99.0 + nan_threshold: 99.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 99.0 + inf_threshold: 99.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 99.0 + fill_threshold: 99.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 99.0 + total_invalid_threshold: 99.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + quiver_plots: + + # The vmin and vmax values to generate the quiver plot. + # The magnitude of the offsets is clipped to this interval, + # which (in turn) is used for the interval of the colorbar. + # If None, the interval is computed using the min and max + # magnitudes of along track and slant range offset. + # Default: None + colorbar_min_max: + + # Number of arrows (vectors) to plot per the + # longest edge of the raster. + # Default: 20 + arrow_density: 20 + + # Scales the length of the arrow inversely. + # Number of data units per arrow length unit, e.g., m/s per plot width; + # a smaller scale parameter makes the arrow longer. + # If None, a simple autoscaling algorithm is used, based on the average + # vector length and the number of vectors. + # See: The `scaling` parameter for `matplotlib.axes.Axes.quiver()`. + # Default: None + arrow_scaling: + + az_and_rg_variance: + + # The vmin and vmax values to generate the plots for the + # az and slant range offset variance layers for ROFF and GOFF. + # The square root of these layers (i.e. the standard deviation + # of the offsets) is computed, clipped to this interval, and + # then plotted using this interval for the colorbar. + # If None, the interval is computed by: + # If the variance layers' units are 'meters^2': [0.0, 10.0] + # Else-if the variance layers' units are 'pixels^2': [0.0, 0.1] + # Otherwise: [0.0, max(sqrt(), sqrt())] + # Default: None + colorbar_min_max: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 99.0 + nan_threshold: 99.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 99.0 + inf_threshold: 99.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 99.0 + fill_threshold: 99.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 99.0 + total_invalid_threshold: 99.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + cross_offset_variance: + + # The vmin and vmax values to generate the plots + # for the cross offset variance layer for ROFF and GOFF. + # If None, then the colorbar range will be computed based + # on `percentile_for_clipping`. + # Default: None + colorbar_min_max: + + # Percentile range that the cross offset variance raster + # will be clipped to, which determines the colormap interval. + # Must be in range [0.0, 100.0]. + # Superseded by `cbar_min_max` parameter. + # Default: (1.0, 99.0) + percentile_for_clipping: [1.0, 99.0] + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 99.0 + nan_threshold: 99.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 99.0 + inf_threshold: 99.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 99.0 + fill_threshold: 99.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: -1 + near_zero_threshold: -1 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 99.0 + total_invalid_threshold: 99.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: False + zero_is_invalid: false + + correlation_surface_peak: + + # *** Threshold Percentage Parameters *** + # Check for malformed datasets using these threshold percentages. + # An exception is raised if the percent of pixels with the corresponding value + # per the total raster area exceeds the threshold percentage provided. + # Thresholds must be in the interval [0, 100], or -1 to be ignored. + # (-1 logs the percentages as INFO, but will not trigger a QA failure.) + + # Percent of NaN pixels per total raster area allowed + # Default: 99.0 + nan_threshold: 99.0 + + # Percent of +/- Inf pixels per total raster area allowed + # Default: 99.0 + inf_threshold: 99.0 + + # Percent of Fill-valued pixels per total raster area allowed + # Default: 99.0 + fill_threshold: 99.0 + + # Percent of Near-zero pixels per total raster area allowed + # Default: 99.0 + near_zero_threshold: 99.0 + + # Percent of total invalid pixels per total raster area allowed + # 'Total invalid pixels' is the sum of NaN, Inf, fill, & (if `zero_is_invalid` is True) near-zero pixels. + # Default: 99.0 + total_invalid_threshold: 99.0 + + # Absolute tolerance for determining if a raster pixel is 'near zero'. + # Default: 1e-06 + epsilon: 1e-06 + + # True if near-zero pixels should be counted towards the + # total number of invalid pixels. False to exclude them. + # If False, consider setting `near_zero_threshold` to -1. + # Note: fill values are always considered invalid, so if a raster's + # fill value is zero, then zeros will still be included in the total. + # Default: True + zero_is_invalid: true \ No newline at end of file diff --git a/share/stack/invert_offsets.py b/share/stack/invert_offsets.py new file mode 100644 index 0000000000..a5c7bd2413 --- /dev/null +++ b/share/stack/invert_offsets.py @@ -0,0 +1,884 @@ +from pathlib import Path + +import numpy as np +from osgeo import gdal + + +def _solve_offset_block(A, L): + """ + Solve: + + L_ij = u_j - u_i + + for one block of pixels. + + Parameters + ---------- + A : (n_pairs, n_dates - 1) + Network design matrix with reference date removed. + + L : (n_pairs, n_pixels) + Pairwise offsets. + + Returns + ------- + U : (n_dates - 1, n_pixels) + Per-date offsets relative to reference. + + Notes + ----- + Handles NaNs. If all pairs are valid, uses one matrix + multiplication. For pixels containing NaNs, uses the + available pair network if it remains full rank. + """ + + n_unknown = A.shape[1] + n_pix = L.shape[1] + + U = np.full( + (n_unknown, n_pix), + np.nan, + dtype=np.float64, + ) + + # -------------------------------------------------- + # Fast path: all pairs valid + # -------------------------------------------------- + + full_valid = np.all(np.isfinite(L), axis=0) + + if np.any(full_valid): + Ainv = np.linalg.pinv(A) + + U[:, full_valid] = ( + Ainv @ L[:, full_valid] + ) + + # -------------------------------------------------- + # Pixels with missing pair observations + # -------------------------------------------------- + + bad_pixels = np.where(~full_valid)[0] + + # Cache pseudoinverses for repeated validity patterns + pinv_cache = {} + + for k in bad_pixels: + + valid_pairs = np.isfinite(L[:, k]) + + key = valid_pairs.tobytes() + + if key not in pinv_cache: + + Av = A[valid_pairs, :] + + # Need a connected/full-rank network + if ( + Av.shape[0] >= n_unknown + and np.linalg.matrix_rank(Av) == n_unknown + ): + pinv_cache[key] = np.linalg.pinv(Av) + else: + pinv_cache[key] = None + + Ainv = pinv_cache[key] + + if Ainv is not None: + U[:, k] = ( + Ainv + @ L[valid_pairs, k] + ) + + return U + + +def invert_offsets_isce3( + pair_offsets, + output_dir, + reference_date, + block_rows=128, +): + """ + Network-invert pairwise ISCE3 dense offsets. + + Parameters + ---------- + pair_offsets : dict + Dictionary mapping: + + (reference_date, secondary_date) -> dense_offsets path + + Example: + + { + ("20260617", "20260629"): + "/.../20260617_20260629/dense_offsets", + + ("20260629", "20260711"): + "/.../20260629_20260711/dense_offsets", + + ("20260617", "20260711"): + "/.../20260617_20260711/dense_offsets", + } + + Each raster must contain: + + band 1 = azimuth offset + band 2 = range offset + + All rasters must have identical dimensions AND + represent identical Ampcor window locations. + + output_dir : str or Path + Output directory. + + reference_date : str + Stack reference date. Its inverted residual offset + is defined as zero. + + block_rows : int + Number of offset-grid rows processed at once. + + Returns + ------- + dict + Mapping: + + date -> output dense_offsets path + + Output + ------ + output_dir/ + 20260617/dense_offsets + 20260629/dense_offsets + 20260711/dense_offsets + ... + + Each output is a 2-band ENVI Float32 raster: + + band 1 = azimuth residual relative to reference + band 2 = range residual relative to reference + """ + + output_dir = Path(output_dir) + output_dir.mkdir( + parents=True, + exist_ok=True, + ) + + # Normalize inputs + pair_offsets = { + (str(d1), str(d2)): Path(path) + for (d1, d2), path in pair_offsets.items() + } + + pairs = list(pair_offsets.keys()) + + if len(pairs) == 0: + raise ValueError( + "No pair offset rasters supplied." + ) + + # -------------------------------------------------- + # Build date list + # -------------------------------------------------- + + dates = sorted( + { + date + for pair in pairs + for date in pair + } + ) + + reference_date = str(reference_date) + + if reference_date not in dates: + raise ValueError( + f"Reference date {reference_date} " + f"is not in network: {dates}" + ) + + print("Dates:") + for date in dates: + marker = " <-- reference" if date == reference_date else "" + print(f" {date}{marker}") + + print() + print("Pairs:") + for pair in pairs: + print( + f" {pair[0]} -> {pair[1]}" + ) + + # -------------------------------------------------- + # Construct full incidence matrix + # + # pair measurement: + # + # offset_ij = offset_j - offset_i + # + # -------------------------------------------------- + + n_pairs = len(pairs) + n_dates = len(dates) + + date_index = { + date: i + for i, date in enumerate(dates) + } + + A_full = np.zeros( + (n_pairs, n_dates), + dtype=np.float64, + ) + + for k, (d1, d2) in enumerate(pairs): + + A_full[ + k, + date_index[d1] + ] = -1.0 + + A_full[ + k, + date_index[d2] + ] = +1.0 + + # -------------------------------------------------- + # Fix reference date = 0 + # -------------------------------------------------- + + unknown_dates = [ + date + for date in dates + if date != reference_date + ] + + unknown_columns = [ + date_index[date] + for date in unknown_dates + ] + + A = A_full[:, unknown_columns] + + print() + print("Design matrix:") + print(A) + + # -------------------------------------------------- + # Check global network connectivity + # -------------------------------------------------- + + rank = np.linalg.matrix_rank(A) + + if rank != len(unknown_dates): + raise ValueError( + "Offset network is not fully connected to the " + f"reference. rank={rank}, " + f"required={len(unknown_dates)}" + ) + + print( + f"Network rank: {rank}/" + f"{len(unknown_dates)}" + ) + + # -------------------------------------------------- + # Open input rasters + # -------------------------------------------------- + + input_ds = [] + + width = None + length = None + + for pair in pairs: + + path = pair_offsets[pair] + + ds = gdal.Open( + str(path), + gdal.GA_ReadOnly, + ) + + if ds is None: + raise RuntimeError( + f"Could not open {path}" + ) + + if ds.RasterCount < 2: + raise ValueError( + f"{path} has {ds.RasterCount} bands; " + "expected at least 2." + ) + + if width is None: + width = ds.RasterXSize + length = ds.RasterYSize + + elif ( + ds.RasterXSize != width + or ds.RasterYSize != length + ): + raise ValueError( + "All pair offsets must have identical " + "dimensions.\n" + f"Expected: {length} x {width}\n" + f"{path}: " + f"{ds.RasterYSize} x " + f"{ds.RasterXSize}" + ) + + input_ds.append(ds) + + print() + print( + f"Offset grid: {length} x {width}" + ) + + # -------------------------------------------------- + # Create output rasters + # -------------------------------------------------- + + driver = gdal.GetDriverByName("ENVI") + + output_ds = {} + output_paths = {} + + for date in dates: + + date_dir = output_dir / date + date_dir.mkdir( + parents=True, + exist_ok=True, + ) + + path = date_dir / "dense_offsets" + + ds = driver.Create( + str(path), + width, + length, + 2, + gdal.GDT_Float32, + options=[ + "INTERLEAVE=BIP", + ], + ) + + ds.GetRasterBand(1).SetDescription( + "azimuth_offset" + ) + + ds.GetRasterBand(2).SetDescription( + "range_offset" + ) + + output_ds[date] = ds + output_paths[date] = path + + # -------------------------------------------------- + # Process by row blocks + # -------------------------------------------------- + + for y0 in range( + 0, + length, + block_rows, + ): + + nrows = min( + block_rows, + length - y0, + ) + + print( + f"Rows {y0}:{y0 + nrows} " + f"of {length}" + ) + + # Pair stacks: + # + # (n_pairs, nrows, width) + # + + az_pairs = np.empty( + (n_pairs, nrows, width), + dtype=np.float64, + ) + + rg_pairs = np.empty_like( + az_pairs + ) + + for k, ds in enumerate(input_ds): + + az_pairs[k] = ( + ds.GetRasterBand(1) + .ReadAsArray( + 0, + y0, + width, + nrows, + ) + ) + + rg_pairs[k] = ( + ds.GetRasterBand(2) + .ReadAsArray( + 0, + y0, + width, + nrows, + ) + ) + + # Flatten spatial dimensions + # + # (n_pairs, n_pixels) + # + + az_L = az_pairs.reshape( + n_pairs, + -1, + ) + + rg_L = rg_pairs.reshape( + n_pairs, + -1, + ) + + # ---------------------------------------------- + # Network inversion + # ---------------------------------------------- + + az_U = _solve_offset_block( + A, + az_L, + ) + + rg_U = _solve_offset_block( + A, + rg_L, + ) + + # ---------------------------------------------- + # Reference = zero + # ---------------------------------------------- + + ref_az = np.zeros( + (nrows, width), + dtype=np.float32, + ) + + ref_rg = np.zeros_like( + ref_az + ) + + output_ds[ + reference_date + ].GetRasterBand(1).WriteArray( + ref_az, + 0, + y0, + ) + + output_ds[ + reference_date + ].GetRasterBand(2).WriteArray( + ref_rg, + 0, + y0, + ) + + # ---------------------------------------------- + # Other dates + # ---------------------------------------------- + + for i, date in enumerate( + unknown_dates + ): + + az = az_U[i].reshape( + nrows, + width, + ).astype(np.float32) + + rg = rg_U[i].reshape( + nrows, + width, + ).astype(np.float32) + + output_ds[ + date + ].GetRasterBand(1).WriteArray( + az, + 0, + y0, + ) + + output_ds[ + date + ].GetRasterBand(2).WriteArray( + rg, + 0, + y0, + ) + + # -------------------------------------------------- + # Close datasets + # -------------------------------------------------- + + for ds in output_ds.values(): + band1 = ds.GetRasterBand(1) + band2 = ds.GetRasterBand(2) + band1.FlushCache() + band2.FlushCache() + ds.FlushCache() + output_ds.clear() + input_ds.clear() + + print() + print("Offset inversion complete.") + + for date in dates: + print( + f"{date}: {output_paths[date]}" + ) + + return output_paths + + +def _add_one_offset( + offset1, + offset2, + output, + block_rows=1024, + invalid_threshold=1.0e5, + invalid_value=-1.0e6, +): + """ + Add two single-band GDAL-readable offset rasters and + write a Float64 ISCE raster. + """ + + offset1 = Path(offset1) + offset2 = Path(offset2) + output = Path(output) + + ds1 = gdal.Open( + str(offset1), + gdal.GA_ReadOnly, + ) + + if ds1 is None: + raise RuntimeError( + f"Could not open {offset1}" + ) + + ds2 = gdal.Open( + str(offset2), + gdal.GA_ReadOnly, + ) + + if ds2 is None: + ds1 = None + + raise RuntimeError( + f"Could not open {offset2}" + ) + + if ds1.RasterCount != 1: + raise ValueError( + f"{offset1} has {ds1.RasterCount} bands; " + "expected 1." + ) + + if ds2.RasterCount != 1: + raise ValueError( + f"{offset2} has {ds2.RasterCount} bands; " + "expected 1." + ) + + width1 = ds1.RasterXSize + length1 = ds1.RasterYSize + + width2 = ds2.RasterXSize + length2 = ds2.RasterYSize + + if ( + width1 != width2 + or length1 != length2 + ): + raise ValueError( + "Offset dimensions do not match:\n" + f" {offset1}: {length1} x {width1}\n" + f" {offset2}: {length2} x {width2}" + ) + + width = width1 + length = length1 + + band1 = ds1.GetRasterBand(1) + band2 = ds2.GetRasterBand(1) + + nodata1 = band1.GetNoDataValue() + nodata2 = band2.GetNoDataValue() + + dtype1 = gdal.GetDataTypeName( + band1.DataType + ) + + dtype2 = gdal.GetDataTypeName( + band2.DataType + ) + + print() + print("Adding:") + print(" input 1:", offset1) + print(" dtype:", dtype1) + print(" input 2:", offset2) + print(" dtype:", dtype2) + print(" output :", output) + print(" dtype: Float64") + print( + " shape :", + f"{length} x {width}", + ) + + output.parent.mkdir( + parents=True, + exist_ok=True, + ) + + driver = gdal.GetDriverByName( + "ISCE" + ) + + if driver is None: + raise RuntimeError( + "GDAL ISCE driver is not available." + ) + + out_ds = driver.Create( + str(output), + width, + length, + 1, + gdal.GDT_Float64, + ) + + if out_ds is None: + raise RuntimeError( + f"Could not create {output}" + ) + + out_band = out_ds.GetRasterBand(1) + + for y0 in range( + 0, + length, + block_rows, + ): + + nrows = min( + block_rows, + length - y0, + ) + + y1 = y0 + nrows + + print( + f" rows {y0}:{y1} / {length}" + ) + + a = band1.ReadAsArray( + 0, + y0, + width, + nrows, + ).astype( + np.float64, + copy=False, + ) + + b = band2.ReadAsArray( + 0, + y0, + width, + nrows, + ).astype( + np.float64, + copy=False, + ) + + invalid = ( + ~np.isfinite(a) + | ~np.isfinite(b) + ) + + if nodata1 is not None: + invalid |= ( + a == nodata1 + ) + + if nodata2 is not None: + invalid |= ( + b == nodata2 + ) + + if invalid_threshold is not None: + + invalid |= ( + np.abs(a) + >= invalid_threshold + ) + + invalid |= ( + np.abs(b) + >= invalid_threshold + ) + + result = a + b + + result[invalid] = ( + invalid_value + ) + + out_band.WriteArray( + result, + 0, + y0, + ) + + out_band.FlushCache() + out_ds.FlushCache() + + ds1 = None + ds2 = None + out_ds = None + + +def add_offset_folders( + folder1, + folder2, + outfolder, + block_rows=1024, + invalid_threshold=1.0e5, + invalid_value=-1.0e6, +): + """ + Add azimuth.off and range.off from two folders. + + Expected input: + + folder1/ + azimuth.off + range.off + + folder2/ + azimuth.off + range.off + + Output: + + outfolder/ + azimuth.off + azimuth.off.xml + range.off + range.off.xml + + Operation: + + out/azimuth.off = + folder1/azimuth.off + + + folder2/azimuth.off + + out/range.off = + folder1/range.off + + + folder2/range.off + """ + + folder1 = Path(folder1) + folder2 = Path(folder2) + outfolder = Path(outfolder) + + outfolder.mkdir( + parents=True, + exist_ok=True, + ) + + az1 = folder1 / "azimuth.off" + rg1 = folder1 / "range.off" + + az2 = folder2 / "azimuth.off" + rg2 = folder2 / "range.off" + + out_az = outfolder / "azimuth.off" + out_rg = outfolder / "range.off" + + required = [ + az1, + rg1, + az2, + rg2, + ] + + missing = [ + path + for path in required + if not path.exists() + ] + + if missing: + raise FileNotFoundError( + "Missing required offset files:\n" + + "\n".join( + f" {path}" + for path in missing + ) + ) + + print("Folder 1:", folder1) + print("Folder 2:", folder2) + print("Output :", outfolder) + + # -------------------------------------------- + # Azimuth + # -------------------------------------------- + + _add_one_offset( + offset1=az1, + offset2=az2, + output=out_az, + block_rows=block_rows, + invalid_threshold=invalid_threshold, + invalid_value=invalid_value, + ) + + # -------------------------------------------- + # Range + # -------------------------------------------- + + _add_one_offset( + offset1=rg1, + offset2=rg2, + output=out_rg, + block_rows=block_rows, + invalid_threshold=invalid_threshold, + invalid_value=invalid_value, + ) + + print() + print("Done.") + print("Azimuth:", out_az) + print("Range :", out_rg) + + return { + "azimuth": out_az, + "range": out_rg, + } diff --git a/share/stack/make_stack.py b/share/stack/make_stack.py new file mode 100644 index 0000000000..9b34d085cb --- /dev/null +++ b/share/stack/make_stack.py @@ -0,0 +1,122 @@ +#!/usr/bin/env python3 + +import argparse + +from StackProcessor import StackProcessor + + +STEPS = [ + "coarse", + "dense", + "invert", + "rubbersheet", + "resample", +] + + +def parse_args(): + parser = argparse.ArgumentParser( + description="Run NISAR StackProcessor pipeline." + ) + + parser.add_argument( + "--slc-folder", + required=True, + help="Folder containing NISAR RSLC HDF5 files.", + ) + + parser.add_argument( + "--out-folder", + required=True, + help="Output stack folder.", + ) + + parser.add_argument( + "--dem-file", + required=True, + help="DEM file.", + ) + + parser.add_argument( + "--pair-level", + type=int, + default=2, + help="Dense-offset pair level (default: 2).", + ) + + parser.add_argument( + "--freq", + default="A", + help="Frequency, e.g. A (default: A).", + ) + + parser.add_argument( + "--pol", + default="HH", + help="Polarization (default: HH).", + ) + + parser.add_argument( + "--overwrite", + action="store_true", + help="Overwrite existing outputs.", + ) + + parser.add_argument( + "--debug", + action="store_true", + help="Enable StackProcessor debug mode.", + ) + + parser.add_argument( + "--steps", + nargs="+", + choices=STEPS, + default=STEPS, + help=( + "Pipeline stages to run. " + "Default: coarse dense invert rubbersheet resample" + ), + ) + + return parser.parse_args() + + +def main(): + args = parse_args() + processor = StackProcessor( + slc_folder=args.slc_folder, + out_folder=args.out_folder, + dem_file=args.dem_file, + pair_level=args.pair_level, + freq=args.freq, + pol=args.pol, + overwrite=args.overwrite, + debug=args.debug, + ) + + if "coarse" in args.steps: + print("\n===== COARSE REGISTER =====") + processor.coarse_register_scans() + + if "dense" in args.steps: + print("\n===== DENSE OFFSETS =====") + processor.dense_offset_pairs() + + if "invert" in args.steps: + print("\n===== INVERT DENSE OFFSETS =====") + processor.invert_dense_offset_pairs() + + if "rubbersheet" in args.steps: + print("\n===== RUBBERSHEET =====") + processor.rubbersheet() + + if "resample" in args.steps: + print("\n===== RESAMPLE SLCs =====") + processor.resample_slcs() + + print("\n===== DONE =====") + + +if __name__ == "__main__": + main() diff --git a/share/stack/offset_util.py b/share/stack/offset_util.py new file mode 100644 index 0000000000..0d18f3a6d8 --- /dev/null +++ b/share/stack/offset_util.py @@ -0,0 +1,265 @@ +from pathlib import Path +import isce3 +from osgeo import gdal +from nisar.workflows.dense_offsets import ( + set_optional_attributes, + create_empty_dataset, +) +from nisar.workflows.helpers import (get_cfg_freq_pols, + get_offset_radar_grid) +import h5py +import numpy as np +from nisar.products.readers import SLC +from nisar.products.insar.product_paths import RIFGGroupsPaths + +def create_minimal_rifg_for_rubbersheet( + out_h5, + cfg, + overwrite=True, +): + """ + Create only the HDF5 datasets required by nisar.workflows.rubbersheet. + + Parameters + ---------- + out_h5 : str or Path + Output minimal RIFG HDF5. + cfg : dict + cfg = ref_config["runconfig"]["groups"] + overwrite : bool + Replace existing file. + """ + + out_h5 = Path(out_h5) + out_h5.parent.mkdir(parents=True, exist_ok=True) + + if out_h5.exists(): + if overwrite: + out_h5.unlink() + else: + raise FileExistsError(out_h5) + + ref_file = cfg["input_file_group"]["reference_rslc_file"] + ref_slc = SLC(hdf5file=ref_file) + + with h5py.File(out_h5, "w", libver="latest") as f: + + for freq, _, pol_list in get_cfg_freq_pols(cfg): + + # Reference radar grid + ref_grid = ref_slc.getRadarGrid(freq) + + # EXACT NISAR pixel-offset grid derived from dense-offset config + off_grid = get_offset_radar_grid(cfg, ref_grid) + + off_length = off_grid.length + off_width = off_grid.width + + # Pixel-offset coordinates. + # These are at the centers of the Ampcor matching windows. + slant_range = ( + off_grid.starting_range + + np.arange(off_width, dtype=np.float64) + * off_grid.range_pixel_spacing + ) + + zero_doppler_time = ( + off_grid.sensing_start + + np.arange(off_length, dtype=np.float64) + / off_grid.prf + ) + + base = ( + f"{RIFGGroupsPaths().SwathsPath}" + f"/frequency{freq}/pixelOffsets" + ) + + g = f.require_group(base) + + # ---------------------------------------------------------- + # Required coordinate vectors + # ---------------------------------------------------------- + + ds = g.create_dataset( + "slantRange", + data=slant_range, + dtype=np.float64, + ) + ds.attrs["units"] = np.bytes_("meters") + + ds = g.create_dataset( + "zeroDopplerTime", + data=zero_doppler_time, + dtype=np.float64, + ) + ds.attrs["units"] = np.bytes_("seconds") + + # Not strictly required by rubbersheet, but useful + ds = g.create_dataset( + "slantRangeSpacing", + data=np.float64(off_grid.range_pixel_spacing), + ) + ds.attrs["units"] = np.bytes_("meters") + + ds = g.create_dataset( + "zeroDopplerTimeSpacing", + data=np.float64(1.0 / off_grid.prf), + ) + ds.attrs["units"] = np.bytes_("seconds") + + # ---------------------------------------------------------- + # Optional mask + # + # rubbersheet reads this only when: + # subswath_mask_apply_enabled == True + # + # 255 matches the native writer's default fill value. + # ---------------------------------------------------------- + + g.create_dataset( + "mask", + shape=(off_length, off_width), + dtype=np.uint32, + fillvalue=np.uint32(255), + ) + + # ---------------------------------------------------------- + # Datasets that rubbersheet writes into + # ---------------------------------------------------------- + + for pol in pol_list: + + pg = g.require_group(pol) + + ds = pg.create_dataset( + "alongTrackOffset", + shape=(off_length, off_width), + dtype=np.float32, + fillvalue=np.nan, + ) + ds.attrs["units"] = np.bytes_("meters") + + ds = pg.create_dataset( + "slantRangeOffset", + shape=(off_length, off_width), + dtype=np.float32, + fillvalue=np.nan, + ) + ds.attrs["units"] = np.bytes_("meters") + + ds = pg.create_dataset( + "correlationSurfacePeak", + shape=(off_length, off_width), + dtype=np.float32, + fillvalue=np.nan, + ) + ds.attrs["units"] = np.bytes_("1") + + print( + f"frequency{freq}: " + f"{off_length} x {off_width}, " + f"pols={pol_list}" + ) + + print(f"Created minimal rubbersheet HDF5: {out_h5}") + + return str(out_h5) + +def set_nested(d, keys, value): + for key in keys[:-1]: + d = d.setdefault(key, {}) + d[keys[-1]] = value + +def dense_offset_coregistered( + reference_slc, + secondary_slc, + out_dir, + dense_cfg, + gpu_id=0, +): + reference_slc = Path(reference_slc) + secondary_slc = Path(secondary_slc) + out_dir = Path(out_dir) + out_dir.mkdir(parents=True, exist_ok=True) + + # GPU + device = isce3.cuda.core.Device(gpu_id) + isce3.cuda.core.set_device(device) + + ampcor = isce3.cuda.matchtemplate.PyCuAmpcor() + ampcor.deviceID = gpu_id + ampcor.useMmap = 1 + + # Read raster dimensions + ref_raster = isce3.io.Raster(str(reference_slc)) + sec_raster = isce3.io.Raster(str(secondary_slc)) + + if ( + ref_raster.length != sec_raster.length + or ref_raster.width != sec_raster.width + ): + raise ValueError( + f"Raster dimensions differ: " + f"ref={ref_raster.length}x{ref_raster.width}, " + f"sec={sec_raster.length}x{sec_raster.width}" + ) + + # Input files + ampcor.referenceImageName = str(reference_slc) + ampcor.referenceImageHeight = ref_raster.length + ampcor.referenceImageWidth = ref_raster.width + + ampcor.secondaryImageName = str(secondary_slc) + ampcor.secondaryImageHeight = sec_raster.length + ampcor.secondaryImageWidth = sec_raster.width + + # Apply normal NISAR dense-offset parameters + ampcor = set_optional_attributes( + ampcor, + dense_cfg, + ref_raster.length, + ref_raster.width, + ) + + # Outputs + ampcor.offsetImageName = str(out_dir / "dense_offsets") + ampcor.grossOffsetImageName = str(out_dir / "gross_offset") + ampcor.snrImageName = str(out_dir / "snr") + ampcor.covImageName = str(out_dir / "covariance") + ampcor.corrImageName = str(out_dir / "correlation_peak") + + nx = ampcor.numberWindowAcross + ny = ampcor.numberWindowDown + + create_empty_dataset( + str(out_dir / "dense_offsets"), + nx, ny, 2, gdal.GDT_Float32, + ) + + create_empty_dataset( + str(out_dir / "gross_offset"), + nx, ny, 2, gdal.GDT_Float32, + ) + + create_empty_dataset( + str(out_dir / "snr"), + nx, ny, 1, gdal.GDT_Float32, + ) + + create_empty_dataset( + str(out_dir / "covariance"), + nx, ny, 3, gdal.GDT_Float32, + ) + + create_empty_dataset( + str(out_dir / "correlation_peak"), + nx, ny, 1, gdal.GDT_Float32, + ) + + print( + f"Dense offset grid: {ny} × {nx}" + ) + + ampcor.runAmpcor() + + return out_dir / "dense_offsets" \ No newline at end of file diff --git a/share/stack/stack_utils.py b/share/stack/stack_utils.py new file mode 100644 index 0000000000..7bb8d6fa3d --- /dev/null +++ b/share/stack/stack_utils.py @@ -0,0 +1,161 @@ +import os +from pathlib import Path +from nisar.workflows.resample_slc_v2 import ( + resample_secondary_rslc_onto_reference +) +import shutil +from nisar.products.insar import ROFFWriter +import yaml +from argparse import Namespace +from ScanList import * +from nisar.workflows.runconfig import RunConfig + +MODULE_DIR = Path(__file__).resolve().parent +template = f'{MODULE_DIR}/insar.yaml' +with open(template) as f: + template_yaml = yaml.safe_load(f) + + +def load_config(template=template): + args = Namespace( + run_config_path=template, + log_file=False, + restart=True, + ) + + insar_runcfg = RunConfig(args, workflow_name="insar") + insar_runcfg.load_yaml_to_dict() + return insar_runcfg + + +def set_nested(d, keys, value): + for key in keys[:-1]: + d = d.setdefault(key, {}) + d[keys[-1]] = value + + +def resample( + out_dir, + offset_dir, + ref_path, + secondary_path, + out_tag="coarse", + freq="A", + pols=["HH"], + block_size_az=1024, + block_size_rg=0, + with_gpu=True, +): + os.makedirs(out_dir, exist_ok=True) + resample_secondary_rslc_onto_reference( + ref_file_path=ref_path, + sec_file_path=secondary_path, + out_path=out_dir, + az_off_file=offset_dir / "azimuth.off", + rg_off_file=offset_dir / "range.off", + freq=freq, + pols=pols, + block_size_az=block_size_az, + block_size_rg=block_size_rg, + with_gpu=with_gpu, + ) + shutil.move( + out_dir / + Path('HH/coregistered_secondary.hdr'), + out_dir / + Path(f'{out_tag}.hdr')) + shutil.move( + out_dir / + Path('HH/coregistered_secondary.slc'), + out_dir / + Path(f'{out_tag}.slc')) + shutil.rmtree(out_dir / 'HH') + + +def create_roff_h5(output_hdf5, cfg): + with ROFFWriter(name=output_hdf5, mode='w', + runconfig_dict=cfg, + runconfig_path="None") as roff: + roff.save_to_hdf5() + + +def relative_symlink_dir(target_dir, link_dir): + target_dir = Path(target_dir).resolve() + link_dir = Path(link_dir) + + link_dir.parent.mkdir(parents=True, exist_ok=True) + + if link_dir.is_symlink(): + link_dir.unlink() + elif link_dir.exists(): + raise FileExistsError( + f"{link_dir} already exists and is not a symlink") + + relative_target = os.path.relpath( + target_dir, + start=link_dir.parent.resolve() + ) + + link_dir.symlink_to( + relative_target, + target_is_directory=True + ) + + return link_dir + + +def relative_symlink_contents(target_dir, link_dir): + """ + Symlink all top-level contents of target_dir into link_dir. + + Existing destination files/directories/symlinks with the same + names are overwritten. + + Symlinks are relative. + """ + + target_dir = Path(target_dir).resolve() + link_dir = Path(link_dir) + + if not target_dir.is_dir(): + raise NotADirectoryError(target_dir) + + # If link_dir itself is currently a symlink or file, + # remove it and create a real directory. + if link_dir.is_symlink() or link_dir.is_file(): + link_dir.unlink() + + link_dir.mkdir( + parents=True, + exist_ok=True, + ) + + for src in target_dir.iterdir(): + + dst = link_dir / src.name + + # ---------------------------------------- + # Remove existing destination + # ---------------------------------------- + + if dst.is_symlink() or dst.is_file(): + dst.unlink() + + elif dst.is_dir(): + shutil.rmtree(dst) + + # ---------------------------------------- + # Create relative symlink + # ---------------------------------------- + + relative_target = os.path.relpath( + src, + start=link_dir.resolve(), + ) + + dst.symlink_to( + relative_target, + target_is_directory=src.is_dir(), + ) + + return link_dir diff --git a/tests/python/extensions/pybind/geometry/rdr2geo.py b/tests/python/extensions/pybind/geometry/rdr2geo.py index daceb6c5d3..afa29ba7f5 100644 --- a/tests/python/extensions/pybind/geometry/rdr2geo.py +++ b/tests/python/extensions/pybind/geometry/rdr2geo.py @@ -1,6 +1,5 @@ #!/usr/bin/env python3 -import itertools import os import types diff --git a/tests/python/extensions/pybind/unwrap/phass.py b/tests/python/extensions/pybind/unwrap/phass.py index ea93149952..24bde34e38 100644 --- a/tests/python/extensions/pybind/unwrap/phass.py +++ b/tests/python/extensions/pybind/unwrap/phass.py @@ -2,7 +2,6 @@ Unit test for Phass unwrapper ''' -import os import numpy as np import numpy.testing as npt diff --git a/tests/python/packages/estimate_dop.py b/tests/python/packages/estimate_dop.py index c9a371a5ac..8755eafe89 100644 --- a/tests/python/packages/estimate_dop.py +++ b/tests/python/packages/estimate_dop.py @@ -13,7 +13,6 @@ """ import os -import sys class EstimateDoptest(object): def __init__(self): diff --git a/tests/python/packages/isce3/core/gpu_check.py b/tests/python/packages/isce3/core/gpu_check.py index 5525851869..6df2a7f4d5 100644 --- a/tests/python/packages/isce3/core/gpu_check.py +++ b/tests/python/packages/isce3/core/gpu_check.py @@ -1,4 +1,3 @@ -import os import numpy.testing as npt diff --git a/tests/python/packages/isce3/io/gdal/gdal_raster.py b/tests/python/packages/isce3/io/gdal/gdal_raster.py index eff9e86209..faadbf7fd4 100644 --- a/tests/python/packages/isce3/io/gdal/gdal_raster.py +++ b/tests/python/packages/isce3/io/gdal/gdal_raster.py @@ -1,7 +1,6 @@ import os import pytest import numpy as np -from dataclasses import dataclass from osgeo import gdal from pathlib import Path from tempfile import TemporaryDirectory diff --git a/tools/imagesets/imgset.py b/tools/imagesets/imgset.py index 1c2acf0b4e..340147ab99 100644 --- a/tools/imagesets/imgset.py +++ b/tools/imagesets/imgset.py @@ -1,4 +1,10 @@ -import os, subprocess, sys, shutil, stat, logging, shlex, getpass +import os +import subprocess +import sys +import shutil +import logging +import shlex +import getpass from datetime import datetime, timezone from pathlib import Path from textwrap import dedent diff --git a/tools/run.py b/tools/run.py index 1bf855f183..285715c6c2 100755 --- a/tools/run.py +++ b/tools/run.py @@ -1,5 +1,5 @@ #!/usr/bin/env python3 -import argparse, os +import argparse from imagesets import imagesets, projsrcdir def run(*, steps, imgset, **kwargs): From 016341ced49abe3449ef61764defde450b5f98b6 Mon Sep 17 00:00:00 2001 From: william Wu Date: Wed, 2 Sep 2026 10:23:00 +0800 Subject: [PATCH 2/2] 1.added date and baseline threshold dependent network formation similar with isce2 2.standardized all folder and file path definition in StackProcessor --- share/stack/Igram.py | 570 ++++++++++++++ share/stack/StackProcessor.py | 576 ++++++++++---- share/stack/iono.py | 1358 +++++++++++++++++++++++++++++++++ 3 files changed, 2358 insertions(+), 146 deletions(-) create mode 100644 share/stack/Igram.py create mode 100644 share/stack/iono.py diff --git a/share/stack/Igram.py b/share/stack/Igram.py new file mode 100644 index 0000000000..6be8b84831 --- /dev/null +++ b/share/stack/Igram.py @@ -0,0 +1,570 @@ +""" +isce3_pair_igram +================ + +Pair-level interferogram workflow for two already-aligned SLC rasters. + +Provides: +- ISCE3 Crossmul +- wrapped interferogram +- coherence +- SNAPHU, ICU, or Whirlwind unwrapping +- connected components + +Example +------- +from isce3_pair_igram import ( + PairIgramWorkflow, + IgramConfig, + RadarMeta, +) + +cfg = IgramConfig( + range_looks=11, + azimuth_looks=11, + gpu=True, + overwrite=True, +) + +wf = PairIgramWorkflow( + ref_slc="/path/reference.slc", + sec_slc="/path/coregistered_secondary.slc", + outdir="/path/pair_igram", + metadata_rslc="/path/reference_rslc.h5", + ref_rslc="/path/reference_rslc.h5", + sec_rslc="/path/secondary_rslc.h5", + config=cfg, +) + +products = wf.run() +""" + +from __future__ import annotations + +import shutil +from dataclasses import dataclass +from pathlib import Path +from typing import Optional + +import numpy as np +from osgeo import gdal + +import isce3 +from isce3.splitspectrum.splitspectrum import BandpassMetaData +from nisar.products.readers import RSLC + +gdal.UseExceptions() + +C = float(isce3.core.speed_of_light) + + +@dataclass +class RadarMeta: + center_frequency: float + range_bandwidth: float + range_sample_frequency: float + range_pixel_spacing: float + starting_range: float + wavelength: float + + @classmethod + def from_nisar_rslc(cls, rslc_path: str | Path, frequency: str = "A"): + slc = RSLC(hdf5file=str(rslc_path)) + meta = BandpassMetaData.load_from_slc( + slc_product=slc, + freq=frequency, + ) + grid = slc.getRadarGrid(frequency) + + return cls( + center_frequency=float(meta.center_freq), + range_bandwidth=float(meta.rg_bandwidth), + range_sample_frequency=float(meta.rg_sample_freq), + range_pixel_spacing=float(meta.rg_pxl_spacing), + starting_range=float(grid.starting_range), + wavelength=float(meta.wavelength), + ) + + @classmethod + def from_values( + cls, + *, + center_frequency: float, + range_bandwidth: float, + range_sample_frequency: float, + starting_range: float = 0.0, + ): + fs = float(range_sample_frequency) + dr = C / (2.0 * fs) + f0 = float(center_frequency) + + return cls( + center_frequency=f0, + range_bandwidth=float(range_bandwidth), + range_sample_frequency=fs, + range_pixel_spacing=dr, + starting_range=float(starting_range), + wavelength=C / f0, + ) + + def slant_range(self, index: int) -> float: + return self.starting_range + float(index) * self.range_pixel_spacing + + +@dataclass +class IgramConfig: + # Crossmul + range_looks: int = 11 + azimuth_looks: int = 11 + oversample: int = 2 + lines_per_block: int = 1024 + + gpu: bool = False + gpu_id: int = 0 + + # Optional flattening + flatten_range_offset: Optional[str] = None + starting_range_shift: float = 0.0 + + # Unwrapping: "snaphu", "icu", or "whirlwind" + unwrap_algorithm: str = "snaphu" + + # Whirlwind (package: whirlwind-insar; import: whirlwind) + whirlwind_nlooks: Optional[float] = None + whirlwind_downsample: Optional[int] = None + + snaphu_nlooks: Optional[float] = None + snaphu_cost: str = "smooth" + snaphu_init: str = "mcf" + snaphu_min_conncomp_frac: float = 0.01 + snaphu_phase_grad_window: tuple[int, int] = (7, 7) + snaphu_ntiles: tuple[int, int] = (1, 1) + snaphu_tile_overlap: tuple[int, int] = (0, 0) + snaphu_nproc: int = 1 + snaphu_tile_cost_thresh: float = 500.0 + snaphu_min_region_size: int = 300 + snaphu_single_tile_reoptimize: bool = True + + overwrite: bool = False + + +@dataclass +class IgramProducts: + wrapped_interferogram: Path + coherence: Path + unwrapped_phase: Path + connected_components: Path + + +class PairIgramWorkflow: + """ + Crossmul + coherence + unwrap for an already-aligned SLC pair. + """ + + def __init__( + self, + *, + ref_slc: str | Path, + sec_slc: str | Path, + outdir: str | Path, + metadata_rslc: str | Path | None = None, + radar_meta: RadarMeta | None = None, + ref_rslc: str | Path | None = None, + sec_rslc: str | Path | None = None, + frequency: str = "A", + config: IgramConfig | None = None, + ): + self.ref_slc = Path(ref_slc).resolve() + self.sec_slc = Path(sec_slc).resolve() + self.outdir = Path(outdir).resolve() + self.outdir.mkdir(parents=True, exist_ok=True) + + self.frequency = frequency + self.config = config or IgramConfig() + + self.ref_rslc = Path(ref_rslc).resolve() if ref_rslc else None + self.sec_rslc = Path(sec_rslc).resolve() if sec_rslc else None + + if metadata_rslc is not None and radar_meta is not None: + raise ValueError("Provide either metadata_rslc or radar_meta, not both.") + + if radar_meta is not None: + self.meta = radar_meta + elif metadata_rslc is not None: + self.meta = RadarMeta.from_nisar_rslc( + metadata_rslc, + frequency=frequency, + ) + else: + raise ValueError("Provide metadata_rslc or radar_meta.") + + self.width, self.length = self._validate_aligned() + self.dopplers = self._load_dopplers() + + def run(self) -> IgramProducts: + ifg, coh = self.crossmul() + unw, cc = self.unwrap() + + return IgramProducts( + wrapped_interferogram=ifg, + coherence=coh, + unwrapped_phase=unw, + connected_components=cc, + ) + + def crossmul(self) -> tuple[Path, Path]: + cfg = self.config + + ifg_path = self.outdir / "wrappedInterferogram.int" + coh_path = self.outdir / "coherenceMagnitude.cor" + + self._prepare_output(ifg_path) + self._prepare_output(coh_path) + + ref = isce3.io.Raster(str(self.ref_slc)) + sec = isce3.io.Raster(str(self.sec_slc)) + + if ref.width != sec.width or ref.length != sec.length: + raise ValueError("Reference and secondary SLC grids differ.") + + out_width = ref.width // cfg.range_looks + out_length = ref.length // cfg.azimuth_looks + + ifg = isce3.io.Raster( + str(ifg_path), + out_width, + out_length, + 1, + gdal.GDT_CFloat32, + "ENVI", + ) + coh = isce3.io.Raster( + str(coh_path), + out_width, + out_length, + 1, + gdal.GDT_Float32, + "ENVI", + ) + + cm = self._new_crossmul() + + flatten = ( + isce3.io.Raster(str(cfg.flatten_range_offset)) + if cfg.flatten_range_offset + else None + ) + + cm.crossmul( + ref, + sec, + ifg, + coh, + flatten, + ) + + del ifg + del coh + del ref + del sec + + if flatten is not None: + del flatten + + return ifg_path, coh_path + + def unwrap(self) -> tuple[Path, Path]: + algo = self.config.unwrap_algorithm.lower() + + ifg = self.outdir / "wrappedInterferogram.int" + coh = self.outdir / "coherenceMagnitude.cor" + + if algo == "snaphu": + return self._unwrap_snaphu(ifg, coh) + + if algo == "icu": + return self._unwrap_icu(ifg, coh) + + if algo == "whirlwind": + return self._unwrap_whirlwind(ifg, coh) + + raise ValueError( + f"Unsupported unwrap algorithm: {algo!r}. " + "Choose 'snaphu', 'icu', or 'whirlwind'." + ) + + def _new_crossmul(self): + cfg = self.config + + if cfg.gpu: + device = isce3.cuda.core.Device(cfg.gpu_id) + isce3.cuda.core.set_device(device) + cm = isce3.cuda.signal.Crossmul() + else: + cm = isce3.signal.Crossmul() + + cm.range_looks = cfg.range_looks + cm.az_looks = cfg.azimuth_looks + cm.oversample_factor = cfg.oversample + cm.lines_per_block = cfg.lines_per_block + + cm.range_pixel_spacing = self.meta.range_pixel_spacing + cm.wavelength = self.meta.wavelength + + if self.dopplers is not None: + cm.set_dopplers(*self.dopplers) + + if cfg.flatten_range_offset: + cm.ref_sec_offset_starting_range_shift = cfg.starting_range_shift + + return cm + + def _load_dopplers(self): + if self.ref_rslc is None or self.sec_rslc is None: + return None + + ref = RSLC(hdf5file=str(self.ref_rslc)) + sec = RSLC(hdf5file=str(self.sec_rslc)) + + ref_dopp = isce3.core.avg_lut2d_to_lut1d( + ref.getDopplerCentroid(frequency=self.frequency) + ) + sec_dopp = isce3.core.avg_lut2d_to_lut1d( + sec.getDopplerCentroid(frequency=self.frequency) + ) + + return ref_dopp, sec_dopp + + def _unwrap_snaphu( + self, + ifg_path: Path, + coh_path: Path, + ) -> tuple[Path, Path]: + import snaphu + + cfg = self.config + + unw_path = self.outdir / "unwrappedPhase.unw" + cc_path = self.outdir / "connectedComponents.conncomp" + + self._prepare_output(unw_path) + self._prepare_output(cc_path) + + igram = self.read_array(ifg_path, np.complex64) + coh = self.read_array(coh_path, np.float32) + + nlooks = ( + cfg.snaphu_nlooks + if cfg.snaphu_nlooks is not None + else float(cfg.range_looks * cfg.azimuth_looks) + ) + + scratch = self.outdir / "snaphu_scratch" + + if scratch.exists() and cfg.overwrite: + shutil.rmtree(scratch) + + scratch.mkdir(parents=True, exist_ok=True) + + unw, conncomp = snaphu.unwrap( + igram, + coh, + nlooks=nlooks, + cost=cfg.snaphu_cost, + init=cfg.snaphu_init, + min_conncomp_frac=cfg.snaphu_min_conncomp_frac, + phase_grad_window=cfg.snaphu_phase_grad_window, + ntiles=cfg.snaphu_ntiles, + tile_overlap=cfg.snaphu_tile_overlap, + nproc=cfg.snaphu_nproc, + tile_cost_thresh=cfg.snaphu_tile_cost_thresh, + min_region_size=cfg.snaphu_min_region_size, + single_tile_reoptimize=cfg.snaphu_single_tile_reoptimize, + scratchdir=scratch, + ) + + self.write_array( + unw_path, + np.asarray(unw, dtype=np.float32), + gdal.GDT_Float32, + ) + self.write_array( + cc_path, + np.asarray(conncomp, dtype=np.uint32), + gdal.GDT_UInt32, + ) + + return unw_path, cc_path + + def _unwrap_whirlwind( + self, + ifg_path: Path, + coh_path: Path, + ) -> tuple[Path, Path]: + """Unwrap with the ``whirlwind-insar`` Python package.""" + try: + import whirlwind as ww + except ImportError as exc: + raise ImportError( + "Whirlwind unwrapping requires the 'whirlwind-insar' package. " + "Install it with `pip install whirlwind-insar` or " + "`conda install -c conda-forge whirlwind-insar`." + ) from exc + + cfg = self.config + unw_path = self.outdir / "unwrappedPhase.unw" + cc_path = self.outdir / "connectedComponents.conncomp" + + self._prepare_output(unw_path) + self._prepare_output(cc_path) + + igram = self.read_array(ifg_path, np.complex64) + coh = self.read_array(coh_path, np.float32) + + nlooks = ( + cfg.whirlwind_nlooks + if cfg.whirlwind_nlooks is not None + else float(cfg.range_looks * cfg.azimuth_looks) + ) + + kwargs = {"nlooks": float(nlooks)} + if cfg.whirlwind_downsample is not None: + if cfg.whirlwind_downsample < 1: + raise ValueError("whirlwind_downsample must be >= 1") + kwargs["downsample"] = int(cfg.whirlwind_downsample) + + unw, conncomp = ww.unwrap(igram, coh, **kwargs) + + self.write_array( + unw_path, + np.asarray(unw, dtype=np.float32), + gdal.GDT_Float32, + ) + self.write_array( + cc_path, + np.asarray(conncomp, dtype=np.uint32), + gdal.GDT_UInt32, + ) + + return unw_path, cc_path + + def _unwrap_icu( + self, + ifg_path: Path, + coh_path: Path, + ) -> tuple[Path, Path]: + unw_path = self.outdir / "unwrappedPhase.unw" + cc_path = self.outdir / "connectedComponents.conncomp" + + self._prepare_output(unw_path) + self._prepare_output(cc_path) + + igram = isce3.io.Raster(str(ifg_path)) + coh = isce3.io.Raster(str(coh_path)) + + unw = isce3.io.Raster( + str(unw_path), + igram.width, + igram.length, + 1, + gdal.GDT_Float32, + "ENVI", + ) + cc = isce3.io.Raster( + str(cc_path), + igram.width, + igram.length, + 1, + gdal.GDT_Byte, + "ENVI", + ) + + icu = isce3.unwrap.ICU() + icu.unwrap(unw, cc, igram, coh) + + del unw + del cc + del igram + del coh + + return unw_path, cc_path + + def _validate_aligned(self): + ref = self._open_gdal(self.ref_slc) + sec = self._open_gdal(self.sec_slc) + + a = (ref.RasterXSize, ref.RasterYSize) + b = (sec.RasterXSize, sec.RasterYSize) + + if a != b: + raise ValueError(f"SLC dimensions differ: REF={a}, SEC={b}") + + ref = None + sec = None + + return a + + @staticmethod + def _open_gdal(path): + ds = gdal.Open(str(path), gdal.GA_ReadOnly) + if ds is None: + raise RuntimeError(f"GDAL cannot open {path}") + return ds + + @classmethod + def read_array(cls, path, dtype=None): + ds = cls._open_gdal(path) + arr = ds.GetRasterBand(1).ReadAsArray() + ds = None + + if dtype is not None: + arr = arr.astype(dtype, copy=False) + + return arr + + def write_array(self, path, array, gdal_type): + path = Path(path) + self._prepare_output(path) + + arr = np.asarray(array) + length, width = arr.shape + + ds = gdal.GetDriverByName("ENVI").Create( + str(path), + width, + length, + 1, + gdal_type, + ) + ds.GetRasterBand(1).WriteArray(arr) + ds.FlushCache() + ds = None + + def _prepare_output(self, path: Path): + path.parent.mkdir(parents=True, exist_ok=True) + + if path.exists(): + if not self.config.overwrite: + raise FileExistsError( + f"{path} exists. Set overwrite=True." + ) + self.remove_envi(path) + + @staticmethod + def remove_envi(path): + path = Path(path) + for target in ( + path, + Path(str(path) + ".hdr"), + Path(str(path) + ".aux.xml"), + ): + if target.exists(): + target.unlink() + + +__all__ = [ + "IgramConfig", + "IgramProducts", + "PairIgramWorkflow", + "RadarMeta", +] diff --git a/share/stack/StackProcessor.py b/share/stack/StackProcessor.py index 9d8d7837b2..86c8c3da82 100644 --- a/share/stack/StackProcessor.py +++ b/share/stack/StackProcessor.py @@ -1,139 +1,331 @@ -from nisar.workflows.rubbersheet import run_rubbersheet_with_interpolation +from itertools import combinations import os from pathlib import Path -import shutil import pickle +import shutil + +from nisar.workflows.rubbersheet import run_rubbersheet_with_interpolation +from nisar.workflows.helpers import copy_raster from coarse_alignment import align_secondary, extract_ref_topo from invert_offsets import invert_offsets_isce3 -from nisar.workflows.helpers import copy_raster -from offset_util import (create_minimal_rifg_for_rubbersheet, - dense_offset_coregistered) +from offset_util import ( + create_minimal_rifg_for_rubbersheet, + dense_offset_coregistered, +) from ScanList import ScanList -from stack_utils import resample, relative_symlink_contents, set_nested, load_config +from stack_utils import ( + baseline, + load_config, + relative_symlink_contents, + resample, + set_nested, +) + -pol_freq = {'A': ['HH']} +pol_freq = {"A": ["HH"]} runcfg = load_config() cfg = runcfg.cfg["runconfig"]["groups"] dense_cfg = cfg["processing"]["dense_offsets"] -rdr2geo_cfg = cfg['processing']['rdr2geo'] -rubbersheet_cfg = cfg['processing']['rubbersheet'] +rdr2geo_cfg = cfg["processing"]["rdr2geo"] +rubbersheet_cfg = cfg["processing"]["rubbersheet"] class StackProcessor(ScanList): def __init__( - self, - slc_folder, - out_folder, - dem_file, - pair_level=2, - freq='A', - pol='HH', - overwrite=False, - debug=False): + self, + slc_folder, + out_folder, + dem_file, + max_baseline_distance=500, + max_date_apart=30, + freq="A", + pol="HH", + overwrite=False, + debug=False, + ): self.slc_folder = Path(slc_folder) self.out_folder = Path(out_folder) - self.config_save = self.out_folder / 'processor.pkl' self.dem_file = Path(dem_file) self.overwrite = overwrite self.freq = freq self.pol = pol self.debug = debug - h5_files = [str(i) for i in self.slc_folder.glob('*1.h5')] + + self.config_save = self.get_processor_config_path() + + h5_files = [str(i) for i in self.slc_folder.glob("*1.h5")] super().__init__(h5_files) + self.pol_freq = pol_freq - self.pair_level = pair_level - self.init_pairs() - if not os.path.exists(self.config_save) or self.overwrite: - pickle.dump(self, open(self.config_save, 'wb')) + self.max_baseline_distance = max_baseline_distance + self.max_date_apart = max_date_apart + + self.init_pairs_baseline_date() + + if not self.config_save.exists() or self.overwrite: + self.config_save.parent.mkdir(parents=True, exist_ok=True) + with self.config_save.open("wb") as f: + pickle.dump(self, f) + + # ------------------------------------------------------------------ + # Standardized path getters + # ------------------------------------------------------------------ + + @staticmethod + def _date_str(scan_or_date): + """Return YYYYMMDD string from a Scan-like object or a string.""" + return getattr(scan_or_date, "date_str", str(scan_or_date)) + + def get_processor_config_path(self): + return self.out_folder / "processor.pkl" + + def get_path_folder(self, path): + return self.out_folder / Path(path) + + def get_stack_folder(self, path): + return self.get_path_folder(path) / "stack" + + def get_ref_geom_folder(self, path): + return self.get_path_folder(path) / "ref_geom" + + def get_coarse_offsets_folder(self, path): + return self.get_path_folder(path) / "coarse_offsets" + + def get_coarse_offset_folder(self, path, date): + return self.get_coarse_offsets_folder(path) / self._date_str(date) + + def get_geo2rdr_folder(self, path, date): + return ( + self.get_coarse_offset_folder(path, date) + / "geo2rdr" + / f"freq{self.freq}" + ) + + def get_coarse_offset_file(self, path, date, component): + if component not in {"azimuth", "range"}: + raise ValueError( + "component must be either 'azimuth' or 'range'" + ) + return self.get_geo2rdr_folder(path, date) / f"{component}.off" + + def get_stack_slc(self, path, date): + date_str = self._date_str(date) + return self.get_stack_folder(path) / f"{date_str}.slc" + + def get_pair_name(self, ref, sec): + return f"{self._date_str(ref)}_{self._date_str(sec)}" + + def get_dense_pairwise_folder(self, path): + return self.get_path_folder(path) / "dense_pairwise" + + def get_dense_pair_folder(self, path, ref, sec): + return ( + self.get_dense_pairwise_folder(path) + / self.get_pair_name(ref, sec) + ) + + def get_dense_offsets_file(self, path, ref, sec): + return self.get_dense_pair_folder(path, ref, sec) / "dense_offsets" + + def get_inverted_offsets_folder(self, path): + return self.get_path_folder(path) / "inverted_offsets" + + def get_inverted_offset_folder(self, path, date): + return self.get_inverted_offsets_folder(path) / self._date_str(date) + + def get_rubber_folder(self, path): + return self.get_path_folder(path) / "rubber" + + def get_rubber_date_folder(self, path, date): + return self.get_rubber_folder(path) / self._date_str(date) - def init_pairs(self): + def get_rubber_pol_folder(self, path, date): + return ( + self.get_rubber_date_folder(path, date) + / f"freq{self.freq}" + / self.pol + ) + + def get_h5_stack_folder(self, path): + return self.get_path_folder(path) / "h5_stack" + + def get_h5_stack_date_folder(self, path, date): + return self.get_h5_stack_folder(path) / self._date_str(date) + + def get_dense_offsets_link_folder(self, path): + return ( + self.get_path_folder(path) + / "dense_offsets" + / f"freq{self.freq}" + / self.pol + ) + + def get_rubbersheet_offsets_folder(self, path): + return self.get_path_folder(path) / "rubbersheet_offsets" + + def get_merged_folder(self, path): + return self.get_path_folder(path) / "merged" + + def get_merged_date_folder(self, path, date): + return self.get_merged_folder(path) / self._date_str(date) + + def get_merged_slc(self, path, date): + date_str = self._date_str(date) + return self.get_merged_date_folder(path, date) / f"{date_str}.slc" + + def get_aligned_slc(self, path, date): + """Alias for the final aligned SLC used by downstream workflows.""" + return self.get_merged_slc(path, date) + + # ------------------------------------------------------------------ + # Pair selection + # ------------------------------------------------------------------ + + def init_pairs_baseline_date(self): self.dense_pairs = {} + for path, scans in self.scan_by_id.items(): - nscans = len(scans) path_pairs = [] - for i in range(self.pair_level): - if nscans - i > 1: - for j in range(nscans - i - 1): - path_pairs.append((scans[j], scans[j + i + 1])) + + for scan_a, scan_b in combinations(scans, 2): + b = baseline( + scan_a.path_str, + scan_b.path_str, + freq=self.freq, + ) + days_apart = abs((scan_a.date - scan_b.date).days) + + if self.debug: + print( + f"{self.get_pair_name(scan_a, scan_b)} " + f"{b['B']:12.2f} " + f"{b['Bperp']:12.2f} " + f"{b['Bpar']:14.2f} " + f"{b['Balong']:12.2f}" + f"{days_apart:12.2f}" + ) + + baseline_check = b["B"] <= self.max_baseline_distance + date_check = days_apart <= self.max_date_apart + + if baseline_check and date_check: + path_pairs.append((scan_a, scan_b)) + self.dense_pairs[path] = path_pairs def iter_path_dense_pairs(self, path): yield from self.dense_pairs[path] + # ------------------------------------------------------------------ + # Coarse registration + # ------------------------------------------------------------------ + def coarse_register_scans(self, rdr2geo_cfg=rdr2geo_cfg): for path in self.scan_by_id: main_ref = self.ref_scans[path] - path_folder = self.out_folder / Path(path) - stack_folder = path_folder / Path('stack') - ref_geom_folder = path_folder / Path('ref_geom') - coarse_offset_folder = path_folder / Path('coarse_offsets') - os.makedirs(path_folder, exist_ok=True) - os.makedirs(stack_folder, exist_ok=True) - os.makedirs(coarse_offset_folder, exist_ok=True) - reference_slc = stack_folder / Path(f'{main_ref.date_str}.slc') - if not os.path.exists(ref_geom_folder) or self.overwrite: - os.makedirs(ref_geom_folder, exist_ok=True) + + path_folder = self.get_path_folder(path) + stack_folder = self.get_stack_folder(path) + ref_geom_folder = self.get_ref_geom_folder(path) + coarse_offsets_folder = self.get_coarse_offsets_folder(path) + + path_folder.mkdir(parents=True, exist_ok=True) + stack_folder.mkdir(parents=True, exist_ok=True) + coarse_offsets_folder.mkdir(parents=True, exist_ok=True) + + reference_slc = self.get_stack_slc(path, main_ref) + + if not ref_geom_folder.exists() or self.overwrite: + ref_geom_folder.mkdir(parents=True, exist_ok=True) extract_ref_topo( main_ref.path_str, self.dem_file, rdr2geo_cfg, - ref_geom_folder) - if not os.path.exists(reference_slc) or self.overwrite: + ref_geom_folder, + ) + + if not reference_slc.exists() or self.overwrite: copy_raster( main_ref.path_str, - self.freq, # frequency - self.pol, # polarization - 1024, # lines per block + self.freq, + self.pol, + 1024, reference_slc, file_type="ENVI", ) + pairs = {} + for ref, sec in self.iter_path_ref_sec(path): - coarse_offset_folderi = coarse_offset_folder / \ - Path(sec.date_str) - val = {} - val['azimuth'] = coarse_offset_folderi / \ - Path(f'geo2rdr/freq{self.freq}/azimuth.off') - val['range'] = coarse_offset_folderi / \ - Path(f'geo2rdr/freq{self.freq}/range.off') - pairs[(ref.date_str, sec.date_str)] = val - if not os.path.exists(coarse_offset_folderi) or self.overwrite: - os.makedirs(coarse_offset_folderi, exist_ok=True) + coarse_offset_folder = self.get_coarse_offset_folder( + path, + sec, + ) + + pairs[(ref.date_str, sec.date_str)] = { + "azimuth": self.get_coarse_offset_file( + path, + sec, + "azimuth", + ), + "range": self.get_coarse_offset_file( + path, + sec, + "range", + ), + } + + if not coarse_offset_folder.exists() or self.overwrite: + coarse_offset_folder.mkdir( + parents=True, + exist_ok=True, + ) align_secondary( sec.path_str, self.dem_file, ref_geom_folder, - coarse_offset_folderi) + coarse_offset_folder, + ) + + sec_slc = self.get_stack_slc(path, sec) + + if not sec_slc.exists() or self.overwrite: + print( + f"calculating coarse reg for " + f"{path} {sec.date_str}" + ) - sec_slc = stack_folder / Path(f'{sec.date_str}.slc') - if not os.path.exists(sec_slc) or self.overwrite: - print(f'calculating coarse reg for {path} {sec.date_str}') - offset_dir = coarse_offset_folderi / "geo2rdr" / "freqA" resample( stack_folder, - offset_dir, + self.get_geo2rdr_folder(path, sec), main_ref.path_str, sec.path_str, - sec.date_str) + sec.date_str, + ) + + # ------------------------------------------------------------------ + # Dense offsets + # ------------------------------------------------------------------ - def dense_offset_pairs(self, dense_cfg): + def dense_offset_pairs(self, dense_cfg=dense_cfg): for path in self.scan_by_id: + dense_folder = self.get_dense_pairwise_folder(path) + dense_folder.mkdir(parents=True, exist_ok=True) + for ref, sec in self.iter_path_dense_pairs(path): - path_folder = self.out_folder / Path(path) - dense_folder = path_folder / Path('dense_pairwise') - stack_folder = path_folder / Path('stack') - os.makedirs(dense_folder, exist_ok=True) - starting_pair = f'{ref.date_str}_{sec.date_str}' - out_dir = dense_folder / Path(starting_pair) - if not os.path.exists(out_dir) or self.overwrite: - print(f'calculating dense offset {path} {starting_pair}') - os.makedirs(out_dir, exist_ok=True) - reference_slc = stack_folder / Path(f'{ref.date_str}.slc') - secondary_slc = stack_folder / Path(f"{sec.date_str}.slc") + pair_name = self.get_pair_name(ref, sec) + out_dir = self.get_dense_pair_folder(path, ref, sec) + + if not out_dir.exists() or self.overwrite: + print( + f"calculating dense offset {path} {pair_name}" + ) + out_dir.mkdir(parents=True, exist_ok=True) + dense_offset_coregistered( - reference_slc, - secondary_slc, + self.get_stack_slc(path, ref), + self.get_stack_slc(path, sec), out_dir, dense_cfg=dense_cfg, gpu_id=0, @@ -141,103 +333,195 @@ def dense_offset_pairs(self, dense_cfg): def invert_dense_offset_pairs(self): for path in self.scan_by_id: - pairs = {} - path_folder = self.out_folder / Path(path) - dense_folder = path_folder / Path('dense_pairwise') - inverted_offsets_folder = path_folder / Path('inverted_offsets') - if not os.path.exists(inverted_offsets_folder) or self.debug: - os.makedirs(inverted_offsets_folder, exist_ok=True) + inverted_offsets_folder = self.get_inverted_offsets_folder(path) + + # Preserve original behavior: debug=True forces this block. + if not inverted_offsets_folder.exists() or self.debug: + inverted_offsets_folder.mkdir( + parents=True, + exist_ok=True, + ) + + pairs = {} + for ref, sec in self.iter_path_dense_pairs(path): - starting_pair = f'{ref.date_str}_{sec.date_str}' - out_dir = dense_folder / Path(starting_pair) - pairs[(ref.date_str, sec.date_str) - ] = out_dir / 'dense_offsets' - if len(pairs) > 0: - outputs = invert_offsets_isce3( + pairs[(ref.date_str, sec.date_str)] = ( + self.get_dense_offsets_file(path, ref, sec) + ) + + if pairs: + invert_offsets_isce3( pair_offsets=pairs, output_dir=str(inverted_offsets_folder), reference_date=self.ref_scans[path].date_str, ) + # ------------------------------------------------------------------ + # Rubbersheet + # ------------------------------------------------------------------ + def rubbersheet(self, rubbersheet_cfg=rubbersheet_cfg): for path in self.scan_by_id: - path_folder = self.out_folder / Path(path) - coarse_offset_folder = path_folder / Path('coarse_offsets') - inverted_offsets_folder = path_folder / Path('inverted_offsets') - dense_folder = path_folder / Path('dense_pairwise') - rubber_folder = path_folder / Path('rubber') - h5_stack_folder = path_folder / Path('h5_stack') - os.makedirs(h5_stack_folder, exist_ok=True) - os.makedirs(rubber_folder, exist_ok=True) + path_folder = self.get_path_folder(path) + h5_stack_folder = self.get_h5_stack_folder(path) + rubber_folder = self.get_rubber_folder(path) + + h5_stack_folder.mkdir(parents=True, exist_ok=True) + rubber_folder.mkdir(parents=True, exist_ok=True) + main_ref = self.ref_scans[path] + for ref, sec in self.iter_path_ref_sec(path): - coarse_offset_folderi = coarse_offset_folder / \ - Path(sec.date_str) - out_dir = h5_stack_folder / sec.date_str - rubberi = rubber_folder / sec.date_str - if not os.path.exists(rubberi) or self.overwrite: - print(f'calculating rubbersheet {path} {sec.date_str}') - cfg['processing']['rubbersheet'] = rubbersheet_cfg - inverted_offset_pathi = inverted_offsets_folder / sec.date_str - starting_pair = f'{ref.date_str}_{sec.date_str}' - densei = dense_folder / Path(starting_pair) - link_to = path_folder / Path('dense_offsets/freqA/HH') - relative_symlink_contents(densei, link_to) - relative_symlink_contents(inverted_offset_pathi, link_to) - cfg['dynamic_ancillary_file_group']['dem_file'] = str( - self.dem_file) + coarse_offset_folder = self.get_coarse_offset_folder( + path, + sec, + ) + out_dir = self.get_h5_stack_date_folder(path, sec) + rubber_date_folder = self.get_rubber_date_folder( + path, + sec, + ) + + if not rubber_date_folder.exists() or self.overwrite: + print( + f"calculating rubbersheet " + f"{path} {sec.date_str}" + ) + + cfg["processing"]["rubbersheet"] = rubbersheet_cfg + + inverted_offset_folder = ( + self.get_inverted_offset_folder(path, sec) + ) + dense_pair_folder = self.get_dense_pair_folder( + path, + ref, + sec, + ) + link_to = self.get_dense_offsets_link_folder(path) + + relative_symlink_contents( + dense_pair_folder, + link_to, + ) + relative_symlink_contents( + inverted_offset_folder, + link_to, + ) + + cfg["dynamic_ancillary_file_group"]["dem_file"] = str( + self.dem_file + ) + set_nested( - cfg, [ - 'input_file_group', 'reference_rslc_file'], main_ref.path_str) + cfg, + [ + "input_file_group", + "reference_rslc_file", + ], + main_ref.path_str, + ) set_nested( - cfg, [ - 'product_path_group', 'sas_output_file'], out_dir) + cfg, + [ + "product_path_group", + "sas_output_file", + ], + out_dir, + ) set_nested( - cfg, [ - 'product_path_group', 'scratch_path'], path_folder) - cfg['processing']['rubbersheet']['geo2rdr_offsets_path'] = coarse_offset_folderi - cfg['processing']['rubbersheet']['dense_offsets_path'] = path_folder - cfg['processing']['input_subset']['list_of_frequencies'] = { - self.freq: [self.pol]} + cfg, + [ + "product_path_group", + "scratch_path", + ], + path_folder, + ) + + cfg["processing"]["rubbersheet"][ + "geo2rdr_offsets_path" + ] = coarse_offset_folder + + cfg["processing"]["rubbersheet"][ + "dense_offsets_path" + ] = path_folder + + cfg["processing"]["input_subset"][ + "list_of_frequencies" + ] = { + self.freq: [self.pol] + } + create_minimal_rifg_for_rubbersheet( - out_dir, cfg) + out_dir, + cfg, + ) run_rubbersheet_with_interpolation( - cfg, out_dir) + cfg, + out_dir, + ) + shutil.move( - path_folder / - Path('rubbersheet_offsets'), - rubberi) + self.get_rubbersheet_offsets_folder(path), + rubber_date_folder, + ) + + # ------------------------------------------------------------------ + # Final aligned SLC stack + # ------------------------------------------------------------------ def resample_slcs(self): for path in self.scan_by_id: - path_folder = self.out_folder / Path(path) - merged_folder = path_folder / Path('merged') - rubber_folder = path_folder / Path('rubber') - os.makedirs(merged_folder, exist_ok=True) + merged_folder = self.get_merged_folder(path) + merged_folder.mkdir(parents=True, exist_ok=True) + main_ref = self.ref_scans[path] - ref_slc_folder = merged_folder / main_ref.date_str - reference_slc = ref_slc_folder / Path(f'{main_ref.date_str}.slc') - if not os.path.exists(reference_slc) or self.overwrite: - print(f'copying reference slc for {path}') - os.makedirs(ref_slc_folder, exist_ok=True) + ref_slc_folder = self.get_merged_date_folder( + path, + main_ref, + ) + reference_slc = self.get_merged_slc( + path, + main_ref, + ) + + if not reference_slc.exists() or self.overwrite: + print(f"copying reference slc for {path}") + ref_slc_folder.mkdir(parents=True, exist_ok=True) + copy_raster( main_ref.path_str, - self.freq, # frequency - self.pol, # polarization - 1024, # lines per block + self.freq, + self.pol, + 1024, reference_slc, file_type="ENVI", ) + for _, sec in self.iter_path_ref_sec(path): - print(f'calculating final offset for {path} {sec.date_str}') - rubberi = rubber_folder / sec.date_str / \ - f'freq{self.freq}' / self.pol - pairs_out = merged_folder / sec.date_str - if not os.path.exists(pairs_out) or self.overwrite: - print(f'resampling {path} {sec.date_str}') + print( + f"calculating final offset for " + f"{path} {sec.date_str}" + ) + + rubber_pol_folder = self.get_rubber_pol_folder( + path, + sec, + ) + pairs_out = self.get_merged_date_folder( + path, + sec, + ) + + if not pairs_out.exists() or self.overwrite: + print( + f"resampling {path} {sec.date_str}" + ) + resample( pairs_out, - rubberi, + rubber_pol_folder, main_ref.path_str, sec.path_str, - out_tag=sec.date_str) + out_tag=sec.date_str, + ) diff --git a/share/stack/iono.py b/share/stack/iono.py new file mode 100644 index 0000000000..eb31c00e2e --- /dev/null +++ b/share/stack/iono.py @@ -0,0 +1,1358 @@ +""" +isce3_pair_iono +=============== + +Importable pair-level ISCE3 interferogram + split-spectrum ionosphere workflow +for two already-aligned SLC rasters. + +Typical use +----------- +from isce3_pair_iono import ( + AlignedPairIonoWorkflow, + WorkflowConfig, +) + +cfg = WorkflowConfig( + range_looks=11, + azimuth_looks=11, + gpu=True, +) + +wf = AlignedPairIonoWorkflow( + ref_slc="/path/reference.slc", + sec_slc="/path/coregistered_secondary.slc", + outdir="/path/pair_iono", + metadata_rslc="/path/reference_rslc.h5", + ref_rslc="/path/reference_rslc.h5", + sec_rslc="/path/secondary_rslc.h5", + config=cfg, +) + +products = wf.run() + +Notes +----- +The official full NISAR ISCE3 ionosphere workflow splits the ORIGINAL RSLCs +before resampling/coregistering the low/high secondary subbands. + +This module instead starts from two already-aligned SLC rasters, so the split +spectrum operation necessarily happens after alignment. From Crossmul onward, +it uses the same ISCE3 module families as the current NISAR workflow: + +- isce3.signal.Crossmul +- isce3.splitspectrum.SplitSpectrum +- SNAPHU (default), isce3.unwrap.ICU, or Whirlwind +- LowHighSubbandIonosphereEstimation +- IonosphereFilter +- split-band unwrapping-error estimation/correction +""" + +from __future__ import annotations + +import json +import shutil +from dataclasses import asdict, dataclass +from pathlib import Path +from typing import Callable, Optional + +import numpy as np +from osgeo import gdal +from scipy.fft import next_fast_len + +import isce3 +from isce3.atmosphere.ionosphere_filter import IonosphereFilter +from isce3.atmosphere.split_band_estimation import ( + LowHighSubbandIonosphereEstimation, +) +from isce3.splitspectrum.splitspectrum import BandpassMetaData, SplitSpectrum +from isce3.unwrap.bridge_phase import bridge_unwrapped_phase +from nisar.products.readers import RSLC + +gdal.UseExceptions() + +C = float(isce3.core.speed_of_light) + + +@dataclass +class WorkflowConfig: + """Processing configuration for :class:`AlignedPairIonoWorkflow`.""" + + # Crossmul + range_looks: int = 11 + azimuth_looks: int = 11 + oversample: int = 2 + crossmul_lines_per_block: int = 1024 + gpu: bool = False + gpu_id: int = 0 + + # Optional flattening + flatten_range_offset: Optional[str] = None + starting_range_shift: float = 0.0 + + # Split spectrum + split_lines_per_block: int = 2048 + subband_fraction: float = 1.0 / 3.0 + window_function: str = "tukey" + window_shape: float = 0.25 + + # Unwrapping: "snaphu", "icu", or "whirlwind" + unwrap_algorithm: str = "snaphu" + + # Whirlwind (package: whirlwind-insar; import: whirlwind) + whirlwind_nlooks: Optional[float] = None + whirlwind_downsample: Optional[int] = None + + # SNAPHU defaults close to current NISAR ISCE3 defaults + snaphu_nlooks: Optional[float] = None + snaphu_cost: str = "smooth" + snaphu_init: str = "mcf" + snaphu_min_conncomp_frac: float = 0.01 + snaphu_phase_grad_window: tuple[int, int] = (7, 7) + snaphu_ntiles: tuple[int, int] = (1, 1) + snaphu_tile_overlap: tuple[int, int] = (0, 0) + snaphu_nproc: int = 1 + snaphu_tile_cost_thresh: float = 500.0 + snaphu_min_region_size: int = 300 + snaphu_single_tile_reoptimize: bool = True + + # Bridge disconnected unwrapped regions + bridge: bool = True + bridge_radius: int = 500 + bridge_erosion_size: int = 2 + bridge_minimum_samples: int = 14 + bridge_ramp_type: Optional[str] = None + bridge_ramp_maximum_pixel: int = 1_000_000 + + # Ionosphere mask + iono_mask_types: tuple[str, ...] = ("coherence",) + iono_coherence_threshold: float = 0.5 + + # Ionosphere filtering + iono_filter: bool = True + iono_kernel_range: int = 100 + iono_kernel_azimuth: int = 100 + iono_sigma_range: int = 33 + iono_sigma_azimuth: int = 33 + iono_filter_iterations: int = 1 + iono_filling_method: str = "nearest" + iono_min_cluster_pixels: int = 2 + + # Current official-style split-band unwrap-error correction + iono_unwrap_correction: bool = True + + # File handling + overwrite: bool = False + + +@dataclass +class RadarMeta: + """Radar metadata needed by split-spectrum and Crossmul.""" + + center_frequency: float + range_bandwidth: float + range_sample_frequency: float + range_pixel_spacing: float + starting_range: float + wavelength: float + + @classmethod + def from_nisar_rslc(cls, rslc_path: str | Path, frequency: str = "A"): + slc = RSLC(hdf5file=str(rslc_path)) + meta = BandpassMetaData.load_from_slc( + slc_product=slc, + freq=frequency, + ) + grid = slc.getRadarGrid(frequency) + + return cls( + center_frequency=float(meta.center_freq), + range_bandwidth=float(meta.rg_bandwidth), + range_sample_frequency=float(meta.rg_sample_freq), + range_pixel_spacing=float(meta.rg_pxl_spacing), + starting_range=float(grid.starting_range), + wavelength=float(meta.wavelength), + ) + + @classmethod + def from_values( + cls, + *, + center_frequency: float, + range_bandwidth: float, + range_sample_frequency: float, + starting_range: float = 0.0, + ): + dr = C / (2.0 * float(range_sample_frequency)) + f0 = float(center_frequency) + + return cls( + center_frequency=f0, + range_bandwidth=float(range_bandwidth), + range_sample_frequency=float(range_sample_frequency), + range_pixel_spacing=dr, + starting_range=float(starting_range), + wavelength=C / f0, + ) + + def slant_range(self, index: int) -> float: + return self.starting_range + float(index) * self.range_pixel_spacing + + +@dataclass +class WorkflowProducts: + """Paths to the principal products generated by the workflow.""" + + main_wrapped_interferogram: Path + main_coherence: Path + main_unwrapped_phase: Path + main_connected_components: Path + + low_wrapped_interferogram: Path + low_coherence: Path + low_unwrapped_phase: Path + low_connected_components: Path + + high_wrapped_interferogram: Path + high_coherence: Path + high_unwrapped_phase: Path + high_connected_components: Path + + ionosphere_phase_screen: Path + ionosphere_uncertainty: Path + ionosphere_mask: Path + nondispersive_phase_screen: Path + + main_iono_corrected_unwrapped_phase: Optional[Path] = None + + +class AlignedPairIonoWorkflow: + """ + Pair-level ISCE3 workflow for two already-aligned SLC rasters. + + Parameters + ---------- + ref_slc, sec_slc + GDAL-readable complex SLC rasters on the same radar grid. + outdir + Output directory. + metadata_rslc + Optional NISAR RSLC HDF5 used to derive frequency/range metadata. + radar_meta + Explicit :class:`RadarMeta`. Use this instead of ``metadata_rslc``. + ref_rslc, sec_rslc + Optional original NISAR RSLC HDF5 files. If supplied, their Doppler + LUTs are passed to Crossmul like the official NISAR workflow. + frequency + NISAR frequency, normally ``"A"``. + config + :class:`WorkflowConfig`. + """ + + def __init__( + self, + *, + ref_slc: str | Path, + sec_slc: str | Path, + outdir: str | Path, + metadata_rslc: str | Path | None = None, + radar_meta: RadarMeta | None = None, + ref_rslc: str | Path | None = None, + sec_rslc: str | Path | None = None, + frequency: str = "A", + config: WorkflowConfig | None = None, + ): + self.ref_slc = Path(ref_slc).resolve() + self.sec_slc = Path(sec_slc).resolve() + self.outdir = Path(outdir).resolve() + + self.frequency = frequency + self.config = config or WorkflowConfig() + + self.ref_rslc = Path(ref_rslc).resolve() if ref_rslc else None + self.sec_rslc = Path(sec_rslc).resolve() if sec_rslc else None + + if radar_meta is not None and metadata_rslc is not None: + raise ValueError( + "Provide either radar_meta or metadata_rslc, not both." + ) + + if radar_meta is not None: + self.meta = radar_meta + elif metadata_rslc is not None: + self.meta = RadarMeta.from_nisar_rslc( + metadata_rslc, + frequency=frequency, + ) + else: + raise ValueError( + "Provide either metadata_rslc or an explicit RadarMeta." + ) + + self.outdir.mkdir(parents=True, exist_ok=True) + + self.width, self.length = self._validate_aligned() + self.dopplers = self._load_dopplers() + + ( + self.low_limits, + self.high_limits, + self.low_center_frequency, + self.high_center_frequency, + ) = self._compute_subbands() + + self._write_workflow_metadata() + + # ------------------------------------------------------------------ + # Public high-level API + # ------------------------------------------------------------------ + + def run(self) -> WorkflowProducts: + """Run the complete pair workflow.""" + + print("[1/6] Main-band Crossmul") + main_ifg, main_coh = self.crossmul_main() + + print("[2/6] Main-band unwrap") + main_unw, main_cc = self.unwrap_main() + + print("[3/6] Split spectrum") + self.split_aligned_slcs() + + print("[4/6] Low/high Crossmul") + low_ifg, low_coh, high_ifg, high_coh = self.crossmul_subbands() + + print("[5/6] Low/high unwrap") + low_unw, low_cc, high_unw, high_cc = self.unwrap_subbands() + + print("[6/6] Ionosphere estimation") + iono, iono_sigma, iono_mask, nondisp = self.estimate_ionosphere() + + corrected = self.correct_main_unwrapped_phase() + + return WorkflowProducts( + main_wrapped_interferogram=main_ifg, + main_coherence=main_coh, + main_unwrapped_phase=main_unw, + main_connected_components=main_cc, + low_wrapped_interferogram=low_ifg, + low_coherence=low_coh, + low_unwrapped_phase=low_unw, + low_connected_components=low_cc, + high_wrapped_interferogram=high_ifg, + high_coherence=high_coh, + high_unwrapped_phase=high_unw, + high_connected_components=high_cc, + ionosphere_phase_screen=iono, + ionosphere_uncertainty=iono_sigma, + ionosphere_mask=iono_mask, + nondispersive_phase_screen=nondisp, + main_iono_corrected_unwrapped_phase=corrected, + ) + + def crossmul_main(self) -> tuple[Path, Path]: + """Form the main-band wrapped interferogram and coherence.""" + return self._crossmul_pair( + self.ref_slc, + self.sec_slc, + self.outdir / "main", + ) + + def unwrap_main(self) -> tuple[Path, Path]: + """Unwrap the main-band interferogram.""" + main_dir = self.outdir / "main" + return self._unwrap_pair( + main_dir / "wrappedInterferogram.int", + main_dir / "coherenceMagnitude.cor", + main_dir, + ) + + def split_aligned_slcs(self) -> dict[str, Path]: + """Split the two aligned SLCs into low/high range subbands.""" + + split_dir = self.outdir / "split_spectrum" + split_dir.mkdir(parents=True, exist_ok=True) + + paths = { + "ref_low": split_dir / "ref_low.slc", + "ref_high": split_dir / "ref_high.slc", + "sec_low": split_dir / "sec_low.slc", + "sec_high": split_dir / "sec_high.slc", + } + + self._split_slc( + self.ref_slc, + paths["ref_low"], + paths["ref_high"], + ) + self._split_slc( + self.sec_slc, + paths["sec_low"], + paths["sec_high"], + ) + + return paths + + def crossmul_subbands( + self, + ) -> tuple[Path, Path, Path, Path]: + """Form low/high subband interferograms and coherences.""" + + split_dir = self.outdir / "split_spectrum" + + low_ifg, low_coh = self._crossmul_pair( + split_dir / "ref_low.slc", + split_dir / "sec_low.slc", + self.outdir / "low", + ) + + high_ifg, high_coh = self._crossmul_pair( + split_dir / "ref_high.slc", + split_dir / "sec_high.slc", + self.outdir / "high", + ) + + return low_ifg, low_coh, high_ifg, high_coh + + def unwrap_subbands( + self, + ) -> tuple[Path, Path, Path, Path]: + """Unwrap low/high subband interferograms.""" + + low_dir = self.outdir / "low" + high_dir = self.outdir / "high" + + low_unw, low_cc = self._unwrap_pair( + low_dir / "wrappedInterferogram.int", + low_dir / "coherenceMagnitude.cor", + low_dir, + ) + + high_unw, high_cc = self._unwrap_pair( + high_dir / "wrappedInterferogram.int", + high_dir / "coherenceMagnitude.cor", + high_dir, + ) + + return low_unw, low_cc, high_unw, high_cc + + def estimate_ionosphere( + self, + ) -> tuple[Path, Path, Path, Path]: + """ + Estimate, mask, optionally unwrap-correct, and filter the ionosphere. + """ + + cfg = self.config + low_dir = self.outdir / "low" + high_dir = self.outdir / "high" + iono_dir = self.outdir / "ionosphere" + iono_dir.mkdir(parents=True, exist_ok=True) + + phi_low = self._read_array( + low_dir / "unwrappedPhase.unw", + np.float64, + ) + phi_high = self._read_array( + high_dir / "unwrappedPhase.unw", + np.float64, + ) + + low_coh = self._read_array( + low_dir / "coherenceMagnitude.cor", + np.float32, + ) + high_coh = self._read_array( + high_dir / "coherenceMagnitude.cor", + np.float32, + ) + + low_cc = self._read_array( + low_dir / "connectedComponents.conncomp" + ) + high_cc = self._read_array( + high_dir / "connectedComponents.conncomp" + ) + + if cfg.bridge: + phi_low = self._bridge_phase(phi_low) + phi_high = self._bridge_phase(phi_high) + + estimator = LowHighSubbandIonosphereEstimation( + main_center_freq=self.meta.center_frequency, + low_center_freq=self.low_center_frequency, + high_center_freq=self.high_center_frequency, + ) + + dispersive, nondispersive = estimator.compute_disp_nondisp( + phi_sub_low=phi_low, + phi_sub_high=phi_high, + no_data=0, + ) + + nlooks = cfg.range_looks * cfg.azimuth_looks + + iono_sigma, nondisp_sigma = estimator.estimate_iono_std( + low_band_coh=low_coh, + high_band_coh=high_coh, + number_looks=nlooks, + resample_flag=False, + ) + + mask = self._make_iono_mask( + estimator, + phi_low, + phi_high, + low_coh, + high_coh, + low_cc, + high_cc, + ) + + raw_disp = iono_dir / "dispersive.raw" + raw_non = iono_dir / "nonDispersive.raw" + raw_sig = iono_dir / "dispersive.sig" + raw_non_sig = iono_dir / "nonDispersive.sig" + mask_path = iono_dir / "ionoValidMask.mask" + + self._write_array( + raw_disp, + dispersive.astype(np.float32), + gdal.GDT_Float32, + ) + self._write_array( + raw_non, + nondispersive.astype(np.float32), + gdal.GDT_Float32, + ) + self._write_array( + raw_sig, + iono_sigma.astype(np.float32), + gdal.GDT_Float32, + ) + self._write_array( + raw_non_sig, + nondisp_sigma.astype(np.float32), + gdal.GDT_Float32, + ) + self._write_array( + mask_path, + mask.astype(np.uint8), + gdal.GDT_Byte, + ) + + final_disp = iono_dir / "ionospherePhaseScreen" + final_sig = iono_dir / "ionospherePhaseScreenUncertainty" + final_non = iono_dir / "nonDispersivePhaseScreen" + final_non_sig = iono_dir / "nonDispersivePhaseScreenUncertainty" + + if not cfg.iono_filter: + self._write_array( + final_disp, + np.where(mask, dispersive, 0).astype(np.float32), + gdal.GDT_Float32, + ) + self._write_array( + final_sig, + np.where(mask, iono_sigma, 0).astype(np.float32), + gdal.GDT_Float32, + ) + self._write_array( + final_non, + np.where(mask, nondispersive, 0).astype(np.float32), + gdal.GDT_Float32, + ) + self._write_array( + final_non_sig, + np.where(mask, nondisp_sigma, 0).astype(np.float32), + gdal.GDT_Float32, + ) + return final_disp, final_sig, mask_path, final_non + + if cfg.iono_unwrap_correction: + prelim_disp = iono_dir / "dispersive.prelim_filt" + prelim_disp_sig = iono_dir / "dispersive.prelim_filt.sig" + prelim_non = iono_dir / "nonDispersive.prelim_filt" + prelim_non_sig = iono_dir / "nonDispersive.prelim_filt.sig" + + self._run_iono_filter( + raw_disp, + raw_sig, + mask_path, + prelim_disp, + prelim_disp_sig, + iono_dir / "filter_prelim_disp", + ) + self._run_iono_filter( + raw_non, + raw_non_sig, + mask_path, + prelim_non, + prelim_non_sig, + iono_dir / "filter_prelim_non", + ) + + filt_disp = self._read_array(prelim_disp, np.float64) + filt_non = self._read_array(prelim_non, np.float64) + + common_coef, diff_coef = estimator.compute_unwrapp_error( + disp_array=filt_disp, + nondisp_array=filt_non, + low_sub_runw=phi_low, + high_sub_runw=phi_high, + ) + + self._write_array( + iono_dir / "commonUnwrapErrorCoefficient", + common_coef.astype(np.int32), + gdal.GDT_Int32, + ) + self._write_array( + iono_dir / "differentialUnwrapErrorCoefficient", + diff_coef.astype(np.int32), + gdal.GDT_Int32, + ) + + corrected_disp, corrected_non = estimator.compute_disp_nondisp( + phi_sub_low=phi_low, + phi_sub_high=phi_high, + comm_unwcor_coef=common_coef, + diff_unwcor_coef=diff_coef, + no_data=0, + ) + + corr_disp = iono_dir / "dispersive.unwrap_corrected" + corr_non = iono_dir / "nonDispersive.unwrap_corrected" + + self._write_array( + corr_disp, + corrected_disp.astype(np.float32), + gdal.GDT_Float32, + ) + self._write_array( + corr_non, + corrected_non.astype(np.float32), + gdal.GDT_Float32, + ) + + self._run_iono_filter( + corr_disp, + raw_sig, + mask_path, + final_disp, + final_sig, + iono_dir / "filter_final_disp", + ) + self._run_iono_filter( + corr_non, + raw_non_sig, + mask_path, + final_non, + final_non_sig, + iono_dir / "filter_final_non", + ) + + else: + self._run_iono_filter( + raw_disp, + raw_sig, + mask_path, + final_disp, + final_sig, + iono_dir / "filter_final_disp", + ) + self._run_iono_filter( + raw_non, + raw_non_sig, + mask_path, + final_non, + final_non_sig, + iono_dir / "filter_final_non", + ) + + return final_disp, final_sig, mask_path, final_non + + def correct_main_unwrapped_phase(self) -> Optional[Path]: + """ + Subtract the estimated ionosphere screen from the main unwrapped phase. + """ + + main_path = self.outdir / "main" / "unwrappedPhase.unw" + iono_path = self.outdir / "ionosphere" / "ionospherePhaseScreen" + + if not main_path.exists() or not iono_path.exists(): + return None + + main = self._read_array(main_path, np.float32) + iono = self._read_array(iono_path, np.float32) + + if main.shape != iono.shape: + return None + + valid = np.isfinite(main) & np.isfinite(iono) + + corrected = np.zeros_like(main, dtype=np.float32) + corrected[valid] = main[valid] - iono[valid] + + out = ( + self.outdir + / "main" + / "unwrappedPhase_ionoCorrected.unw" + ) + + self._write_array( + out, + corrected, + gdal.GDT_Float32, + ) + + return out + + # ------------------------------------------------------------------ + # Internal utilities + # ------------------------------------------------------------------ + + def _validate_aligned(self) -> tuple[int, int]: + ref = self._open_gdal(self.ref_slc) + sec = self._open_gdal(self.sec_slc) + + ref_shape = (ref.RasterXSize, ref.RasterYSize) + sec_shape = (sec.RasterXSize, sec.RasterYSize) + + if ref_shape != sec_shape: + raise ValueError( + f"SLC dimensions differ: REF={ref_shape}, SEC={sec_shape}" + ) + + ref = None + sec = None + return ref_shape + + def _load_dopplers(self): + if self.ref_rslc is None or self.sec_rslc is None: + return None + + ref = RSLC(hdf5file=str(self.ref_rslc)) + sec = RSLC(hdf5file=str(self.sec_rslc)) + + ref_dopp = isce3.core.avg_lut2d_to_lut1d( + ref.getDopplerCentroid(frequency=self.frequency) + ) + sec_dopp = isce3.core.avg_lut2d_to_lut1d( + sec.getDopplerCentroid(frequency=self.frequency) + ) + + return ref_dopp, sec_dopp + + def _compute_subbands(self): + cfg = self.config + f0 = self.meta.center_frequency + bw = self.meta.range_bandwidth + sub_bw = bw * cfg.subband_fraction + + if not (0 < sub_bw <= bw / 2): + raise ValueError( + "subband_fraction produces an invalid subband bandwidth" + ) + + f_min = f0 - bw / 2 + f_max = f0 + bw / 2 + + low_limits = (f_min, f_min + sub_bw) + high_limits = (f_max - sub_bw, f_max) + + f_low = 0.5 * sum(low_limits) + f_high = 0.5 * sum(high_limits) + + return low_limits, high_limits, f_low, f_high + + def _write_workflow_metadata(self): + metadata = { + "ref_slc": str(self.ref_slc), + "sec_slc": str(self.sec_slc), + "width": self.width, + "length": self.length, + "frequency": self.frequency, + "radar_meta": asdict(self.meta), + "config": asdict(self.config), + "low_band_limits_hz": self.low_limits, + "high_band_limits_hz": self.high_limits, + "low_center_frequency_hz": self.low_center_frequency, + "high_center_frequency_hz": self.high_center_frequency, + "note": ( + "Pair-level aligned-SLC workflow. Official full NISAR ISCE3 " + "splits original RSLCs before subband resampling/coregistration." + ), + } + + (self.outdir / "workflow_metadata.json").write_text( + json.dumps(metadata, indent=2) + ) + + def _new_crossmul(self): + cfg = self.config + + if cfg.gpu: + device = isce3.cuda.core.Device(cfg.gpu_id) + isce3.cuda.core.set_device(device) + cm = isce3.cuda.signal.Crossmul() + else: + cm = isce3.signal.Crossmul() + + cm.range_looks = cfg.range_looks + cm.az_looks = cfg.azimuth_looks + cm.oversample_factor = cfg.oversample + cm.lines_per_block = cfg.crossmul_lines_per_block + + cm.range_pixel_spacing = self.meta.range_pixel_spacing + cm.wavelength = self.meta.wavelength + + if self.dopplers is not None: + cm.set_dopplers(*self.dopplers) + + if cfg.flatten_range_offset: + cm.ref_sec_offset_starting_range_shift = ( + cfg.starting_range_shift + ) + + return cm + + def _crossmul_pair( + self, + ref_path: Path, + sec_path: Path, + outdir: Path, + ) -> tuple[Path, Path]: + cfg = self.config + outdir.mkdir(parents=True, exist_ok=True) + + ifg_path = outdir / "wrappedInterferogram.int" + coh_path = outdir / "coherenceMagnitude.cor" + + self._prepare_output(ifg_path) + self._prepare_output(coh_path) + + ref = isce3.io.Raster(str(ref_path)) + sec = isce3.io.Raster(str(sec_path)) + + if ref.width != sec.width or ref.length != sec.length: + raise ValueError( + f"Grid mismatch: {ref_path} vs {sec_path}" + ) + + out_width = ref.width // cfg.range_looks + out_length = ref.length // cfg.azimuth_looks + + ifg = isce3.io.Raster( + str(ifg_path), + out_width, + out_length, + 1, + gdal.GDT_CFloat32, + "ENVI", + ) + coh = isce3.io.Raster( + str(coh_path), + out_width, + out_length, + 1, + gdal.GDT_Float32, + "ENVI", + ) + + cm = self._new_crossmul() + + flatten = ( + isce3.io.Raster(str(cfg.flatten_range_offset)) + if cfg.flatten_range_offset + else None + ) + + cm.crossmul( + ref, + sec, + ifg, + coh, + flatten, + ) + + del ifg + del coh + del ref + del sec + + if flatten is not None: + del flatten + + return ifg_path, coh_path + + def _split_slc( + self, + src_path: Path, + low_path: Path, + high_path: Path, + ): + cfg = self.config + + self._prepare_output(low_path) + self._prepare_output(high_path) + + src = self._open_gdal(src_path) + width = src.RasterXSize + length = src.RasterYSize + + driver = gdal.GetDriverByName("ENVI") + + low_ds = driver.Create( + str(low_path), + width, + length, + 1, + gdal.GDT_CFloat32, + ) + high_ds = driver.Create( + str(high_path), + width, + length, + 1, + gdal.GDT_CFloat32, + ) + + splitter = SplitSpectrum( + rg_sample_freq=self.meta.range_sample_frequency, + rg_bandwidth=self.meta.range_bandwidth, + center_frequency=self.meta.center_frequency, + slant_range=self.meta.slant_range, + freq=self.frequency, + ) + + fft_size = next_fast_len(width) + + src_band = src.GetRasterBand(1) + low_band = low_ds.GetRasterBand(1) + high_band = high_ds.GetRasterBand(1) + + for y0 in range(0, length, cfg.split_lines_per_block): + nlines = min( + cfg.split_lines_per_block, + length - y0, + ) + + block = src_band.ReadAsArray( + 0, + y0, + width, + nlines, + ).astype( + np.complex64, + copy=False, + ) + + low, _ = splitter.bandpass_shift_spectrum( + slc_raster=block, + low_frequency=self.low_limits[0], + high_frequency=self.low_limits[1], + new_center_frequency=self.low_center_frequency, + window_function=cfg.window_function, + window_shape=cfg.window_shape, + fft_size=fft_size, + resampling=False, + ) + + high, _ = splitter.bandpass_shift_spectrum( + slc_raster=block, + low_frequency=self.high_limits[0], + high_frequency=self.high_limits[1], + new_center_frequency=self.high_center_frequency, + window_function=cfg.window_function, + window_shape=cfg.window_shape, + fft_size=fft_size, + resampling=False, + ) + + low_band.WriteArray( + np.asarray(low, dtype=np.complex64), + 0, + y0, + ) + high_band.WriteArray( + np.asarray(high, dtype=np.complex64), + 0, + y0, + ) + + low_ds.FlushCache() + high_ds.FlushCache() + + low_ds = None + high_ds = None + src = None + + def _unwrap_pair( + self, + ifg_path: Path, + coh_path: Path, + outdir: Path, + ) -> tuple[Path, Path]: + algo = self.config.unwrap_algorithm.lower() + + if algo == "snaphu": + return self._unwrap_snaphu(ifg_path, coh_path, outdir) + + if algo == "icu": + return self._unwrap_icu(ifg_path, coh_path, outdir) + + if algo == "whirlwind": + return self._unwrap_whirlwind(ifg_path, coh_path, outdir) + + raise ValueError( + f"Unsupported unwrap_algorithm={self.config.unwrap_algorithm!r}. " + "Choose 'snaphu', 'icu', or 'whirlwind'." + ) + + def _unwrap_snaphu( + self, + ifg_path: Path, + coh_path: Path, + outdir: Path, + ) -> tuple[Path, Path]: + import snaphu + + cfg = self.config + outdir.mkdir(parents=True, exist_ok=True) + + unw_path = outdir / "unwrappedPhase.unw" + cc_path = outdir / "connectedComponents.conncomp" + + self._prepare_output(unw_path) + self._prepare_output(cc_path) + + igram = self._read_array( + ifg_path, + np.complex64, + ) + coh = self._read_array( + coh_path, + np.float32, + ) + + nlooks = ( + cfg.snaphu_nlooks + if cfg.snaphu_nlooks is not None + else float(cfg.range_looks * cfg.azimuth_looks) + ) + + scratch = outdir / "snaphu_scratch" + + if scratch.exists() and cfg.overwrite: + shutil.rmtree(scratch) + + scratch.mkdir(parents=True, exist_ok=True) + + unw, conncomp = snaphu.unwrap( + igram, + coh, + nlooks=nlooks, + cost=cfg.snaphu_cost, + init=cfg.snaphu_init, + min_conncomp_frac=cfg.snaphu_min_conncomp_frac, + phase_grad_window=cfg.snaphu_phase_grad_window, + ntiles=cfg.snaphu_ntiles, + tile_overlap=cfg.snaphu_tile_overlap, + nproc=cfg.snaphu_nproc, + tile_cost_thresh=cfg.snaphu_tile_cost_thresh, + min_region_size=cfg.snaphu_min_region_size, + single_tile_reoptimize=cfg.snaphu_single_tile_reoptimize, + scratchdir=scratch, + ) + + self._write_array( + unw_path, + np.asarray(unw, dtype=np.float32), + gdal.GDT_Float32, + ) + self._write_array( + cc_path, + np.asarray(conncomp, dtype=np.uint32), + gdal.GDT_UInt32, + ) + + return unw_path, cc_path + + def _unwrap_whirlwind( + self, + ifg_path: Path, + coh_path: Path, + outdir: Path, + ) -> tuple[Path, Path]: + """Unwrap with the ``whirlwind-insar`` Python package.""" + try: + import whirlwind as ww + except ImportError as exc: + raise ImportError( + "Whirlwind unwrapping requires the 'whirlwind-insar' package. " + "Install it with `pip install whirlwind-insar` or " + "`conda install -c conda-forge whirlwind-insar`." + ) from exc + + cfg = self.config + outdir.mkdir(parents=True, exist_ok=True) + + unw_path = outdir / "unwrappedPhase.unw" + cc_path = outdir / "connectedComponents.conncomp" + + self._prepare_output(unw_path) + self._prepare_output(cc_path) + + igram = self._read_array(ifg_path, np.complex64) + coh = self._read_array(coh_path, np.float32) + + nlooks = ( + cfg.whirlwind_nlooks + if cfg.whirlwind_nlooks is not None + else float(cfg.range_looks * cfg.azimuth_looks) + ) + + kwargs = {"nlooks": float(nlooks)} + if cfg.whirlwind_downsample is not None: + if cfg.whirlwind_downsample < 1: + raise ValueError("whirlwind_downsample must be >= 1") + kwargs["downsample"] = int(cfg.whirlwind_downsample) + + unw, conncomp = ww.unwrap(igram, coh, **kwargs) + + self._write_array( + unw_path, + np.asarray(unw, dtype=np.float32), + gdal.GDT_Float32, + ) + self._write_array( + cc_path, + np.asarray(conncomp, dtype=np.uint32), + gdal.GDT_UInt32, + ) + + return unw_path, cc_path + + def _unwrap_icu( + self, + ifg_path: Path, + coh_path: Path, + outdir: Path, + ) -> tuple[Path, Path]: + unw_path = outdir / "unwrappedPhase.unw" + cc_path = outdir / "connectedComponents.conncomp" + + self._prepare_output(unw_path) + self._prepare_output(cc_path) + + igram = isce3.io.Raster(str(ifg_path)) + coh = isce3.io.Raster(str(coh_path)) + + unw = isce3.io.Raster( + str(unw_path), + igram.width, + igram.length, + 1, + gdal.GDT_Float32, + "ENVI", + ) + cc = isce3.io.Raster( + str(cc_path), + igram.width, + igram.length, + 1, + gdal.GDT_Byte, + "ENVI", + ) + + icu = isce3.unwrap.ICU() + icu.unwrap( + unw, + cc, + igram, + coh, + ) + + del unw + del cc + del igram + del coh + + return unw_path, cc_path + + def _bridge_phase(self, phase): + cfg = self.config + + return bridge_unwrapped_phase( + phase, + radius=cfg.bridge_radius, + min_num_pixel=cfg.bridge_minimum_samples, + erosion_size=cfg.bridge_erosion_size, + ramp_type=cfg.bridge_ramp_type, + deramp_max_num_sample=cfg.bridge_ramp_maximum_pixel, + ) + + def _make_iono_mask( + self, + estimator, + phi_low, + phi_high, + low_coh, + high_coh, + low_cc, + high_cc, + ): + cfg = self.config + mask = np.ones(phi_low.shape, dtype=bool) + + if "coherence" in cfg.iono_mask_types: + mask &= estimator.get_coherence_mask_array( + low_band_array=low_coh, + high_band_array=high_coh, + threshold=cfg.iono_coherence_threshold, + ) + + if "connected_components" in cfg.iono_mask_types: + mask &= estimator.get_conn_component_mask_array( + low_band_array=low_cc, + high_band_array=high_cc, + ) + + mask &= estimator.get_valid_area( + low_band_array=phi_low, + high_band_array=phi_high, + invalid_value=0, + ) + + mask &= estimator.get_valid_area( + low_band_array=low_coh, + high_band_array=high_coh, + invalid_value=0, + ) + + mask &= np.isfinite(phi_low) + mask &= np.isfinite(phi_high) + mask &= np.isfinite(low_coh) + mask &= np.isfinite(high_coh) + + return mask + + def _run_iono_filter( + self, + data_path: Path, + sigma_path: Path, + mask_path: Path, + out_path: Path, + out_sigma_path: Path, + filter_dir: Path, + ): + cfg = self.config + filter_dir.mkdir(parents=True, exist_ok=True) + + filt = IonosphereFilter( + x_kernel=cfg.iono_kernel_range, + y_kernel=cfg.iono_kernel_azimuth, + sig_x=cfg.iono_sigma_range, + sig_y=cfg.iono_sigma_azimuth, + iteration=cfg.iono_filter_iterations, + filling_method=cfg.iono_filling_method, + guide_filter_method="median_gaussian", + guide_median_size=3, + outlier_threshold=3.5, + outlier_min_scale=0.0, + mad_scale_factor=1.4826, + outputdir=str(filter_dir), + ) + + self._remove_envi(out_path) + self._remove_envi(out_sigma_path) + + filt.low_pass_filter( + input_data=str(data_path), + input_std_dev=str(sigma_path), + mask_path=str(mask_path), + filtered_output=str(out_path), + filtered_std_dev=str(out_sigma_path), + lines_per_block=1000, + min_cluster_pixels=cfg.iono_min_cluster_pixels, + ) + + @staticmethod + def _open_gdal(path: str | Path): + ds = gdal.Open(str(path), gdal.GA_ReadOnly) + + if ds is None: + raise RuntimeError( + f"GDAL cannot open {path}" + ) + + return ds + + @classmethod + def _read_array(cls, path: str | Path, dtype=None): + ds = cls._open_gdal(path) + arr = ds.GetRasterBand(1).ReadAsArray() + ds = None + + if dtype is not None: + arr = arr.astype(dtype, copy=False) + + return arr + + def _write_array( + self, + path: str | Path, + array, + gdal_type, + ): + path = Path(path) + self._prepare_output(path) + + path.parent.mkdir(parents=True, exist_ok=True) + + arr = np.asarray(array) + length, width = arr.shape + + ds = gdal.GetDriverByName("ENVI").Create( + str(path), + width, + length, + 1, + gdal_type, + ) + + ds.GetRasterBand(1).WriteArray(arr) + ds.FlushCache() + ds = None + + def _prepare_output(self, path: Path): + path.parent.mkdir(parents=True, exist_ok=True) + + if path.exists(): + if not self.config.overwrite: + raise FileExistsError( + f"{path} exists. Set WorkflowConfig(overwrite=True) " + "to replace existing outputs." + ) + + self._remove_envi(path) + + @staticmethod + def _remove_envi(path: str | Path): + path = Path(path) + + for target in ( + path, + Path(str(path) + ".hdr"), + Path(str(path) + ".aux.xml"), + ): + if target.exists(): + target.unlink() + + +__all__ = [ + "AlignedPairIonoWorkflow", + "RadarMeta", + "WorkflowConfig", + "WorkflowProducts", +]