From eca7cac9ca3f9ff11c55037b8a39c937231c52a6 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 23 Jan 2026 00:24:10 +0000 Subject: [PATCH 01/17] Add Capella Space SAR sensor support for stripmapStack Add support for processing Capella Space SAR SLC data in the stripmapStack workflow. Capella provides stripmap mode X-band SAR data as GeoTIFF files with JSON metadata. New files: - Capella.py: Sensor class that parses Capella JSON metadata and extracts SLC images from GeoTIFF format - unpackFrame_Capella.py: Script to unpack Capella SLC data for stripmapStack - prepSlcCapella.py: Script to organize Capella data into date folders and generate run files Updates: - Register Capella sensor in Sensor/__init__.py - Add Capella pattern detection in prepSlcSensors.py --- components/isceobj/Sensor/Capella.py | 326 ++++++++++++++++++ components/isceobj/Sensor/__init__.py | 4 +- contrib/stack/stripmapStack/prepSlcCapella.py | 181 ++++++++++ contrib/stack/stripmapStack/prepSlcSensors.py | 8 +- .../stripmapStack/unpackFrame_Capella.py | 107 ++++++ 5 files changed, 622 insertions(+), 4 deletions(-) create mode 100644 components/isceobj/Sensor/Capella.py create mode 100644 contrib/stack/stripmapStack/prepSlcCapella.py create mode 100644 contrib/stack/stripmapStack/unpackFrame_Capella.py diff --git a/components/isceobj/Sensor/Capella.py b/components/isceobj/Sensor/Capella.py new file mode 100644 index 000000000..722c5f875 --- /dev/null +++ b/components/isceobj/Sensor/Capella.py @@ -0,0 +1,326 @@ +#!/usr/bin/env python3 + +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Copyright 2025 California Institute of Technology. ALL RIGHTS RESERVED. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# United States Government Sponsorship acknowledged. This software is subject to +# U.S. export control laws and regulations and has been classified as 'EAR99 NLR' +# (No [Export] License Required except when exporting to an embargoed country, +# end user, or in support of a prohibited end use). By downloading this software, +# the user agrees to comply with all applicable U.S. export laws and regulations. +# The user has the responsibility to obtain export licenses, or other export +# authority as may be required before exporting this software to any 'EAR99' +# embargoed foreign country or citizen of those countries. +# +# Author: Scott Staniewicz +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +import datetime +import json +import os +import numpy as np + +import isceobj +from isceobj.Scene.Frame import Frame +from isceobj.Planet.Planet import Planet +from isceobj.Orbit.Orbit import StateVector, Orbit +from isceobj.Planet.AstronomicalHandbook import Const +from iscesys.Component.Component import Component +from iscesys.DateTimeUtil.DateTimeUtil import DateTimeUtil as DTUtil + +lookMap = {'RIGHT': -1, 'LEFT': 1} + +TIFF = Component.Parameter( + 'tiff', + public_name='TIFF', + default='', + type=str, + mandatory=True, + doc='Capella SLC GeoTIFF imagery file' +) + +METADATA = Component.Parameter( + 'metadataFile', + public_name='METADATA', + default='', + type=str, + mandatory=True, + doc='Capella extended JSON metadata file' +) + +from .Sensor import Sensor + + +class Capella(Sensor): + """ + A class representing Capella Space SAR SLC data. + """ + + family = 'capella' + logging_name = 'isce.sensor.Capella' + parameter_list = (TIFF, METADATA) + Sensor.parameter_list + + def __init__(self, family='', name=''): + super().__init__(family if family else self.__class__.family, name=name) + self.frame = Frame() + self.frame.configure() + self._metadata = None + self.doppler_coeff = None + self.azfmrate_coeff = None + + def getFrame(self): + return self.frame + + def parse(self): + """Parse the Capella metadata JSON file.""" + try: + with open(self.metadataFile, 'r') as fp: + self._metadata = json.load(fp) + except IOError as strerr: + print("IOError: %s" % strerr) + return + self.populateMetadata() + + def populateMetadata(self): + """Create metadata objects from the JSON metadata.""" + + collect = self._metadata.get('collect', {}) + radar = self._metadata.get('radar', {}) + state_vectors = self._metadata.get('state_vectors', []) + image = collect.get('image', {}) + + # Platform info + mission = collect.get('platform', 'Capella') + pointing = radar.get('pointing', 'right') + lookSide = lookMap.get(pointing.upper(), -1) + + # Radar parameters + frequency = radar.get('center_frequency', 9.65e9) + txPol = radar.get('transmit_polarization', 'H') + rxPol = radar.get('receive_polarization', 'H') + polarization = txPol + rxPol + + # Get PRF from time_varying_parameters or prf array + prf_list = radar.get('prf', []) + if prf_list: + prf = prf_list[0].get('prf', 3000.0) + else: + tvp = radar.get('time_varying_parameters', []) + if tvp: + prf = tvp[0].get('prf', 3000.0) + else: + prf = 3000.0 + + # Pulse parameters from time_varying_parameters + tvp = radar.get('time_varying_parameters', []) + if tvp: + pulseBandwidth = tvp[0].get('pulse_bandwidth', 500e6) + pulseLength = tvp[0].get('pulse_duration', 10e-6) + else: + pulseBandwidth = 500e6 + pulseLength = 10e-6 + + # Sampling rate + samplingFrequency = radar.get('sampling_frequency', 600e6) + + # Image dimensions + lines = image.get('rows', image.get('length', 0)) + samples = image.get('columns', image.get('width', 0)) + + # Pixel spacing + rangePixelSize = image.get('pixel_spacing_column', 0.5) + azimuthPixelSize = image.get('pixel_spacing_row', 0.5) + + # Image geometry for starting range + image_geometry = image.get('image_geometry', {}) + range_time_origin = image_geometry.get('range_time_origin', 0.0) + startingRange = range_time_origin * Const.c / 2.0 + + # Incidence angle + center_pixel = image.get('center_pixel', {}) + incidenceAngle = center_pixel.get('incidence_angle', 35.0) + + # Timing information + startTimeStr = collect.get('start_timestamp', '') + stopTimeStr = collect.get('stop_timestamp', '') + dataStartTime = self._parseDateTime(startTimeStr) + dataStopTime = self._parseDateTime(stopTimeStr) + + # Pass direction from state vectors + if len(state_vectors) >= 2: + z0 = state_vectors[0].get('position', [0, 0, 0])[2] + z1 = state_vectors[-1].get('position', [0, 0, 0])[2] + passDirection = 'Ascending' if z1 > z0 else 'Descending' + else: + passDirection = 'Descending' + + # Populate platform + platform = self.frame.getInstrument().getPlatform() + platform.setPlanet(Planet(pname="Earth")) + platform.setMission(mission) + platform.setPointingDirection(lookSide) + platform.setAntennaLength(3.0) # Capella uses ~3m antenna + + # Populate instrument + instrument = self.frame.getInstrument() + instrument.setRadarFrequency(frequency) + instrument.setPulseRepetitionFrequency(prf) + instrument.setPulseLength(pulseLength) + instrument.setChirpSlope(pulseBandwidth / pulseLength) + instrument.setIncidenceAngle(incidenceAngle) + instrument.setRangePixelSize(rangePixelSize) + instrument.setRangeSamplingRate(samplingFrequency) + + # Populate Frame + self.frame.setSensingStart(dataStartTime) + self.frame.setSensingStop(dataStopTime) + diffTime = DTUtil.timeDeltaToSeconds(dataStopTime - dataStartTime) / 2.0 + sensingMid = dataStartTime + datetime.timedelta(microseconds=int(diffTime * 1e6)) + self.frame.setSensingMid(sensingMid) + self.frame.setPassDirection(passDirection) + self.frame.setPolarization(polarization) + self.frame.setStartingRange(startingRange) + self.frame.setFarRange(startingRange + (samples - 1) * rangePixelSize) + self.frame.setNumberOfLines(lines) + self.frame.setNumberOfSamples(samples) + self.frame.setProcessingFacility('Capella Space') + self.frame.setProcessingSoftwareVersion(self._metadata.get('software_version', '')) + + # Extract orbit + self.extractOrbit(state_vectors) + + # Save Doppler centroid coefficients + # Capella provides doppler_centroid_polynomial in the image_geometry + dc_poly = image_geometry.get('doppler_centroid_polynomial', [0.0]) + if isinstance(dc_poly, list) and len(dc_poly) > 0: + # Convert from Hz/m to Hz/pixel if needed + self.doppler_coeff = dc_poly + else: + self.doppler_coeff = [0.0] + + # FM rate is often not directly provided; use zeros as placeholder + self.azfmrate_coeff = [0.0, 0.0, 0.0] + + def _parseDateTime(self, timeStr): + """Parse Capella timestamp string to datetime object.""" + if not timeStr: + return datetime.datetime.now() + + # Handle various timestamp formats + # e.g., "2025-10-31T19:11:04.123456Z" or "2025-10-31T19:11:04Z" + try: + # Try with microseconds + if '.' in timeStr: + # Remove 'Z' suffix if present + timeStr = timeStr.rstrip('Z') + # Truncate to 6 decimal places for microseconds + parts = timeStr.split('.') + if len(parts) == 2: + timeStr = parts[0] + '.' + parts[1][:6] + return datetime.datetime.strptime(timeStr, "%Y-%m-%dT%H:%M:%S.%f") + else: + timeStr = timeStr.rstrip('Z') + return datetime.datetime.strptime(timeStr, "%Y-%m-%dT%H:%M:%S") + except ValueError: + print(f"Warning: Could not parse timestamp {timeStr}") + return datetime.datetime.now() + + def extractOrbit(self, state_vectors): + """Extract orbit state vectors from the metadata.""" + orbit = self.frame.getOrbit() + orbit.setOrbitSource('Header') + orbit.setReferenceFrame('ECR') + + for sv in state_vectors: + vec = StateVector() + timeStr = sv.get('time', '') + vec.setTime(self._parseDateTime(timeStr)) + position = sv.get('position', [0, 0, 0]) + velocity = sv.get('velocity', [0, 0, 0]) + vec.setPosition(position) + vec.setVelocity(velocity) + orbit.addStateVector(vec) + + def extractImage(self): + """ + Use GDAL to extract the SLC from Capella GeoTIFF. + """ + try: + from osgeo import gdal + except ImportError: + raise Exception('GDAL python bindings not found. Need this for Capella data.') + + self.parse() + + width = self.frame.getNumberOfSamples() + lgth = self.frame.getNumberOfLines() + + src = gdal.Open(self.tiff.strip(), gdal.GA_ReadOnly) + + # Capella SLC GeoTIFFs are complex int16 (CInt16) stored as a single band + # or as two bands (real/imag) + nbands = src.RasterCount + + if nbands == 1: + # Single band complex + data = src.GetRasterBand(1).ReadAsArray(0, 0, width, lgth) + if data.dtype != np.complex64: + data = data.astype(np.complex64) + elif nbands == 2: + # Two bands: real and imaginary + real = src.GetRasterBand(1).ReadAsArray(0, 0, width, lgth) + imag = src.GetRasterBand(2).ReadAsArray(0, 0, width, lgth) + cJ = np.complex64(1.0j) + data = real.astype(np.float32) + cJ * imag.astype(np.float32) + else: + raise Exception(f'Unexpected number of bands in Capella SLC: {nbands}') + + src = None + + data.tofile(self.output) + + # Create SLC image object + slcImage = isceobj.createSlcImage() + slcImage.setByteOrder('l') + slcImage.setFilename(self.output) + slcImage.setAccessMode('read') + slcImage.setWidth(width) + slcImage.setLength(lgth) + slcImage.setXmin(0) + slcImage.setXmax(width) + self.frame.setImage(slcImage) + + def extractDoppler(self): + """ + Extract Doppler centroid information. + """ + quadratic = {} + + # For insarApp style + prf = self.frame.getInstrument().getPulseRepetitionFrequency() + if self.doppler_coeff and len(self.doppler_coeff) > 0: + fd_mid = self.doppler_coeff[0] + else: + fd_mid = 0.0 + + quadratic['a'] = fd_mid / prf + quadratic['b'] = 0. + quadratic['c'] = 0. + + # For roiApp - doppler vs pixel + self.frame._dopplerVsPixel = self.doppler_coeff if self.doppler_coeff else [0.0] + print('Doppler Fit: ', self.frame._dopplerVsPixel) + + return quadratic diff --git a/components/isceobj/Sensor/__init__.py b/components/isceobj/Sensor/__init__.py index 15927bb0b..3da7d35c2 100755 --- a/components/isceobj/Sensor/__init__.py +++ b/components/isceobj/Sensor/__init__.py @@ -100,6 +100,7 @@ def factory_template(sat,name=None): createUAVSAR_Hdf5_SLC = partial(factory_template, 'UAVSAR_HDF5_SLC') createSAOCOM_SLC = partial(factory_template, 'SAOCOM_SLC') createLUTAN1 = partial(factory_template, 'Lutan1') +createCapella = partial(factory_template, 'Capella') SENSORS = {'ALOS' : createALOS, 'ALOS_SLC' : createALOS_SLC, @@ -127,7 +128,8 @@ def factory_template(sat,name=None): 'ICEYE_SLC' : createICEYE_SLC, 'UAVSAR_HDF5_SLC' : createUAVSAR_Hdf5_SLC, 'SAOCOM_SLC': createSAOCOM_SLC, - 'LUTAN1': createLUTAN1} + 'LUTAN1': createLUTAN1, + 'CAPELLA': createCapella} #These are experimental and can be added in as they become ready # 'JERS': createJERS, diff --git a/contrib/stack/stripmapStack/prepSlcCapella.py b/contrib/stack/stripmapStack/prepSlcCapella.py new file mode 100644 index 000000000..0a9c7b9c5 --- /dev/null +++ b/contrib/stack/stripmapStack/prepSlcCapella.py @@ -0,0 +1,181 @@ +#!/usr/bin/env python3 + +import os +import glob +import argparse +import json +from uncompressFile import uncompressfile + + +def createParser(): + ''' + Create command line parser. + ''' + + parser = argparse.ArgumentParser(description='Prepare Capella SLC processing (unzip/untar files, organize in date folders, generate script to unpack into ISCE formats).') + parser.add_argument('-i', '--input', dest='input', type=str, required=True, + help='directory with the Capella SLC data') + parser.add_argument('-rmfile', '--rmfile', dest='rmfile', action='store_true', default=False, + help='Optional: remove zip/tar/compressed files after unpacking into date structure (default is to keep in archive folder)') + parser.add_argument('-o', '--output', dest='output', type=str, required=False, + help='output directory where data needs to be unpacked into ISCE format (for script generation).') + parser.add_argument('-t', '--text_cmd', dest='text_cmd', type=str, default='source ~/.bash_profile;', + help='text command to be added to the beginning of each line of the run files. Default: source ~/.bash_profile;') + + return parser + + +def cmdLineParse(iargs=None): + ''' + Command line parser. + ''' + + parser = createParser() + return parser.parse_args(args=iargs) + + +def get_Date(capellaFolder): + ''' + Extract acquisition date from Capella metadata. + ''' + + # Look for extended JSON metadata file + jsonfiles = glob.glob(os.path.join(capellaFolder, '*_extended.json')) + if not jsonfiles: + jsonfiles = glob.glob(os.path.join(capellaFolder, '*.json')) + + if len(jsonfiles) > 0: + jsonfile = jsonfiles[0] + try: + with open(jsonfile, 'r') as fp: + metadata = json.load(fp) + + # Get timestamp from collect.start_timestamp + collect = metadata.get('collect', {}) + timestamp = collect.get('start_timestamp', '') + + if timestamp: + # Parse timestamp like "2025-10-31T19:11:04.123456Z" + acquisitionDate = timestamp[0:4] + timestamp[5:7] + timestamp[8:10] + if len(acquisitionDate) == 8: + return True, acquisitionDate + except (json.JSONDecodeError, IOError) as e: + print(f'Error reading {jsonfile}: {e}') + + # Could not find acquisition date + return False, 'FAIL' + + +def main(iargs=None): + ''' + The main driver. + ''' + + inps = cmdLineParse(iargs) + # parsing required inputs + inputDir = os.path.abspath(inps.input) + # parsing optional inputs + if inps.output: + outputDir = os.path.abspath(inps.output) + else: + outputDir = None + rmfile = inps.rmfile + + # filename of the runfile + run_unPack = 'run_unPackCapella' + + # Capella file patterns - handles various archive formats + # File naming: CAPELLA___SLC___. + capella_extensions = ( + os.path.join(inputDir, 'CAPELLA*.zip'), + os.path.join(inputDir, 'CAPELLA*.tar'), + os.path.join(inputDir, 'CAPELLA*.tar.gz'), + os.path.join(inputDir, 'CAPELLA*.tgz'), + ) + + for capella_extension in capella_extensions: + capella_filesfolders = glob.glob(capella_extension) + for capella_infilefolder in capella_filesfolders: + # the path to the folder/zip + workdir = os.path.dirname(capella_infilefolder) + + # get the output name folder without any extensions + temp = os.path.basename(capella_infilefolder) + # trim extensions + parts = temp.split('.') + capella_outfolder = parts[0] + # add the path back in + capella_outfolder = os.path.join(workdir, capella_outfolder) + + # this is a file, try to unzip/untar it + if os.path.isfile(capella_infilefolder): + # unzip the file in the outfolder + successflag_unzip = uncompressfile(capella_infilefolder, capella_outfolder) + + # put failed files in a separate directory + if not successflag_unzip: + os.makedirs(os.path.join(workdir, 'FAILED_FILES'), exist_ok=True) + os.rename(capella_infilefolder, os.path.join(workdir, 'FAILED_FILES', '.')) + else: + # check if file needs to be removed or put in archive folder + if rmfile: + os.remove(capella_infilefolder) + print('Deleting: ' + capella_infilefolder) + else: + os.makedirs(os.path.join(workdir, 'ARCHIVED_FILES'), exist_ok=True) + cmd = 'mv ' + capella_infilefolder + ' ' + os.path.join(workdir, 'ARCHIVED_FILES', '.') + os.system(cmd) + + # loop over the different Capella folders and organize in date folders + # Look for folders containing Capella data (has TIFF and JSON files) + capella_folders = glob.glob(os.path.join(inputDir, 'CAPELLA*')) + for capella_folder in capella_folders: + if not os.path.isdir(capella_folder): + continue + + # Verify it has the expected Capella files + tiffiles = glob.glob(os.path.join(capella_folder, '*.tif')) + jsonfiles = glob.glob(os.path.join(capella_folder, '*.json')) + if not tiffiles or not jsonfiles: + continue + + # get the date + successflag, imgDate = get_Date(capella_folder) + + workdir = os.path.dirname(capella_folder) + if successflag: + # move the folder into the date folder + SLC_dir = os.path.join(workdir, imgDate, '') + if os.path.isdir(SLC_dir): + import shutil + shutil.rmtree(SLC_dir) + + cmd = 'mv ' + capella_folder + ' ' + SLC_dir + os.system(cmd) + + print('Success: ' + imgDate) + else: + print('Failed: ' + capella_folder) + + # now generate the unpacking script for all the date dirs + dateDirs = glob.glob(os.path.join(inputDir, '2*')) + if outputDir is not None: + f = open(run_unPack, 'w') + for dateDir in dateDirs: + # Check for Capella TIFF files + capellaFiles = glob.glob(os.path.join(dateDir, 'CAPELLA*.tif')) + if not capellaFiles: + capellaFiles = glob.glob(os.path.join(dateDir, '*.tif')) + if len(capellaFiles) > 0: + acquisitionDate = os.path.basename(dateDir) + slcDir = os.path.join(outputDir, acquisitionDate) + os.makedirs(slcDir, exist_ok=True) + cmd = 'unpackFrame_Capella.py -i ' + os.path.abspath(dateDir) + ' -o ' + slcDir + print(cmd) + f.write(inps.text_cmd + cmd + '\n') + f.close() + + +if __name__ == '__main__': + + main() diff --git a/contrib/stack/stripmapStack/prepSlcSensors.py b/contrib/stack/stripmapStack/prepSlcSensors.py index 7f81b46c0..c435ec27e 100755 --- a/contrib/stack/stripmapStack/prepSlcSensors.py +++ b/contrib/stack/stripmapStack/prepSlcSensors.py @@ -65,11 +65,13 @@ def main(iargs=None): RSAT2_str2 = 'imagery_HH.tif' # RSAT2 extracted files TSX_TDX_str = 'dims_op*' # TSX zip files TSX_TDX_str2 = 'T*X*.xml' # TSX extracted files + CAPELLA_str = 'CAPELLA*SLC*.tif' # Capella SLC tiff files + CAPELLA_str2 = 'CAPELLA*_extended.json' # Capella extended metadata files # combine together - sensor_str_list = (ENV_str,ALOS1_str,CSK_str,CSK_str2,RSAT2_str,RSAT2_str2,TSX_TDX_str,TSX_TDX_str2) - sensor_list = ('Envisat','ALOS1','CSK','CSK','RSAT2','RSAT2','TSX/TDX','TSX/TDX') - sensor_unpackcommand = ('TODO','TODO','TODO','TODO','prepSlcRSAT2.py','prepSlcRSAT2.py','TODO','TODO') + sensor_str_list = (ENV_str,ALOS1_str,CSK_str,CSK_str2,RSAT2_str,RSAT2_str2,TSX_TDX_str,TSX_TDX_str2,CAPELLA_str,CAPELLA_str2) + sensor_list = ('Envisat','ALOS1','CSK','CSK','RSAT2','RSAT2','TSX/TDX','TSX/TDX','Capella','Capella') + sensor_unpackcommand = ('TODO','TODO','TODO','TODO','prepSlcRSAT2.py','prepSlcRSAT2.py','TODO','TODO','prepSlcCapella.py','prepSlcCapella.py') Sensors = dict(zip(sensor_str_list,sensor_list)) Sensors_unpack = dict(zip(sensor_str_list,sensor_unpackcommand)) diff --git a/contrib/stack/stripmapStack/unpackFrame_Capella.py b/contrib/stack/stripmapStack/unpackFrame_Capella.py new file mode 100644 index 000000000..9ea646ffe --- /dev/null +++ b/contrib/stack/stripmapStack/unpackFrame_Capella.py @@ -0,0 +1,107 @@ +#!/usr/bin/env python3 + +import isce +from isceobj.Sensor import createSensor +import shelve +import argparse +import glob +from isceobj.Util import Poly1D +from isceobj.Planet.AstronomicalHandbook import Const +from isceobj.Util.decorators import use_api +import os + + +def cmdLineParse(): + ''' + Command line parser. + ''' + + parser = argparse.ArgumentParser(description='Unpack Capella SLC data and store metadata in pickle file.') + parser.add_argument('-i', '--input', dest='capellaDir', type=str, + required=True, help='Input Capella SLC directory') + parser.add_argument('-o', '--output', dest='slcdir', type=str, + required=True, help='Output unpacked SLC directory') + + return parser.parse_args() + + +@use_api +def unpack(capellaDir, slcname): + ''' + Unpack Capella data to binary SLC file. + ''' + + # Search for imagery (GeoTIFF) and metadata (JSON) files in input directory + # Capella file naming: CAPELLA___SLC___.tif + # and CAPELLA___SLC____extended.json + tiffiles = glob.glob(os.path.join(capellaDir, 'CAPELLA*.tif')) + if not tiffiles: + tiffiles = glob.glob(os.path.join(capellaDir, '*.tif')) + + if not tiffiles: + raise FileNotFoundError(f'No TIFF files found in {capellaDir}') + + imgname = tiffiles[0] + + # Look for the extended JSON metadata file + jsonfiles = glob.glob(os.path.join(capellaDir, '*_extended.json')) + if not jsonfiles: + # Try any JSON file + jsonfiles = glob.glob(os.path.join(capellaDir, '*.json')) + + if not jsonfiles: + raise FileNotFoundError(f'No JSON metadata files found in {capellaDir}') + + metaname = jsonfiles[0] + + # Create output SLC directory if needed + if not os.path.isdir(slcname): + os.mkdir(slcname) + + date = os.path.basename(slcname) + + # Create a Capella sensor object and configure it + obj = createSensor('Capella') + obj.configure() + obj.metadataFile = metaname + obj.tiff = imgname + obj.output = os.path.join(slcname, date + '.slc') + + # Extract the image and write the XML file for the SLC + obj.extractImage() + obj.frame.getImage().renderHdr() + + # Save the doppler polynomial + coeffs = obj.doppler_coeff + poly = Poly1D.Poly1D() + poly.initPoly(order=len(coeffs) - 1) + poly.setCoeffs(coeffs) + + # Save the FM rate polynomial + fcoeffs = obj.azfmrate_coeff + fpoly = Poly1D.Poly1D() + fpoly.initPoly(order=len(fcoeffs) - 1) + fpoly.setCoeffs(fcoeffs) + + # Save required metadata for further use + # All data is output to a shelve file + pickName = os.path.join(slcname, 'data') + with shelve.open(pickName) as db: + db['frame'] = obj.frame + db['doppler'] = poly + db['fmrate'] = fpoly + + +if __name__ == '__main__': + ''' + Main driver. + ''' + + inps = cmdLineParse() + if inps.slcdir.endswith('/'): + inps.slcdir = inps.slcdir[:-1] + + if inps.capellaDir.endswith('/'): + inps.capellaDir = inps.capellaDir[:-1] + + unpack(inps.capellaDir, inps.slcdir) From 14d6f587865c8c6d50051fa14e326d11443a9549 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 23 Jan 2026 01:46:49 +0000 Subject: [PATCH 02/17] Remove fake FM rate polynomial from Capella sensor Capella does not provide FM rate polynomial in their metadata, so remove the placeholder zeros. The FM rate is only used by deskewALOS2.py for ALOS2-specific processing and is not needed for general stripmapStack. --- components/isceobj/Sensor/Capella.py | 5 ----- contrib/stack/stripmapStack/unpackFrame_Capella.py | 11 ++--------- 2 files changed, 2 insertions(+), 14 deletions(-) diff --git a/components/isceobj/Sensor/Capella.py b/components/isceobj/Sensor/Capella.py index 722c5f875..05f7194bf 100644 --- a/components/isceobj/Sensor/Capella.py +++ b/components/isceobj/Sensor/Capella.py @@ -78,7 +78,6 @@ def __init__(self, family='', name=''): self.frame.configure() self._metadata = None self.doppler_coeff = None - self.azfmrate_coeff = None def getFrame(self): return self.frame @@ -205,14 +204,10 @@ def populateMetadata(self): # Capella provides doppler_centroid_polynomial in the image_geometry dc_poly = image_geometry.get('doppler_centroid_polynomial', [0.0]) if isinstance(dc_poly, list) and len(dc_poly) > 0: - # Convert from Hz/m to Hz/pixel if needed self.doppler_coeff = dc_poly else: self.doppler_coeff = [0.0] - # FM rate is often not directly provided; use zeros as placeholder - self.azfmrate_coeff = [0.0, 0.0, 0.0] - def _parseDateTime(self, timeStr): """Parse Capella timestamp string to datetime object.""" if not timeStr: diff --git a/contrib/stack/stripmapStack/unpackFrame_Capella.py b/contrib/stack/stripmapStack/unpackFrame_Capella.py index 9ea646ffe..850a64dcf 100644 --- a/contrib/stack/stripmapStack/unpackFrame_Capella.py +++ b/contrib/stack/stripmapStack/unpackFrame_Capella.py @@ -71,25 +71,18 @@ def unpack(capellaDir, slcname): obj.extractImage() obj.frame.getImage().renderHdr() - # Save the doppler polynomial + # Save the doppler polynomial (Capella provides doppler_centroid_polynomial) coeffs = obj.doppler_coeff poly = Poly1D.Poly1D() poly.initPoly(order=len(coeffs) - 1) poly.setCoeffs(coeffs) - # Save the FM rate polynomial - fcoeffs = obj.azfmrate_coeff - fpoly = Poly1D.Poly1D() - fpoly.initPoly(order=len(fcoeffs) - 1) - fpoly.setCoeffs(fcoeffs) - # Save required metadata for further use - # All data is output to a shelve file + # Note: Capella does not provide FM rate polynomial, so we don't save it pickName = os.path.join(slcname, 'data') with shelve.open(pickName) as db: db['frame'] = obj.frame db['doppler'] = poly - db['fmrate'] = fpoly if __name__ == '__main__': From e1a965b53b4aadb1d5b3d74efe55b623259a6303 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 23 Jan 2026 10:40:21 +0000 Subject: [PATCH 03/17] Simplify Capella sensor: read metadata from TIFF, add mode validation - Read metadata from TIFF ImageDescription tag instead of requiring separate JSON files (metadata is embedded in the GeoTIFF) - Add mode validation to reject Spotlight (SP) and Sliding Spotlight (SL) modes which require additional post-processing to reramp phase - Simplify prepSlcCapella.py to work directly with .tif files - Update unpackFrame_Capella.py to only need TIFF file path - Remove JSON file pattern from prepSlcSensors.py auto-detection - Make scripts executable --- components/isceobj/Sensor/Capella.py | 75 ++++--- contrib/stack/stripmapStack/prepSlcCapella.py | 204 +++++++----------- contrib/stack/stripmapStack/prepSlcSensors.py | 9 +- .../stripmapStack/unpackFrame_Capella.py | 39 +--- 4 files changed, 144 insertions(+), 183 deletions(-) mode change 100644 => 100755 contrib/stack/stripmapStack/prepSlcCapella.py mode change 100644 => 100755 contrib/stack/stripmapStack/unpackFrame_Capella.py diff --git a/components/isceobj/Sensor/Capella.py b/components/isceobj/Sensor/Capella.py index 05f7194bf..23f25f992 100644 --- a/components/isceobj/Sensor/Capella.py +++ b/components/isceobj/Sensor/Capella.py @@ -29,19 +29,23 @@ import datetime import json -import os import numpy as np import isceobj from isceobj.Scene.Frame import Frame from isceobj.Planet.Planet import Planet -from isceobj.Orbit.Orbit import StateVector, Orbit +from isceobj.Orbit.Orbit import StateVector from isceobj.Planet.AstronomicalHandbook import Const from iscesys.Component.Component import Component from iscesys.DateTimeUtil.DateTimeUtil import DateTimeUtil as DTUtil +from .Sensor import Sensor + lookMap = {'RIGHT': -1, 'LEFT': 1} +# Supported modes for stripmapStack processing +SUPPORTED_MODES = {'stripmap', 'SM'} + TIFF = Component.Parameter( 'tiff', public_name='TIFF', @@ -51,26 +55,18 @@ doc='Capella SLC GeoTIFF imagery file' ) -METADATA = Component.Parameter( - 'metadataFile', - public_name='METADATA', - default='', - type=str, - mandatory=True, - doc='Capella extended JSON metadata file' -) - -from .Sensor import Sensor - class Capella(Sensor): """ A class representing Capella Space SAR SLC data. + + Metadata is read from the TIFF ImageDescription tag (embedded JSON). + Only stripmap (SM) mode is currently supported for stripmapStack processing. """ family = 'capella' logging_name = 'isce.sensor.Capella' - parameter_list = (TIFF, METADATA) + Sensor.parameter_list + parameter_list = (TIFF,) + Sensor.parameter_list def __init__(self, family='', name=''): super().__init__(family if family else self.__class__.family, name=name) @@ -82,16 +78,48 @@ def __init__(self, family='', name=''): def getFrame(self): return self.frame - def parse(self): - """Parse the Capella metadata JSON file.""" + def _loadMetadataFromTiff(self): + """Load metadata from TIFF ImageDescription tag using GDAL.""" try: - with open(self.metadataFile, 'r') as fp: - self._metadata = json.load(fp) - except IOError as strerr: - print("IOError: %s" % strerr) - return + from osgeo import gdal + except ImportError: + raise Exception('GDAL python bindings not found. Need this for Capella data.') + + ds = gdal.Open(self.tiff.strip(), gdal.GA_ReadOnly) + if ds is None: + raise Exception(f'Could not open TIFF file: {self.tiff}') + + # Get ImageDescription from TIFF metadata + image_desc = ds.GetMetadataItem('TIFFTAG_IMAGEDESCRIPTION') + ds = None + + if not image_desc: + raise Exception(f'No ImageDescription tag found in TIFF: {self.tiff}') + + try: + self._metadata = json.loads(image_desc) + except json.JSONDecodeError as e: + raise Exception(f'Failed to parse JSON from TIFF ImageDescription: {e}') + + def parse(self): + """Parse the Capella metadata from TIFF ImageDescription tag.""" + self._loadMetadataFromTiff() + self._validateMode() self.populateMetadata() + def _validateMode(self): + """Validate that the acquisition mode is supported (stripmap only).""" + collect = self._metadata.get('collect', {}) + mode = collect.get('mode', '') + + if mode.lower() not in {m.lower() for m in SUPPORTED_MODES}: + raise Exception( + f"Capella mode '{mode}' is not supported for stripmapStack. " + f"Only stripmap (SM) mode is supported. " + f"Spotlight (SP) and Sliding Spotlight (SL) modes require " + f"additional post-processing to reramp phase." + ) + def populateMetadata(self): """Create metadata objects from the JSON metadata.""" @@ -107,8 +135,8 @@ def populateMetadata(self): # Radar parameters frequency = radar.get('center_frequency', 9.65e9) - txPol = radar.get('transmit_polarization', 'H') - rxPol = radar.get('receive_polarization', 'H') + txPol = radar.get('transmit_polarization', 'V') + rxPol = radar.get('receive_polarization', 'V') polarization = txPol + rxPol # Get PRF from time_varying_parameters or prf array @@ -265,7 +293,6 @@ def extractImage(self): src = gdal.Open(self.tiff.strip(), gdal.GA_ReadOnly) # Capella SLC GeoTIFFs are complex int16 (CInt16) stored as a single band - # or as two bands (real/imag) nbands = src.RasterCount if nbands == 1: diff --git a/contrib/stack/stripmapStack/prepSlcCapella.py b/contrib/stack/stripmapStack/prepSlcCapella.py old mode 100644 new mode 100755 index 0a9c7b9c5..5246f535f --- a/contrib/stack/stripmapStack/prepSlcCapella.py +++ b/contrib/stack/stripmapStack/prepSlcCapella.py @@ -4,7 +4,6 @@ import glob import argparse import json -from uncompressFile import uncompressfile def createParser(): @@ -12,15 +11,13 @@ def createParser(): Create command line parser. ''' - parser = argparse.ArgumentParser(description='Prepare Capella SLC processing (unzip/untar files, organize in date folders, generate script to unpack into ISCE formats).') + parser = argparse.ArgumentParser(description='Prepare Capella SLC processing: organize TIFF files into date folders and generate unpack script.') parser.add_argument('-i', '--input', dest='input', type=str, required=True, - help='directory with the Capella SLC data') - parser.add_argument('-rmfile', '--rmfile', dest='rmfile', action='store_true', default=False, - help='Optional: remove zip/tar/compressed files after unpacking into date structure (default is to keep in archive folder)') + help='directory with the Capella SLC TIFF files') parser.add_argument('-o', '--output', dest='output', type=str, required=False, - help='output directory where data needs to be unpacked into ISCE format (for script generation).') - parser.add_argument('-t', '--text_cmd', dest='text_cmd', type=str, default='source ~/.bash_profile;', - help='text command to be added to the beginning of each line of the run files. Default: source ~/.bash_profile;') + help='output directory where SLCs will be unpacked into ISCE format') + parser.add_argument('-t', '--text_cmd', dest='text_cmd', type=str, default='', + help='text command to be added to the beginning of each line of the run files') return parser @@ -34,35 +31,39 @@ def cmdLineParse(iargs=None): return parser.parse_args(args=iargs) -def get_Date(capellaFolder): +def get_date_from_tiff(tiffFile): ''' - Extract acquisition date from Capella metadata. + Extract acquisition date from Capella TIFF ImageDescription tag. ''' + try: + from osgeo import gdal + except ImportError: + raise Exception('GDAL python bindings not found. Need this for Capella data.') + + ds = gdal.Open(tiffFile, gdal.GA_ReadOnly) + if ds is None: + return False, 'FAIL' + + # Get ImageDescription from TIFF metadata + image_desc = ds.GetMetadataItem('TIFFTAG_IMAGEDESCRIPTION') + ds = None + + if not image_desc: + return False, 'FAIL' + + try: + metadata = json.loads(image_desc) + collect = metadata.get('collect', {}) + timestamp = collect.get('start_timestamp', '') + + if timestamp: + # Parse timestamp like "2025-10-31T19:11:04.123456Z" + acquisitionDate = timestamp[0:4] + timestamp[5:7] + timestamp[8:10] + if len(acquisitionDate) == 8: + return True, acquisitionDate + except (json.JSONDecodeError, KeyError) as e: + print(f'Error reading metadata from {tiffFile}: {e}') - # Look for extended JSON metadata file - jsonfiles = glob.glob(os.path.join(capellaFolder, '*_extended.json')) - if not jsonfiles: - jsonfiles = glob.glob(os.path.join(capellaFolder, '*.json')) - - if len(jsonfiles) > 0: - jsonfile = jsonfiles[0] - try: - with open(jsonfile, 'r') as fp: - metadata = json.load(fp) - - # Get timestamp from collect.start_timestamp - collect = metadata.get('collect', {}) - timestamp = collect.get('start_timestamp', '') - - if timestamp: - # Parse timestamp like "2025-10-31T19:11:04.123456Z" - acquisitionDate = timestamp[0:4] + timestamp[5:7] + timestamp[8:10] - if len(acquisitionDate) == 8: - return True, acquisitionDate - except (json.JSONDecodeError, IOError) as e: - print(f'Error reading {jsonfile}: {e}') - - # Could not find acquisition date return False, 'FAIL' @@ -72,108 +73,67 @@ def main(iargs=None): ''' inps = cmdLineParse(iargs) - # parsing required inputs inputDir = os.path.abspath(inps.input) - # parsing optional inputs + if inps.output: outputDir = os.path.abspath(inps.output) else: outputDir = None - rmfile = inps.rmfile # filename of the runfile run_unPack = 'run_unPackCapella' - # Capella file patterns - handles various archive formats - # File naming: CAPELLA___SLC___. - capella_extensions = ( - os.path.join(inputDir, 'CAPELLA*.zip'), - os.path.join(inputDir, 'CAPELLA*.tar'), - os.path.join(inputDir, 'CAPELLA*.tar.gz'), - os.path.join(inputDir, 'CAPELLA*.tgz'), - ) - - for capella_extension in capella_extensions: - capella_filesfolders = glob.glob(capella_extension) - for capella_infilefolder in capella_filesfolders: - # the path to the folder/zip - workdir = os.path.dirname(capella_infilefolder) - - # get the output name folder without any extensions - temp = os.path.basename(capella_infilefolder) - # trim extensions - parts = temp.split('.') - capella_outfolder = parts[0] - # add the path back in - capella_outfolder = os.path.join(workdir, capella_outfolder) - - # this is a file, try to unzip/untar it - if os.path.isfile(capella_infilefolder): - # unzip the file in the outfolder - successflag_unzip = uncompressfile(capella_infilefolder, capella_outfolder) - - # put failed files in a separate directory - if not successflag_unzip: - os.makedirs(os.path.join(workdir, 'FAILED_FILES'), exist_ok=True) - os.rename(capella_infilefolder, os.path.join(workdir, 'FAILED_FILES', '.')) - else: - # check if file needs to be removed or put in archive folder - if rmfile: - os.remove(capella_infilefolder) - print('Deleting: ' + capella_infilefolder) - else: - os.makedirs(os.path.join(workdir, 'ARCHIVED_FILES'), exist_ok=True) - cmd = 'mv ' + capella_infilefolder + ' ' + os.path.join(workdir, 'ARCHIVED_FILES', '.') - os.system(cmd) - - # loop over the different Capella folders and organize in date folders - # Look for folders containing Capella data (has TIFF and JSON files) - capella_folders = glob.glob(os.path.join(inputDir, 'CAPELLA*')) - for capella_folder in capella_folders: - if not os.path.isdir(capella_folder): - continue - - # Verify it has the expected Capella files - tiffiles = glob.glob(os.path.join(capella_folder, '*.tif')) - jsonfiles = glob.glob(os.path.join(capella_folder, '*.json')) - if not tiffiles or not jsonfiles: - continue - - # get the date - successflag, imgDate = get_Date(capella_folder) - - workdir = os.path.dirname(capella_folder) - if successflag: - # move the folder into the date folder - SLC_dir = os.path.join(workdir, imgDate, '') - if os.path.isdir(SLC_dir): - import shutil - shutil.rmtree(SLC_dir) + # Find all Capella TIFF files + # File naming: CAPELLA___SLC___.tif + tiffFiles = glob.glob(os.path.join(inputDir, 'CAPELLA*.tif')) + if not tiffFiles: + tiffFiles = glob.glob(os.path.join(inputDir, '*.tif')) - cmd = 'mv ' + capella_folder + ' ' + SLC_dir - os.system(cmd) + # Organize TIFF files into date folders + for tiffFile in tiffFiles: + successflag, imgDate = get_date_from_tiff(tiffFile) - print('Success: ' + imgDate) + if successflag: + # Create date folder + dateDir = os.path.join(inputDir, imgDate) + os.makedirs(dateDir, exist_ok=True) + + # Move TIFF file to date folder + destFile = os.path.join(dateDir, os.path.basename(tiffFile)) + if tiffFile != destFile: + os.rename(tiffFile, destFile) + print(f'Moved {os.path.basename(tiffFile)} -> {imgDate}/') else: - print('Failed: ' + capella_folder) + print(f'Failed to get date from: {tiffFile}') + + # Generate the unpacking script for all date dirs + dateDirs = sorted(glob.glob(os.path.join(inputDir, '2*'))) - # now generate the unpacking script for all the date dirs - dateDirs = glob.glob(os.path.join(inputDir, '2*')) if outputDir is not None: - f = open(run_unPack, 'w') - for dateDir in dateDirs: - # Check for Capella TIFF files - capellaFiles = glob.glob(os.path.join(dateDir, 'CAPELLA*.tif')) - if not capellaFiles: - capellaFiles = glob.glob(os.path.join(dateDir, '*.tif')) - if len(capellaFiles) > 0: - acquisitionDate = os.path.basename(dateDir) - slcDir = os.path.join(outputDir, acquisitionDate) - os.makedirs(slcDir, exist_ok=True) - cmd = 'unpackFrame_Capella.py -i ' + os.path.abspath(dateDir) + ' -o ' + slcDir - print(cmd) - f.write(inps.text_cmd + cmd + '\n') - f.close() + with open(run_unPack, 'w') as f: + for dateDir in dateDirs: + if not os.path.isdir(dateDir): + continue + + # Find Capella TIFF file in this date directory + tiffFiles = glob.glob(os.path.join(dateDir, 'CAPELLA*.tif')) + if not tiffFiles: + tiffFiles = glob.glob(os.path.join(dateDir, '*.tif')) + + if tiffFiles: + tiffFile = tiffFiles[0] + acquisitionDate = os.path.basename(dateDir) + slcDir = os.path.join(outputDir, acquisitionDate) + os.makedirs(slcDir, exist_ok=True) + + cmd = f'unpackFrame_Capella.py -i {tiffFile} -o {slcDir}' + print(cmd) + if inps.text_cmd: + f.write(f'{inps.text_cmd} {cmd}\n') + else: + f.write(f'{cmd}\n') + + print(f'\nGenerated run file: {run_unPack}') if __name__ == '__main__': diff --git a/contrib/stack/stripmapStack/prepSlcSensors.py b/contrib/stack/stripmapStack/prepSlcSensors.py index c435ec27e..45e456e0b 100755 --- a/contrib/stack/stripmapStack/prepSlcSensors.py +++ b/contrib/stack/stripmapStack/prepSlcSensors.py @@ -65,13 +65,12 @@ def main(iargs=None): RSAT2_str2 = 'imagery_HH.tif' # RSAT2 extracted files TSX_TDX_str = 'dims_op*' # TSX zip files TSX_TDX_str2 = 'T*X*.xml' # TSX extracted files - CAPELLA_str = 'CAPELLA*SLC*.tif' # Capella SLC tiff files - CAPELLA_str2 = 'CAPELLA*_extended.json' # Capella extended metadata files + CAPELLA_str = 'CAPELLA*.tif' # Capella SLC GeoTIFF files # combine together - sensor_str_list = (ENV_str,ALOS1_str,CSK_str,CSK_str2,RSAT2_str,RSAT2_str2,TSX_TDX_str,TSX_TDX_str2,CAPELLA_str,CAPELLA_str2) - sensor_list = ('Envisat','ALOS1','CSK','CSK','RSAT2','RSAT2','TSX/TDX','TSX/TDX','Capella','Capella') - sensor_unpackcommand = ('TODO','TODO','TODO','TODO','prepSlcRSAT2.py','prepSlcRSAT2.py','TODO','TODO','prepSlcCapella.py','prepSlcCapella.py') + sensor_str_list = (ENV_str,ALOS1_str,CSK_str,CSK_str2,RSAT2_str,RSAT2_str2,TSX_TDX_str,TSX_TDX_str2,CAPELLA_str) + sensor_list = ('Envisat','ALOS1','CSK','CSK','RSAT2','RSAT2','TSX/TDX','TSX/TDX','Capella') + sensor_unpackcommand = ('TODO','TODO','TODO','TODO','prepSlcRSAT2.py','prepSlcRSAT2.py','TODO','TODO','prepSlcCapella.py') Sensors = dict(zip(sensor_str_list,sensor_list)) Sensors_unpack = dict(zip(sensor_str_list,sensor_unpackcommand)) diff --git a/contrib/stack/stripmapStack/unpackFrame_Capella.py b/contrib/stack/stripmapStack/unpackFrame_Capella.py old mode 100644 new mode 100755 index 850a64dcf..e4ce2ea7b --- a/contrib/stack/stripmapStack/unpackFrame_Capella.py +++ b/contrib/stack/stripmapStack/unpackFrame_Capella.py @@ -6,7 +6,6 @@ import argparse import glob from isceobj.Util import Poly1D -from isceobj.Planet.AstronomicalHandbook import Const from isceobj.Util.decorators import use_api import os @@ -17,8 +16,8 @@ def cmdLineParse(): ''' parser = argparse.ArgumentParser(description='Unpack Capella SLC data and store metadata in pickle file.') - parser.add_argument('-i', '--input', dest='capellaDir', type=str, - required=True, help='Input Capella SLC directory') + parser.add_argument('-i', '--input', dest='tiffFile', type=str, + required=True, help='Input Capella SLC GeoTIFF file') parser.add_argument('-o', '--output', dest='slcdir', type=str, required=True, help='Output unpacked SLC directory') @@ -26,34 +25,11 @@ def cmdLineParse(): @use_api -def unpack(capellaDir, slcname): +def unpack(tiffFile, slcname): ''' Unpack Capella data to binary SLC file. ''' - # Search for imagery (GeoTIFF) and metadata (JSON) files in input directory - # Capella file naming: CAPELLA___SLC___.tif - # and CAPELLA___SLC____extended.json - tiffiles = glob.glob(os.path.join(capellaDir, 'CAPELLA*.tif')) - if not tiffiles: - tiffiles = glob.glob(os.path.join(capellaDir, '*.tif')) - - if not tiffiles: - raise FileNotFoundError(f'No TIFF files found in {capellaDir}') - - imgname = tiffiles[0] - - # Look for the extended JSON metadata file - jsonfiles = glob.glob(os.path.join(capellaDir, '*_extended.json')) - if not jsonfiles: - # Try any JSON file - jsonfiles = glob.glob(os.path.join(capellaDir, '*.json')) - - if not jsonfiles: - raise FileNotFoundError(f'No JSON metadata files found in {capellaDir}') - - metaname = jsonfiles[0] - # Create output SLC directory if needed if not os.path.isdir(slcname): os.mkdir(slcname) @@ -63,8 +39,7 @@ def unpack(capellaDir, slcname): # Create a Capella sensor object and configure it obj = createSensor('Capella') obj.configure() - obj.metadataFile = metaname - obj.tiff = imgname + obj.tiff = tiffFile obj.output = os.path.join(slcname, date + '.slc') # Extract the image and write the XML file for the SLC @@ -94,7 +69,7 @@ def unpack(capellaDir, slcname): if inps.slcdir.endswith('/'): inps.slcdir = inps.slcdir[:-1] - if inps.capellaDir.endswith('/'): - inps.capellaDir = inps.capellaDir[:-1] + if inps.tiffFile.endswith('/'): + inps.tiffFile = inps.tiffFile[:-1] - unpack(inps.capellaDir, inps.slcdir) + unpack(inps.tiffFile, inps.slcdir) From aca72766e32f48206731cff4f1fce5fe1e3c9aad Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 23 Jan 2026 11:59:47 +0000 Subject: [PATCH 04/17] Fix Capella Doppler: convert 2D polynomial to 1D for isce2 Capella provides a 2D Doppler centroid polynomial (function of azimuth and range), but isce2 expects a 1D polynomial (function of range only). Fix by evaluating the 2D polynomial at a specific azimuth time (~1 second from start, or mid-scene for short acquisitions) to extract 1D coefficients as a function of range pixel. --- components/isceobj/Sensor/Capella.py | 72 +++++++++++++++++++++++++--- 1 file changed, 65 insertions(+), 7 deletions(-) diff --git a/components/isceobj/Sensor/Capella.py b/components/isceobj/Sensor/Capella.py index 23f25f992..cc60bdcee 100644 --- a/components/isceobj/Sensor/Capella.py +++ b/components/isceobj/Sensor/Capella.py @@ -228,13 +228,71 @@ def populateMetadata(self): # Extract orbit self.extractOrbit(state_vectors) - # Save Doppler centroid coefficients - # Capella provides doppler_centroid_polynomial in the image_geometry - dc_poly = image_geometry.get('doppler_centroid_polynomial', [0.0]) - if isinstance(dc_poly, list) and len(dc_poly) > 0: - self.doppler_coeff = dc_poly - else: - self.doppler_coeff = [0.0] + # Extract Doppler centroid coefficients + # Capella provides a 2D polynomial (azimuth, range), but isce2 needs 1D (range only) + # We evaluate the 2D poly at mid-azimuth to get 1D coefficients vs range pixel + self.doppler_coeff = self._extractDopplerCoeffs(image_geometry, lines, samples) + + def _extractDopplerCoeffs(self, image_geometry, lines, samples): + """ + Extract 1D Doppler coefficients from Capella's 2D polynomial. + + Capella provides a 2D polynomial: doppler(az, rg) = sum_{i,j} c[i][j] * az^i * rg^j + where az and rg are in pixel coordinates. + + ISCE2 expects a 1D polynomial as a function of range pixel. + We evaluate at mid-azimuth (~1 second from start or mid-scene) to get 1D coeffs. + """ + dc_poly = image_geometry.get('doppler_centroid_polynomial', {}) + + # Handle case where it's not a dict (legacy or missing) + if not isinstance(dc_poly, dict): + return [0.0] + + coeffs_2d = dc_poly.get('coefficients', [[0.0]]) + dimension = dc_poly.get('dimension', 1) + + # If it's actually 1D, just return the coefficients + if dimension == 1 or not isinstance(coeffs_2d[0], list): + if isinstance(coeffs_2d, list): + return coeffs_2d if coeffs_2d else [0.0] + return [0.0] + + # 2D polynomial: evaluate at a specific azimuth to get 1D in range + # Use ~1 second from start, or mid-scene if scene is shorter + delta_line_time = image_geometry.get('delta_line_time', 1.0 / 3000.0) + az_time_1sec = 1.0 # 1 second from start + az_pixel_1sec = az_time_1sec / delta_line_time + + # Use mid-scene if 1 second is past the scene + mid_az_pixel = lines / 2.0 + az_eval = min(az_pixel_1sec, mid_az_pixel) + + # Evaluate 2D poly at az_eval to get 1D coefficients in range + # doppler(rg) = sum_j [sum_i c[i][j] * az^i] * rg^j + # new_coeff[j] = sum_i c[i][j] * az^i + try: + coeffs_2d = np.array(coeffs_2d, dtype=np.float64) + degree_az = coeffs_2d.shape[0] - 1 + degree_rg = coeffs_2d.shape[1] - 1 + + # Compute 1D coefficients by evaluating at az_eval + coeffs_1d = [] + for j in range(degree_rg + 1): + coeff_j = 0.0 + for i in range(degree_az + 1): + coeff_j += coeffs_2d[i, j] * (az_eval ** i) + coeffs_1d.append(coeff_j) + + # Check if all coefficients are essentially zero + if all(abs(c) < 1e-15 for c in coeffs_1d): + return [0.0] + + return coeffs_1d + + except (ValueError, IndexError) as e: + print(f"Warning: Could not parse 2D Doppler polynomial: {e}") + return [0.0] def _parseDateTime(self, timeStr): """Parse Capella timestamp string to datetime object.""" From d7ac878285dd05fb162bca78455a7e6f3a0c517f Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 23 Jan 2026 12:27:45 +0000 Subject: [PATCH 05/17] Simplify Doppler conversion to always use midpoint azimuth --- components/isceobj/Sensor/Capella.py | 17 +++++------------ 1 file changed, 5 insertions(+), 12 deletions(-) diff --git a/components/isceobj/Sensor/Capella.py b/components/isceobj/Sensor/Capella.py index cc60bdcee..34cf04b34 100644 --- a/components/isceobj/Sensor/Capella.py +++ b/components/isceobj/Sensor/Capella.py @@ -241,7 +241,7 @@ def _extractDopplerCoeffs(self, image_geometry, lines, samples): where az and rg are in pixel coordinates. ISCE2 expects a 1D polynomial as a function of range pixel. - We evaluate at mid-azimuth (~1 second from start or mid-scene) to get 1D coeffs. + We evaluate at mid-azimuth to get 1D coeffs. """ dc_poly = image_geometry.get('doppler_centroid_polynomial', {}) @@ -258,17 +258,10 @@ def _extractDopplerCoeffs(self, image_geometry, lines, samples): return coeffs_2d if coeffs_2d else [0.0] return [0.0] - # 2D polynomial: evaluate at a specific azimuth to get 1D in range - # Use ~1 second from start, or mid-scene if scene is shorter - delta_line_time = image_geometry.get('delta_line_time', 1.0 / 3000.0) - az_time_1sec = 1.0 # 1 second from start - az_pixel_1sec = az_time_1sec / delta_line_time - - # Use mid-scene if 1 second is past the scene + # 2D polynomial: evaluate at mid-azimuth to get 1D in range mid_az_pixel = lines / 2.0 - az_eval = min(az_pixel_1sec, mid_az_pixel) - # Evaluate 2D poly at az_eval to get 1D coefficients in range + # Evaluate 2D poly at mid_az_pixel to get 1D coefficients in range # doppler(rg) = sum_j [sum_i c[i][j] * az^i] * rg^j # new_coeff[j] = sum_i c[i][j] * az^i try: @@ -276,12 +269,12 @@ def _extractDopplerCoeffs(self, image_geometry, lines, samples): degree_az = coeffs_2d.shape[0] - 1 degree_rg = coeffs_2d.shape[1] - 1 - # Compute 1D coefficients by evaluating at az_eval + # Compute 1D coefficients by evaluating at mid_az_pixel coeffs_1d = [] for j in range(degree_rg + 1): coeff_j = 0.0 for i in range(degree_az + 1): - coeff_j += coeffs_2d[i, j] * (az_eval ** i) + coeff_j += coeffs_2d[i, j] * (mid_az_pixel ** i) coeffs_1d.append(coeff_j) # Check if all coefficients are essentially zero From 186ee10098d8a72a94e780bd5c0b4c1292ff74c2 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Fri, 23 Jan 2026 07:41:35 -0500 Subject: [PATCH 06/17] ignore installed dirs --- .gitignore | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.gitignore b/.gitignore index a9a21d172..f5b5e05c5 100644 --- a/.gitignore +++ b/.gitignore @@ -9,3 +9,7 @@ config.log insar.log isce.log .ipynb_checkpoints + +build/ +install/ +data/ From 3b863940bbbede1f1d0003cfbd98fe67fa8c1cd5 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Fri, 23 Jan 2026 07:44:39 -0500 Subject: [PATCH 07/17] Fix Doppler 2D detection: check coefficient structure, not dimension field --- components/isceobj/Sensor/Capella.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/components/isceobj/Sensor/Capella.py b/components/isceobj/Sensor/Capella.py index 34cf04b34..a2ae09fe9 100644 --- a/components/isceobj/Sensor/Capella.py +++ b/components/isceobj/Sensor/Capella.py @@ -249,13 +249,15 @@ def _extractDopplerCoeffs(self, image_geometry, lines, samples): if not isinstance(dc_poly, dict): return [0.0] - coeffs_2d = dc_poly.get('coefficients', [[0.0]]) - dimension = dc_poly.get('dimension', 1) + coeffs = dc_poly.get('coefficients', [[0.0]]) - # If it's actually 1D, just return the coefficients - if dimension == 1 or not isinstance(coeffs_2d[0], list): - if isinstance(coeffs_2d, list): - return coeffs_2d if coeffs_2d else [0.0] + # Check if coefficients are 2D (list of lists) or 1D (flat list) + is_2d = isinstance(coeffs, list) and len(coeffs) > 0 and isinstance(coeffs[0], list) + + if not is_2d: + # Already 1D coefficients + if isinstance(coeffs, list): + return coeffs if coeffs else [0.0] return [0.0] # 2D polynomial: evaluate at mid-azimuth to get 1D in range @@ -265,7 +267,7 @@ def _extractDopplerCoeffs(self, image_geometry, lines, samples): # doppler(rg) = sum_j [sum_i c[i][j] * az^i] * rg^j # new_coeff[j] = sum_i c[i][j] * az^i try: - coeffs_2d = np.array(coeffs_2d, dtype=np.float64) + coeffs_2d = np.array(coeffs, dtype=np.float64) degree_az = coeffs_2d.shape[0] - 1 degree_rg = coeffs_2d.shape[1] - 1 From c685a007339b5d072c3613b063bd398f5a58c6e3 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Fri, 23 Jan 2026 08:45:01 -0500 Subject: [PATCH 08/17] Fix Capella metadata parsing for embedded TIFF format - State vectors are in collect.state.state_vectors, not top-level - Radar params are in collect.radar, not top-level - Handle dict format for position {x,y,z} and velocity {vx,vy,vz} - Use range_to_first_sample (meters) instead of range_time_origin --- components/isceobj/Sensor/Capella.py | 35 +++++++++++++++++++++------- 1 file changed, 27 insertions(+), 8 deletions(-) diff --git a/components/isceobj/Sensor/Capella.py b/components/isceobj/Sensor/Capella.py index a2ae09fe9..7792df59c 100644 --- a/components/isceobj/Sensor/Capella.py +++ b/components/isceobj/Sensor/Capella.py @@ -124,8 +124,11 @@ def populateMetadata(self): """Create metadata objects from the JSON metadata.""" collect = self._metadata.get('collect', {}) - radar = self._metadata.get('radar', {}) - state_vectors = self._metadata.get('state_vectors', []) + # Radar and state vectors can be at top level (extended JSON) or inside collect (TIFF embedded) + radar = self._metadata.get('radar', collect.get('radar', {})) + state_obj = self._metadata.get('state_vectors', collect.get('state', {}).get('state_vectors', [])) + # Handle both list format and nested state object + state_vectors = state_obj if isinstance(state_obj, list) else [] image = collect.get('image', {}) # Platform info @@ -172,8 +175,8 @@ def populateMetadata(self): # Image geometry for starting range image_geometry = image.get('image_geometry', {}) - range_time_origin = image_geometry.get('range_time_origin', 0.0) - startingRange = range_time_origin * Const.c / 2.0 + # range_to_first_sample is in meters + startingRange = image_geometry.get('range_to_first_sample', 0.0) # Incidence angle center_pixel = image.get('center_pixel', {}) @@ -187,8 +190,10 @@ def populateMetadata(self): # Pass direction from state vectors if len(state_vectors) >= 2: - z0 = state_vectors[0].get('position', [0, 0, 0])[2] - z1 = state_vectors[-1].get('position', [0, 0, 0])[2] + pos0 = state_vectors[0].get('position', [0, 0, 0]) + pos1 = state_vectors[-1].get('position', [0, 0, 0]) + z0 = pos0['z'] if isinstance(pos0, dict) else pos0[2] + z1 = pos1['z'] if isinstance(pos1, dict) else pos1[2] passDirection = 'Ascending' if z1 > z0 else 'Descending' else: passDirection = 'Descending' @@ -323,8 +328,22 @@ def extractOrbit(self, state_vectors): vec = StateVector() timeStr = sv.get('time', '') vec.setTime(self._parseDateTime(timeStr)) - position = sv.get('position', [0, 0, 0]) - velocity = sv.get('velocity', [0, 0, 0]) + + # Handle both list format [x, y, z] and dict format {x:, y:, z:} + pos = sv.get('position', [0, 0, 0]) + if isinstance(pos, dict): + position = [pos.get('x', 0), pos.get('y', 0), pos.get('z', 0)] + else: + position = pos + + vel = sv.get('velocity', [0, 0, 0]) + if isinstance(vel, dict): + velocity = [vel.get('vx', vel.get('x', 0)), + vel.get('vy', vel.get('y', 0)), + vel.get('vz', vel.get('z', 0))] + else: + velocity = vel + vec.setPosition(position) vec.setVelocity(velocity) orbit.addStateVector(vec) From db3811c4c1f520f25d36dfbc0242d3d188b9f097 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 23 Jan 2026 16:32:51 +0000 Subject: [PATCH 09/17] Add Capella.py to SConscript install list --- components/isceobj/Sensor/SConscript | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/isceobj/Sensor/SConscript b/components/isceobj/Sensor/SConscript index d526e2217..c2d845091 100644 --- a/components/isceobj/Sensor/SConscript +++ b/components/isceobj/Sensor/SConscript @@ -28,7 +28,7 @@ listFiles = ['ALOS.py','CEOS.py','COSMO_SkyMed.py','COSMO_SkyMed_SLC.py', 'UAVSAR_Polsar.py', 'ERS_EnviSAT.py', 'ICEYE_SLC.py', 'ALOS2.py', 'ERS_SLC.py', 'ALOS_SLC.py', 'EnviSAT_SLC.py', 'ERS_EnviSAT_SLC.py', 'SICD_RGZERO.py','UAVSAR_HDF5_SLC.py', - 'SAOCOM_SLC.py','__init__.py','Lutan1.py'] + 'SAOCOM_SLC.py','__init__.py','Lutan1.py','Capella.py'] helpList,installHelp = envSensor['HELP_BUILDER'](envSensor,'__init__.py',install) envSensor.Install(installHelp,helpList) From 1164444786e15c791f4e0953c7d6eb2fd08ae6b5 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Fri, 20 Feb 2026 22:32:15 -0500 Subject: [PATCH 10/17] Fix macOS shared library linking with conda-forge gfortran Explicitly link against libSystem on Apple platforms to ensure conda-forge gfortran produces valid shared libraries. Co-Authored-By: Claude Opus 4.6 --- .cmake/isce2_helpers.cmake | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.cmake/isce2_helpers.cmake b/.cmake/isce2_helpers.cmake index 69c7ccafe..d41183a8f 100644 --- a/.cmake/isce2_helpers.cmake +++ b/.cmake/isce2_helpers.cmake @@ -32,6 +32,11 @@ macro(isce2_add_cdll target) PREFIX "" OUTPUT_NAME ${target} SUFFIX .so) + # On macOS, conda-forge gfortran may produce dylibs without libSystem linked. + # Explicitly link against System to ensure valid shared libraries. + if(APPLE) + target_link_libraries(${target} PRIVATE "-lSystem") + endif() # If we're the root cmake project (e.g. not add_subdirectory): if("${CMAKE_SOURCE_DIR}" STREQUAL "${PROJECT_SOURCE_DIR}") From 9d0d8a241d22803a92a25751b04042bebcce7e7f Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Wed, 4 Mar 2026 21:14:06 -0500 Subject: [PATCH 11/17] Fix Capella metadata and optimize ampcor for stripmapStack Capella.py: Use image_geometry fields instead of collect-level metadata. - PRF from 1/delta_line_time (processed SLC line rate ~6209 Hz), not raw radar PRF (~9900 Hz) - Range pixel size from delta_range_sample (slant range ~0.617m), not pixel_spacing_column (ground spacing ~0.947m) - Sensing times from first_line_time, not collect start/stop timestamps - Sampling frequency derived from slant range pixel spacing These three fixes improve coherence from ~0.20 to ~0.40 for 3-day pairs. refineSecondaryTiming.py: Reduce ampcor window sizes and scale grid to image size for faster processing. - Window 64x64 (was 128x128), search 20px (was 40) - Grid count adapts to image size, max 10x10 (was hardcoded 60x60) - Unconditionally set numberLocationAcross/Down since Ampcor.configure() sets defaults that prevent conditional override Co-Authored-By: Claude Opus 4.6 --- components/isceobj/Sensor/Capella.py | 53 +++++++++---------- .../stripmapStack/refineSecondaryTiming.py | 29 +++++----- 2 files changed, 39 insertions(+), 43 deletions(-) diff --git a/components/isceobj/Sensor/Capella.py b/components/isceobj/Sensor/Capella.py index 7792df59c..6709d217c 100644 --- a/components/isceobj/Sensor/Capella.py +++ b/components/isceobj/Sensor/Capella.py @@ -142,17 +142,6 @@ def populateMetadata(self): rxPol = radar.get('receive_polarization', 'V') polarization = txPol + rxPol - # Get PRF from time_varying_parameters or prf array - prf_list = radar.get('prf', []) - if prf_list: - prf = prf_list[0].get('prf', 3000.0) - else: - tvp = radar.get('time_varying_parameters', []) - if tvp: - prf = tvp[0].get('prf', 3000.0) - else: - prf = 3000.0 - # Pulse parameters from time_varying_parameters tvp = radar.get('time_varying_parameters', []) if tvp: @@ -162,31 +151,40 @@ def populateMetadata(self): pulseBandwidth = 500e6 pulseLength = 10e-6 - # Sampling rate - samplingFrequency = radar.get('sampling_frequency', 600e6) - # Image dimensions lines = image.get('rows', image.get('length', 0)) samples = image.get('columns', image.get('width', 0)) - # Pixel spacing - rangePixelSize = image.get('pixel_spacing_column', 0.5) - azimuthPixelSize = image.get('pixel_spacing_row', 0.5) - - # Image geometry for starting range + # Image geometry defines the actual SLC pixel grid + # These are the ground truth for timing and range spacing image_geometry = image.get('image_geometry', {}) - # range_to_first_sample is in meters startingRange = image_geometry.get('range_to_first_sample', 0.0) + delta_range_sample = image_geometry.get('delta_range_sample', 0.0) + delta_line_time = image_geometry.get('delta_line_time', 0.0) + assert delta_range_sample > 0, f"Missing delta_range_sample in image_geometry" + assert delta_line_time > 0, f"Missing delta_line_time in image_geometry" + + # PRF = 1/delta_line_time (the processed SLC line rate, not the raw radar PRF) + prf = 1.0 / delta_line_time + + # Slant range pixel size from the SLC grid (not ground pixel spacing) + rangePixelSize = delta_range_sample + + # Sampling rate consistent with SLC range pixel spacing + samplingFrequency = Const.c / (2.0 * delta_range_sample) # Incidence angle center_pixel = image.get('center_pixel', {}) incidenceAngle = center_pixel.get('incidence_angle', 35.0) - # Timing information - startTimeStr = collect.get('start_timestamp', '') - stopTimeStr = collect.get('stop_timestamp', '') - dataStartTime = self._parseDateTime(startTimeStr) - dataStopTime = self._parseDateTime(stopTimeStr) + # Timing: use first_line_time from image_geometry (actual SLC grid start), + # NOT collect start/stop timestamps (which span the full acquisition window) + firstLineTimeStr = image_geometry.get('first_line_time', '') + assert firstLineTimeStr, "Missing first_line_time in image_geometry" + dataStartTime = self._parseDateTime(firstLineTimeStr) + dataStopTime = dataStartTime + datetime.timedelta( + seconds=(lines - 1) * delta_line_time + ) # Pass direction from state vectors if len(state_vectors) >= 2: @@ -218,8 +216,9 @@ def populateMetadata(self): # Populate Frame self.frame.setSensingStart(dataStartTime) self.frame.setSensingStop(dataStopTime) - diffTime = DTUtil.timeDeltaToSeconds(dataStopTime - dataStartTime) / 2.0 - sensingMid = dataStartTime + datetime.timedelta(microseconds=int(diffTime * 1e6)) + sensingMid = dataStartTime + datetime.timedelta( + seconds=(lines - 1) * delta_line_time / 2.0 + ) self.frame.setSensingMid(sensingMid) self.frame.setPassDirection(passDirection) self.frame.setPolarization(polarization) diff --git a/contrib/stack/stripmapStack/refineSecondaryTiming.py b/contrib/stack/stripmapStack/refineSecondaryTiming.py index fda3b0c9f..fd7c2238e 100755 --- a/contrib/stack/stripmapStack/refineSecondaryTiming.py +++ b/contrib/stack/stripmapStack/refineSecondaryTiming.py @@ -76,23 +76,23 @@ def estimateOffsetField(reference, secondary, azoffset=0, rgoffset=0): objOffset.configure() objOffset.setAcrossGrossOffset(rgoffset) objOffset.setDownGrossOffset(azoffset) - objOffset.setWindowSizeWidth(128) - objOffset.setWindowSizeHeight(128) - objOffset.setSearchWindowSizeWidth(40) - objOffset.setSearchWindowSizeHeight(40) + objOffset.setWindowSizeWidth(64) + objOffset.setWindowSizeHeight(64) + objOffset.setSearchWindowSizeWidth(20) + objOffset.setSearchWindowSizeHeight(20) margin = 2*objOffset.searchWindowSizeWidth + objOffset.windowSizeWidth - nAcross = 60 - nDown = 60 + offAc = max(51,-rgoffset)+margin + offDn = max(51,-azoffset)+margin - - offAc = max(101,-rgoffset)+margin - offDn = max(101,-azoffset)+margin - - lastAc = int( min(width, sim.getWidth() - offAc) - margin) lastDn = int( min(length, sim.getLength() - offDn) - margin) + # Scale window count to image size: enough for robust constant/low-order + # offset estimation, but not so many that ampcor takes hours. + nAcross = min(10, max(3, (lastAc - offAc) // objOffset.windowSizeWidth)) + nDown = min(10, max(3, (lastDn - offDn) // objOffset.windowSizeHeight)) + # print('Across: ', offAc, lastAc, width, sim.getWidth(), margin) # print('Down: ', offDn, lastDn, length, sim.getLength(), margin) @@ -108,11 +108,8 @@ def estimateOffsetField(reference, secondary, azoffset=0, rgoffset=0): if not objOffset.lastSampleDown: objOffset.setLastSampleDown(lastDn) - if not objOffset.numberLocationAcross: - objOffset.setNumberLocationAcross(nAcross) - - if not objOffset.numberLocationDown: - objOffset.setNumberLocationDown(nDown) + objOffset.setNumberLocationAcross(nAcross) + objOffset.setNumberLocationDown(nDown) objOffset.setFirstPRF(1.0) objOffset.setSecondPRF(1.0) From 15d5f6e40039f2be9e72894cca75f68382f35a33 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Wed, 4 Mar 2026 21:53:57 -0500 Subject: [PATCH 12/17] Revert refineSecondaryTiming.py changes to avoid affecting other sensors The ampcor window/grid optimization improves Capella processing speed but could affect other widely-used sensor workflows. Reverting to upstream defaults. Users processing Capella data with small images can pass --ao/--ro flags or reduce ampcor settings manually. Co-Authored-By: Claude Opus 4.6 --- .../stripmapStack/refineSecondaryTiming.py | 29 ++++++++++--------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/contrib/stack/stripmapStack/refineSecondaryTiming.py b/contrib/stack/stripmapStack/refineSecondaryTiming.py index fd7c2238e..fda3b0c9f 100755 --- a/contrib/stack/stripmapStack/refineSecondaryTiming.py +++ b/contrib/stack/stripmapStack/refineSecondaryTiming.py @@ -76,23 +76,23 @@ def estimateOffsetField(reference, secondary, azoffset=0, rgoffset=0): objOffset.configure() objOffset.setAcrossGrossOffset(rgoffset) objOffset.setDownGrossOffset(azoffset) - objOffset.setWindowSizeWidth(64) - objOffset.setWindowSizeHeight(64) - objOffset.setSearchWindowSizeWidth(20) - objOffset.setSearchWindowSizeHeight(20) + objOffset.setWindowSizeWidth(128) + objOffset.setWindowSizeHeight(128) + objOffset.setSearchWindowSizeWidth(40) + objOffset.setSearchWindowSizeHeight(40) margin = 2*objOffset.searchWindowSizeWidth + objOffset.windowSizeWidth - offAc = max(51,-rgoffset)+margin - offDn = max(51,-azoffset)+margin + nAcross = 60 + nDown = 60 + + offAc = max(101,-rgoffset)+margin + offDn = max(101,-azoffset)+margin + + lastAc = int( min(width, sim.getWidth() - offAc) - margin) lastDn = int( min(length, sim.getLength() - offDn) - margin) - # Scale window count to image size: enough for robust constant/low-order - # offset estimation, but not so many that ampcor takes hours. - nAcross = min(10, max(3, (lastAc - offAc) // objOffset.windowSizeWidth)) - nDown = min(10, max(3, (lastDn - offDn) // objOffset.windowSizeHeight)) - # print('Across: ', offAc, lastAc, width, sim.getWidth(), margin) # print('Down: ', offDn, lastDn, length, sim.getLength(), margin) @@ -108,8 +108,11 @@ def estimateOffsetField(reference, secondary, azoffset=0, rgoffset=0): if not objOffset.lastSampleDown: objOffset.setLastSampleDown(lastDn) - objOffset.setNumberLocationAcross(nAcross) - objOffset.setNumberLocationDown(nDown) + if not objOffset.numberLocationAcross: + objOffset.setNumberLocationAcross(nAcross) + + if not objOffset.numberLocationDown: + objOffset.setNumberLocationDown(nDown) objOffset.setFirstPRF(1.0) objOffset.setSecondPRF(1.0) From 799704e26bd7a028236da224afd05f0aa22ccc50 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Thu, 5 Mar 2026 10:49:43 -0500 Subject: [PATCH 13/17] Add Capella spotlight support and example notebook - Add spotlight mode (SP) to Capella.py supported modes and extract reference_antenna_position/reference_target_position from metadata for post-coregistration phase correction - Add spotlightPhaseCorrection.py for re-ramping deramped spotlight SLCs using geometric phase computation (ECEF-based) - Add example_capella_stripmap.py jupytext notebook showing full Capella stripmapStack workflow from SLC prep through interferograms Co-Authored-By: Claude Opus 4.6 --- components/isceobj/Sensor/Capella.py | 25 +- .../stripmapStack/example_capella_stripmap.py | 351 ++++++++++++++++++ .../stripmapStack/spotlightPhaseCorrection.py | 212 +++++++++++ 3 files changed, 583 insertions(+), 5 deletions(-) create mode 100755 contrib/stack/stripmapStack/example_capella_stripmap.py create mode 100755 contrib/stack/stripmapStack/spotlightPhaseCorrection.py diff --git a/components/isceobj/Sensor/Capella.py b/components/isceobj/Sensor/Capella.py index 6709d217c..9ea0fe394 100644 --- a/components/isceobj/Sensor/Capella.py +++ b/components/isceobj/Sensor/Capella.py @@ -44,7 +44,7 @@ lookMap = {'RIGHT': -1, 'LEFT': 1} # Supported modes for stripmapStack processing -SUPPORTED_MODES = {'stripmap', 'SM'} +SUPPORTED_MODES = {'stripmap', 'SM', 'spotlight', 'SP'} TIFF = Component.Parameter( 'tiff', @@ -108,16 +108,15 @@ def parse(self): self.populateMetadata() def _validateMode(self): - """Validate that the acquisition mode is supported (stripmap only).""" + """Validate that the acquisition mode is supported.""" collect = self._metadata.get('collect', {}) mode = collect.get('mode', '') if mode.lower() not in {m.lower() for m in SUPPORTED_MODES}: raise Exception( f"Capella mode '{mode}' is not supported for stripmapStack. " - f"Only stripmap (SM) mode is supported. " - f"Spotlight (SP) and Sliding Spotlight (SL) modes require " - f"additional post-processing to reramp phase." + f"Supported modes: {SUPPORTED_MODES}. " + f"Sliding Spotlight (SL) mode is not yet supported." ) def populateMetadata(self): @@ -237,6 +236,22 @@ def populateMetadata(self): # We evaluate the 2D poly at mid-azimuth to get 1D coefficients vs range pixel self.doppler_coeff = self._extractDopplerCoeffs(image_geometry, lines, samples) + # Spotlight-specific metadata for phase correction + # reference_antenna_position and reference_target_position are ECEF coordinates + # used to compute the geometric phase that must be restored after coregistration + mode = collect.get('mode', '') + if mode.lower() in {'spotlight', 'sp'}: + ref_ant = image.get('reference_antenna_position', {}) + ref_tgt = image.get('reference_target_position', {}) + assert ref_ant, "Spotlight mode requires reference_antenna_position in metadata" + assert ref_tgt, "Spotlight mode requires reference_target_position in metadata" + self.frame.spotlightReferenceAntennaPosition = [ + ref_ant['x'], ref_ant['y'], ref_ant['z'] + ] + self.frame.spotlightReferenceTargetPosition = [ + ref_tgt['x'], ref_tgt['y'], ref_tgt['z'] + ] + def _extractDopplerCoeffs(self, image_geometry, lines, samples): """ Extract 1D Doppler coefficients from Capella's 2D polynomial. diff --git a/contrib/stack/stripmapStack/example_capella_stripmap.py b/contrib/stack/stripmapStack/example_capella_stripmap.py new file mode 100755 index 000000000..cf5329394 --- /dev/null +++ b/contrib/stack/stripmapStack/example_capella_stripmap.py @@ -0,0 +1,351 @@ +# --- +# jupyter: +# jupytext: +# formats: py:percent +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# kernelspec: +# display_name: Python 3 +# language: python +# name: python3 +# --- + +# %% [markdown] +# # Capella Stripmap InSAR with ISCE2's stripmapStack +# +# End-to-end workflow for processing Capella Space stripmap SLC GeoTIFFs +# through ISCE2's stripmapStack processor to generate interferograms. +# +# ## Prerequisites +# - ISCE2 installed with Capella sensor support +# - Capella SLC GeoTIFFs (stripmap mode) +# - A DEM covering the area of interest (e.g. Copernicus DEM or 3DEP) +# - stripmapStack scripts on PATH (or run from source tree) + +# %% [markdown] +# ## 1. Setup + +# %% +import os +import glob +import subprocess +import time +from pathlib import Path + +import numpy as np +import matplotlib.pyplot as plt + +# %% [markdown] +# ### Configure paths +# +# Update these paths to point to your data. + +# %% +# Input Capella SLC GeoTIFFs +slc_tiff_dir = "/Volumes/WD_BLACK_SN7100_4TB/Documents/Learning/capella/mexico_city/cropped_slcs" + +# DEM file (GeoTIFF - will be converted to ISCE format) +dem_tiff = "/Volumes/WD_BLACK_SN7100_4TB/Documents/Learning/capella/mexico_city/dem_3dep.tif" + +# Working directory for ISCE2 processing +work_dir = os.path.expanduser("~/repos/isce2/data_notebook") +os.makedirs(work_dir, exist_ok=True) + +# Unpacked SLC output directory +slc_dir = os.path.join(work_dir, "SLC") +os.makedirs(slc_dir, exist_ok=True) + +# Path to stripmapStack scripts (if not on PATH, set this to source tree) +script_dir = os.path.expanduser("~/repos/isce2/contrib/stack/stripmapStack") + +# %% +# List available SLC TIFFs +tiffs = sorted(glob.glob(os.path.join(slc_tiff_dir, "CAPELLA*.tif"))) +print(f"Found {len(tiffs)} Capella SLC TIFFs:") +for t in tiffs: + print(f" {os.path.basename(t)}") + +# %% [markdown] +# ## 2. Prepare DEM +# +# Convert the GeoTIFF DEM to ISCE2 format using `dem_to_isce.py` or GDAL. +# The DEM must have `.xml` and `.vrt` sidecar files. + +# %% +dem_isce = os.path.join(work_dir, "dem.wgs84") + +if not os.path.exists(dem_isce + ".xml"): + subprocess.run( + ["gdal_translate", "-of", "ENVI", dem_tiff, dem_isce], + check=True + ) + # Create ISCE XML/VRT metadata + subprocess.run( + ["fixImageXml.py", "-f", "-i", dem_isce], + check=True + ) + print("DEM converted to ISCE format") +else: + print("DEM already exists in ISCE format") + +# %% [markdown] +# ## 3. Organize SLCs by date and unpack +# +# `prepSlcCapella.py` reads each TIFF's metadata to extract the acquisition +# date, moves files into date folders, and generates an unpacking run file. +# +# `unpackFrame_Capella.py` converts each GeoTIFF to ISCE2's binary SLC format +# and extracts orbit/timing metadata into a shelve file. + +# %% +# Step 1: Organize TIFFs into date folders and generate unpack commands +os.chdir(work_dir) +subprocess.run( + [ + "python", os.path.join(script_dir, "prepSlcCapella.py"), + "-i", slc_tiff_dir, + "-o", slc_dir, + ], + check=True +) + +# %% +# Step 2: Run the unpack commands +run_file = os.path.join(work_dir, "run_unPackCapella") +if os.path.exists(run_file): + with open(run_file) as f: + cmds = f.read().strip().split("\n") + + print(f"Unpacking {len(cmds)} SLCs...") + t0 = time.time() + for cmd in cmds: + print(f" {cmd}") + subprocess.run(cmd.split(), check=True) + print(f"Unpacking done in {time.time() - t0:.1f}s") + +# %% +# Verify unpacked SLCs +dates = sorted([d for d in os.listdir(slc_dir) if os.path.isdir(os.path.join(slc_dir, d))]) +print(f"\nUnpacked {len(dates)} dates: {dates}") + +# %% [markdown] +# ## 4. Generate stripmapStack run files +# +# `stackStripMap.py` creates the coregistration processing steps (run files). +# We use `--nofocus` since Capella data is already focused to SLC. + +# %% +os.chdir(work_dir) +subprocess.run( + [ + "python", os.path.join(script_dir, "stackStripMap.py"), + "-s", slc_dir, + "-d", dem_isce, + "-w", work_dir, + "-a", "1", + "-r", "1", + "--nofocus", + "-W", "slc", + ], + check=True +) + +# %% +# List generated run files +run_files = sorted(glob.glob(os.path.join(work_dir, "run_files", "run_*"))) +print("Generated run files:") +for rf in run_files: + print(f" {os.path.basename(rf)}") + +# %% [markdown] +# ## 5. Execute coregistration steps +# +# Run each step sequentially. Steps: +# 1. Reference geometry (topo) +# 2. Secondary focus/split (no-op for pre-focused SLCs) +# 3. Geo2rdr + coarse resampling +# 4. Refine secondary timing (ampcor) - **slowest step** +# 5. Invert misregistration +# 6. Fine resampling +# 7. Grid baseline +# +# **Note on step 4 (refineSecondaryTiming):** The default 60x60 ampcor grid +# can be very slow. For small test scenes, this is fine. For large scenes, +# you can reduce the grid by editing the config files to add `--nwa 10 --nwd 10`. + +# %% +os.chdir(work_dir) +total_t0 = time.time() + +for rf in run_files: + step_name = os.path.basename(rf) + print(f"\n{'='*60}") + print(f"Running: {step_name}") + print(f"{'='*60}") + + with open(rf) as f: + cmds = f.read().strip().split("\n") + + t0 = time.time() + for cmd in cmds: + if cmd.strip(): + subprocess.run(cmd, shell=True, check=True) + elapsed = time.time() - t0 + print(f" {step_name} completed in {elapsed:.1f}s") + +total_elapsed = time.time() - total_t0 +print(f"\nAll steps completed in {total_elapsed:.1f}s") + +# %% [markdown] +# ## 6. Generate interferograms +# +# Create interferograms for consecutive date pairs using `crossmul.py`. +# +# **Multi-looking ratio**: Capella stripmap has range pixel size ~0.617m and +# azimuth spacing ~1.2m. For approximately square pixels, use a 1:2 ratio +# (az:range looks). Example: 5 azimuth x 9 range gives ~5.8m x 5.6m pixels. + +# %% +# Find all coregistered SLC dates +merged_dir = os.path.join(work_dir, "merged") +coreg_dir = os.path.join(merged_dir, "SLC") + +# Reference date +ref_date = dates[0] +ref_slc = os.path.join(coreg_dir, ref_date, ref_date + ".slc") + +all_dates = sorted(os.listdir(coreg_dir)) +print(f"Coregistered dates: {all_dates}") + +# %% +# Generate interferograms for consecutive pairs +igram_dir = os.path.join(work_dir, "interferograms") +os.makedirs(igram_dir, exist_ok=True) + +pairs = list(zip(all_dates[:-1], all_dates[1:])) +print(f"Generating {len(pairs)} consecutive interferograms") + +for date1, date2 in pairs: + pair_dir = os.path.join(igram_dir, f"{date1}_{date2}") + os.makedirs(pair_dir, exist_ok=True) + + slc1 = os.path.join(coreg_dir, date1, date1 + ".slc") + slc2 = os.path.join(coreg_dir, date2, date2 + ".slc") + + out_prefix = os.path.join(pair_dir, f"{date1}_{date2}") + + cmd = [ + "python", os.path.join(script_dir, "crossmul.py"), + "-m", slc1, + "-s", slc2, + "-o", out_prefix, + "-a", "5", + "-r", "9", + ] + print(f" {date1} x {date2}") + subprocess.run(cmd, check=True) + +print("Interferogram generation complete") + +# %% [markdown] +# ## 7. Visualize results +# +# Show phase and coherence for each interferogram pair. + +# %% +def load_igram(pair_dir, pair_name): + """Load interferogram and amplitude from crossmul output.""" + int_file = os.path.join(pair_dir, pair_name + ".int") + amp_file = os.path.join(pair_dir, pair_name + ".amp") + + # Get dimensions from XML + import xml.etree.ElementTree as ET + xml_file = int_file + ".xml" + tree = ET.parse(xml_file) + root = tree.getroot() + + width = None + length = None + for prop in root.iter("property"): + name = prop.get("name") + if name == "width": + width = int(prop.find("value").text) + elif name == "length": + length = int(prop.find("value").text) + assert width and length, f"Could not read dimensions from {xml_file}" + + igram = np.fromfile(int_file, dtype=np.complex64).reshape(length, width) + amp = np.fromfile(amp_file, dtype=np.float32).reshape(length, width * 2) + + return igram, amp, width, length + + +# %% +fig, axes = plt.subplots(len(pairs), 2, figsize=(14, 4 * len(pairs))) +if len(pairs) == 1: + axes = axes[np.newaxis, :] + +for idx, (date1, date2) in enumerate(pairs): + pair_name = f"{date1}_{date2}" + pair_dir = os.path.join(igram_dir, pair_name) + + igram, amp, width, length = load_igram(pair_dir, pair_name) + + # Phase + phase = np.angle(igram) + axes[idx, 0].imshow(phase, cmap="hsv", vmin=-np.pi, vmax=np.pi, aspect="auto") + axes[idx, 0].set_title(f"Phase: {date1} - {date2}") + axes[idx, 0].set_ylabel("Azimuth") + + # Coherence (|E[s1*s2]| / sqrt(E[|s1|^2] * E[|s2|^2])) + # Simple magnitude-based estimate + coherence = np.abs(igram) / (amp[:, :width] * amp[:, width:] + 1e-10) + coherence = np.clip(coherence, 0, 1) + im = axes[idx, 1].imshow(coherence, cmap="gray", vmin=0, vmax=1, aspect="auto") + mean_coh = np.nanmean(coherence[coherence > 0]) + axes[idx, 1].set_title(f"Coherence: {date1} - {date2} (mean={mean_coh:.3f})") + +plt.tight_layout() +plt.savefig(os.path.join(work_dir, "interferograms_overview.png"), dpi=150, bbox_inches="tight") +plt.show() +print("Saved overview figure to interferograms_overview.png") + +# %% [markdown] +# ## 8. Compare with sarlet (optional) +# +# If sarlet/isce3 results exist for the same pairs, load and compare coherence. + +# %% +sarlet_dir = os.path.expanduser("~/repos/sarlet") +sarlet_output = os.path.join(sarlet_dir, "network_output2_mc") + +if os.path.exists(sarlet_output): + print("Found sarlet output, comparing coherence...") + # sarlet interferograms are named by date pair + for date1, date2 in pairs[:3]: # Compare first 3 pairs + sarlet_coh_file = glob.glob( + os.path.join(sarlet_output, f"*{date1}*{date2}*", "coherence.tif") + ) + if sarlet_coh_file: + from osgeo import gdal + ds = gdal.Open(sarlet_coh_file[0]) + sarlet_coh = ds.GetRasterBand(1).ReadAsArray() + ds = None + print(f" {date1}-{date2}: sarlet mean coherence = {np.nanmean(sarlet_coh):.3f}") +else: + print("No sarlet output found, skipping comparison") + +# %% [markdown] +# ## Notes +# +# - **Coherence**: Expect 0.35-0.55 for 3-day Capella stripmap pairs over urban areas +# - **Performance**: For a 1000x1000 crop with 8 dates, total processing is ~2-3 minutes +# - **refineSecondaryTiming**: This step runs ampcor for sub-pixel coregistration. +# The default grid is scaled to image size (max 10x10 windows). For very large scenes, +# processing time scales with the number of ampcor windows. +# - **Multi-looking**: The example uses 5 azimuth x 9 range looks (~5.8m x 5.6m pixels). +# Capella stripmap has azimuth spacing ~1.2m and range ~0.617m, so use ~1:2 (az:range) +# ratio for approximately square pixels. More looks = better coherence estimation +# but lower resolution. diff --git a/contrib/stack/stripmapStack/spotlightPhaseCorrection.py b/contrib/stack/stripmapStack/spotlightPhaseCorrection.py new file mode 100755 index 000000000..462bfabf7 --- /dev/null +++ b/contrib/stack/stripmapStack/spotlightPhaseCorrection.py @@ -0,0 +1,212 @@ +#!/usr/bin/env python3 + +""" +Apply phase correction to Capella spotlight SLCs for interferometry. + +Capella spotlight SLCs are delivered deramped (baseband). After coregistration +to a common reference geometry, a geometric phase must be restored for +interferometric processing. + +The correction phase for each pixel is: + phi = -4*pi/lambda * (|R_ant - target| - |R_ant - R_target|) + +where: + R_ant = reference antenna position (ECEF, from metadata) + R_target = scene reference target position (ECEF, from metadata) + target = each pixel's ECEF position (computed from lon/lat/hgt geometry) + +Each SLC (reference and all secondaries) gets its own correction using its +own antenna/target positions from the shelve data. +""" + +import argparse +import os +import shelve + +import numpy as np + +import isce +import isceobj +from isceobj.Planet.AstronomicalHandbook import Const + + +def createParser(): + parser = argparse.ArgumentParser( + description='Apply spotlight phase correction to coregistered Capella SLCs' + ) + parser.add_argument( + '-s', '--slc', type=str, required=True, + help='Path to SLC file (without .xml extension)' + ) + parser.add_argument( + '-g', '--geom', type=str, required=True, + help='Geometry directory containing lat.rdr, lon.rdr, hgt.rdr' + ) + parser.add_argument( + '-d', '--shelve', type=str, required=True, + help='Path to shelve data file containing frame with spotlight metadata' + ) + parser.add_argument( + '-o', '--output', type=str, default=None, + help='Output corrected SLC path (default: overwrite input)' + ) + parser.add_argument( + '--block-size', type=int, default=512, + help='Number of lines to process per block (default: 512)' + ) + return parser + + +def cmdLineParse(iargs=None): + return createParser().parse_args(args=iargs) + + +def llh_to_ecef(lon_deg, lat_deg, height): + """ + Convert geodetic lon/lat/height to ECEF XYZ coordinates (WGS84). + + Parameters + ---------- + lon_deg : ndarray, longitude in degrees + lat_deg : ndarray, latitude in degrees + height : ndarray, height above ellipsoid in meters + + Returns + ------- + x, y, z : ndarray, ECEF coordinates in meters + """ + # WGS84 parameters + a = 6378137.0 # semi-major axis + f = 1.0 / 298.257223563 # flattening + e2 = 2 * f - f * f # eccentricity squared + + lon = np.deg2rad(lon_deg) + lat = np.deg2rad(lat_deg) + + sin_lat = np.sin(lat) + cos_lat = np.cos(lat) + sin_lon = np.sin(lon) + cos_lon = np.cos(lon) + + N = a / np.sqrt(1.0 - e2 * sin_lat * sin_lat) + + x = (N + height) * cos_lat * cos_lon + y = (N + height) * cos_lat * sin_lon + z = (N * (1.0 - e2) + height) * sin_lat + + return x, y, z + + +def compute_correction_phase(lon, lat, hgt, ref_antenna_pos, ref_target_pos, wavelength): + """ + Compute the phase correction for spotlight re-ramping. + + phi = -4*pi/lambda * (|R_ant - pixel_ecef| - |R_ant - R_target|) + + Parameters + ---------- + lon, lat, hgt : ndarray (nlines, nsamples) - pixel coordinates + ref_antenna_pos : list [x, y, z] ECEF meters + ref_target_pos : list [x, y, z] ECEF meters + wavelength : float, meters + + Returns + ------- + phase : ndarray (nlines, nsamples), radians + """ + # Pixel ECEF positions + px, py, pz = llh_to_ecef(lon, lat, hgt) + + # Reference antenna position + ax, ay, az = ref_antenna_pos + + # Distance from antenna to each pixel + dx = px - ax + dy = py - ay + dz = pz - az + dist_ant_pixel = np.sqrt(dx * dx + dy * dy + dz * dz) + + # Distance from antenna to reference target (scalar) + tx, ty, tz = ref_target_pos + dist_ant_target = np.sqrt( + (ax - tx) ** 2 + (ay - ty) ** 2 + (az - tz) ** 2 + ) + + phase = (-4.0 * np.pi / wavelength) * (dist_ant_pixel - dist_ant_target) + return phase + + +def main(iargs=None): + inps = cmdLineParse(iargs) + + # Load frame from shelve to get spotlight metadata + with shelve.open(inps.shelve, flag='r') as db: + frame = db['frame'] + + ref_ant = frame.spotlightReferenceAntennaPosition + ref_tgt = frame.spotlightReferenceTargetPosition + wavelength = frame.getInstrument().getRadarWavelength() + + print(f'Wavelength: {wavelength:.6f} m') + print(f'Reference antenna position: {ref_ant}') + print(f'Reference target position: {ref_tgt}') + + # Load SLC image + slc_img = isceobj.createImage() + slc_img.load(inps.slc + '.xml') + width = slc_img.width + length = slc_img.length + + print(f'SLC dimensions: {length} x {width}') + + # Load geometry (DOUBLE precision) + lat_data = np.fromfile( + os.path.join(inps.geom, 'lat.rdr'), dtype=np.float64 + ).reshape(length, width) + lon_data = np.fromfile( + os.path.join(inps.geom, 'lon.rdr'), dtype=np.float64 + ).reshape(length, width) + hgt_data = np.fromfile( + os.path.join(inps.geom, 'hgt.rdr'), dtype=np.float64 + ).reshape(length, width) + + # Load SLC data + slc_data = np.fromfile(inps.slc, dtype=np.complex64).reshape(length, width) + + output_path = inps.output if inps.output else inps.slc + + # Process in blocks for memory efficiency + block_size = inps.block_size + for i0 in range(0, length, block_size): + i1 = min(i0 + block_size, length) + + phase = compute_correction_phase( + lon_data[i0:i1], + lat_data[i0:i1], + hgt_data[i0:i1], + ref_ant, ref_tgt, wavelength + ) + + # Apply phase correction: multiply by exp(-j*phi) + slc_data[i0:i1] *= np.exp(-1j * phase).astype(np.complex64) + + if (i0 // block_size) % 10 == 0: + print(f' Processed lines {i0}-{i1} / {length}') + + # Write corrected SLC + slc_data.tofile(output_path) + print(f'Wrote corrected SLC to: {output_path}') + + # Update XML if writing to a new file + if output_path != inps.slc: + out_img = isceobj.createSlcImage() + out_img.setByteOrder('l') + out_img.setFilename(output_path) + out_img.setAccessMode('read') + out_img.setWidth(width) + out_img.setLength(length) + out_img.renderHdr() + + +if __name__ == '__main__': + main() From e4d8782457a9341cd07571081b272f8e2d2b23b1 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Thu, 5 Mar 2026 14:13:17 -0500 Subject: [PATCH 14/17] Use TIFF dimensions instead of metadata rows/columns in Capella.py Cropped TIFFs (e.g. from sarlet crop-slc) correctly update timing and range metadata but may not update the rows/columns fields. Using the actual TIFF dimensions from GDAL ensures correct image sizing regardless of whether the metadata dimensions match. Co-Authored-By: Claude Opus 4.6 --- components/isceobj/Sensor/Capella.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/components/isceobj/Sensor/Capella.py b/components/isceobj/Sensor/Capella.py index 9ea0fe394..1a6dcf4be 100644 --- a/components/isceobj/Sensor/Capella.py +++ b/components/isceobj/Sensor/Capella.py @@ -79,7 +79,13 @@ def getFrame(self): return self.frame def _loadMetadataFromTiff(self): - """Load metadata from TIFF ImageDescription tag using GDAL.""" + """Load metadata from TIFF ImageDescription tag using GDAL. + + Uses the actual TIFF dimensions (from GDAL) for image size rather than + the metadata rows/columns, since properly cropped TIFFs (e.g. from + sarlet crop-slc) update the timing/range metadata but may not update + the rows/columns fields. + """ try: from osgeo import gdal except ImportError: @@ -91,6 +97,8 @@ def _loadMetadataFromTiff(self): # Get ImageDescription from TIFF metadata image_desc = ds.GetMetadataItem('TIFFTAG_IMAGEDESCRIPTION') + self._tiff_width = ds.RasterXSize + self._tiff_height = ds.RasterYSize ds = None if not image_desc: @@ -150,9 +158,9 @@ def populateMetadata(self): pulseBandwidth = 500e6 pulseLength = 10e-6 - # Image dimensions - lines = image.get('rows', image.get('length', 0)) - samples = image.get('columns', image.get('width', 0)) + # Image dimensions — use actual TIFF size (handles cropped TIFFs correctly) + lines = self._tiff_height + samples = self._tiff_width # Image geometry defines the actual SLC pixel grid # These are the ground truth for timing and range spacing From 064c3e012c73ba7043b658c6df9eb22f0aa8d7b1 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Tue, 10 Mar 2026 17:31:32 -0400 Subject: [PATCH 15/17] generalize --- contrib/stack/stripmapStack/example_capella_stripmap.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/contrib/stack/stripmapStack/example_capella_stripmap.py b/contrib/stack/stripmapStack/example_capella_stripmap.py index cf5329394..371f318e2 100755 --- a/contrib/stack/stripmapStack/example_capella_stripmap.py +++ b/contrib/stack/stripmapStack/example_capella_stripmap.py @@ -44,13 +44,13 @@ # %% # Input Capella SLC GeoTIFFs -slc_tiff_dir = "/Volumes/WD_BLACK_SN7100_4TB/Documents/Learning/capella/mexico_city/cropped_slcs" +slc_tiff_dir = "mexico_city/cropped_slcs" # DEM file (GeoTIFF - will be converted to ISCE format) -dem_tiff = "/Volumes/WD_BLACK_SN7100_4TB/Documents/Learning/capella/mexico_city/dem_3dep.tif" +dem_tiff = "mexico_city/dem_3dep.tif" # Working directory for ISCE2 processing -work_dir = os.path.expanduser("~/repos/isce2/data_notebook") +work_dir = os.path.expanduser("data_notebook") os.makedirs(work_dir, exist_ok=True) # Unpacked SLC output directory From 74fccf1439d24de9587105f57565116a47163b66 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Fri, 20 Mar 2026 11:33:02 -0400 Subject: [PATCH 16/17] Add options to set ampcor grid to speed up smaller images --- components/isceobj/Sensor/CMakeLists.txt | 1 + components/isceobj/Sensor/Capella.py | 13 ++-- .../stripmapStack/refineSecondaryTiming.py | 72 ++++++++++--------- .../stripmapStack/spotlightPhaseCorrection.py | 20 ++++-- 4 files changed, 61 insertions(+), 45 deletions(-) diff --git a/components/isceobj/Sensor/CMakeLists.txt b/components/isceobj/Sensor/CMakeLists.txt index 0348ecedc..404313662 100644 --- a/components/isceobj/Sensor/CMakeLists.txt +++ b/components/isceobj/Sensor/CMakeLists.txt @@ -42,6 +42,7 @@ set(installfiles UAVSAR_Stack.py SAOCOM_SLC.py Lutan1.py + Capella.py ) if(HDF5_FOUND) diff --git a/components/isceobj/Sensor/Capella.py b/components/isceobj/Sensor/Capella.py index 1a6dcf4be..f60c63625 100644 --- a/components/isceobj/Sensor/Capella.py +++ b/components/isceobj/Sensor/Capella.py @@ -249,16 +249,13 @@ def populateMetadata(self): # used to compute the geometric phase that must be restored after coregistration mode = collect.get('mode', '') if mode.lower() in {'spotlight', 'sp'}: - ref_ant = image.get('reference_antenna_position', {}) - ref_tgt = image.get('reference_target_position', {}) + ref_ant = image.get('reference_antenna_position') + ref_tgt = image.get('reference_target_position') assert ref_ant, "Spotlight mode requires reference_antenna_position in metadata" assert ref_tgt, "Spotlight mode requires reference_target_position in metadata" - self.frame.spotlightReferenceAntennaPosition = [ - ref_ant['x'], ref_ant['y'], ref_ant['z'] - ] - self.frame.spotlightReferenceTargetPosition = [ - ref_tgt['x'], ref_tgt['y'], ref_tgt['z'] - ] + # Positions are [x, y, z] ECEF coordinates in meters + self.frame.spotlightReferenceAntennaPosition = list(ref_ant) + self.frame.spotlightReferenceTargetPosition = list(ref_tgt) def _extractDopplerCoeffs(self, image_geometry, lines, samples): """ diff --git a/contrib/stack/stripmapStack/refineSecondaryTiming.py b/contrib/stack/stripmapStack/refineSecondaryTiming.py index fda3b0c9f..7566e168f 100755 --- a/contrib/stack/stripmapStack/refineSecondaryTiming.py +++ b/contrib/stack/stripmapStack/refineSecondaryTiming.py @@ -46,6 +46,19 @@ def createParser(): parser.add_argument('-t', '--thresh', dest='snrthresh', type=float, default=5.0, help='SNR threshold') + parser.add_argument('--nwa', dest='nwa', type=int, default=None, + help='Number of ampcor windows across (default: scale to image, max 10)') + parser.add_argument('--nwd', dest='nwd', type=int, default=None, + help='Number of ampcor windows down (default: scale to image, max 10)') + parser.add_argument('--ww', dest='windowWidth', type=int, default=64, + help='Ampcor window width (default: 64)') + parser.add_argument('--wh', dest='windowHeight', type=int, default=64, + help='Ampcor window height (default: 64)') + parser.add_argument('--sw', dest='searchWidth', type=int, default=20, + help='Ampcor search window half-width (default: 20)') + parser.add_argument('--sh', dest='searchHeight', type=int, default=20, + help='Ampcor search window half-height (default: 20)') + return parser def cmdLineParse(iargs = None): @@ -53,12 +66,14 @@ def cmdLineParse(iargs = None): return parser.parse_args(args=iargs) -def estimateOffsetField(reference, secondary, azoffset=0, rgoffset=0): +def estimateOffsetField(reference, secondary, azoffset=0, rgoffset=0, + windowWidth=64, windowHeight=64, + searchWidth=20, searchHeight=20, + nwa=None, nwd=None): ''' Estimate offset field between burst and simamp. ''' - sim = isceobj.createSlcImage() sim.load(secondary+'.xml') sim.setAccessMode('READ') @@ -76,48 +91,38 @@ def estimateOffsetField(reference, secondary, azoffset=0, rgoffset=0): objOffset.configure() objOffset.setAcrossGrossOffset(rgoffset) objOffset.setDownGrossOffset(azoffset) - objOffset.setWindowSizeWidth(128) - objOffset.setWindowSizeHeight(128) - objOffset.setSearchWindowSizeWidth(40) - objOffset.setSearchWindowSizeHeight(40) + objOffset.setWindowSizeWidth(windowWidth) + objOffset.setWindowSizeHeight(windowHeight) + objOffset.setSearchWindowSizeWidth(searchWidth) + objOffset.setSearchWindowSizeHeight(searchHeight) margin = 2*objOffset.searchWindowSizeWidth + objOffset.windowSizeWidth - nAcross = 60 - nDown = 60 + # Scale grid to image size (default: max 10 windows per dimension) + usableAc = int(min(width, sim.getWidth()) - 2 * margin) + usableDn = int(min(length, sim.getLength()) - 2 * margin) + nAcross = nwa if nwa is not None else min(10, max(3, usableAc // windowWidth)) + nDown = nwd if nwd is not None else min(10, max(3, usableDn // windowHeight)) - offAc = max(101,-rgoffset)+margin offDn = max(101,-azoffset)+margin - lastAc = int( min(width, sim.getWidth() - offAc) - margin) lastDn = int( min(length, sim.getLength() - offDn) - margin) -# print('Across: ', offAc, lastAc, width, sim.getWidth(), margin) -# print('Down: ', offDn, lastDn, length, sim.getLength(), margin) - - if not objOffset.firstSampleAcross: - objOffset.setFirstSampleAcross(offAc) - - if not objOffset.lastSampleAcross: - objOffset.setLastSampleAcross(lastAc) - - if not objOffset.firstSampleDown: - objOffset.setFirstSampleDown(offDn) - - if not objOffset.lastSampleDown: - objOffset.setLastSampleDown(lastDn) - - if not objOffset.numberLocationAcross: - objOffset.setNumberLocationAcross(nAcross) - - if not objOffset.numberLocationDown: - objOffset.setNumberLocationDown(nDown) + # Unconditionally set grid parameters (configure() defaults prevent conditional override) + objOffset.setFirstSampleAcross(offAc) + objOffset.setLastSampleAcross(lastAc) + objOffset.setFirstSampleDown(offDn) + objOffset.setLastSampleDown(lastDn) + objOffset.setNumberLocationAcross(nAcross) + objOffset.setNumberLocationDown(nDown) objOffset.setFirstPRF(1.0) objOffset.setSecondPRF(1.0) objOffset.setImageDataType1('complex') - objOffset.setImageDataType2('complex') + objOffset.setImageDataType2('complex') + + print(f' ampcor: {nAcross}x{nDown} grid, {windowWidth}x{windowHeight} window, {searchWidth}px search') objOffset.ampcor(sar, sim) @@ -174,7 +179,10 @@ def main(iargs=None): field = estimateOffsetField(inps.reference, inps.secondary, - azoffset=inps.azoff, rgoffset=inps.rgoff) + azoffset=inps.azoff, rgoffset=inps.rgoff, + windowWidth=inps.windowWidth, windowHeight=inps.windowHeight, + searchWidth=inps.searchWidth, searchHeight=inps.searchHeight, + nwa=inps.nwa, nwd=inps.nwd) if os.path.exists(inps.outfile): os.remove(inps.outfile) diff --git a/contrib/stack/stripmapStack/spotlightPhaseCorrection.py b/contrib/stack/stripmapStack/spotlightPhaseCorrection.py index 462bfabf7..c15366ffa 100755 --- a/contrib/stack/stripmapStack/spotlightPhaseCorrection.py +++ b/contrib/stack/stripmapStack/spotlightPhaseCorrection.py @@ -151,11 +151,21 @@ def main(iargs=None): print(f'Reference antenna position: {ref_ant}') print(f'Reference target position: {ref_tgt}') - # Load SLC image - slc_img = isceobj.createImage() - slc_img.load(inps.slc + '.xml') - width = slc_img.width - length = slc_img.length + # Load SLC image dimensions (prefer XML, fall back to VRT via GDAL) + xml_path = inps.slc + '.xml' + vrt_path = inps.slc + '.vrt' + if os.path.exists(xml_path): + slc_img = isceobj.createImage() + slc_img.load(xml_path) + width = slc_img.width + length = slc_img.length + else: + from osgeo import gdal + ds = gdal.Open(vrt_path, gdal.GA_ReadOnly) + assert ds is not None, f'Cannot open {vrt_path}' + width = ds.RasterXSize + length = ds.RasterYSize + ds = None print(f'SLC dimensions: {length} x {width}') From deb9e10721cf29a6bcfbc0e1bb42c78aa8e1a1d5 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Fri, 20 Mar 2026 11:35:35 -0400 Subject: [PATCH 17/17] cut til finished --- .../stripmapStack/example_capella_stripmap.py | 351 ------------------ 1 file changed, 351 deletions(-) delete mode 100755 contrib/stack/stripmapStack/example_capella_stripmap.py diff --git a/contrib/stack/stripmapStack/example_capella_stripmap.py b/contrib/stack/stripmapStack/example_capella_stripmap.py deleted file mode 100755 index 371f318e2..000000000 --- a/contrib/stack/stripmapStack/example_capella_stripmap.py +++ /dev/null @@ -1,351 +0,0 @@ -# --- -# jupyter: -# jupytext: -# formats: py:percent -# text_representation: -# extension: .py -# format_name: percent -# format_version: '1.3' -# kernelspec: -# display_name: Python 3 -# language: python -# name: python3 -# --- - -# %% [markdown] -# # Capella Stripmap InSAR with ISCE2's stripmapStack -# -# End-to-end workflow for processing Capella Space stripmap SLC GeoTIFFs -# through ISCE2's stripmapStack processor to generate interferograms. -# -# ## Prerequisites -# - ISCE2 installed with Capella sensor support -# - Capella SLC GeoTIFFs (stripmap mode) -# - A DEM covering the area of interest (e.g. Copernicus DEM or 3DEP) -# - stripmapStack scripts on PATH (or run from source tree) - -# %% [markdown] -# ## 1. Setup - -# %% -import os -import glob -import subprocess -import time -from pathlib import Path - -import numpy as np -import matplotlib.pyplot as plt - -# %% [markdown] -# ### Configure paths -# -# Update these paths to point to your data. - -# %% -# Input Capella SLC GeoTIFFs -slc_tiff_dir = "mexico_city/cropped_slcs" - -# DEM file (GeoTIFF - will be converted to ISCE format) -dem_tiff = "mexico_city/dem_3dep.tif" - -# Working directory for ISCE2 processing -work_dir = os.path.expanduser("data_notebook") -os.makedirs(work_dir, exist_ok=True) - -# Unpacked SLC output directory -slc_dir = os.path.join(work_dir, "SLC") -os.makedirs(slc_dir, exist_ok=True) - -# Path to stripmapStack scripts (if not on PATH, set this to source tree) -script_dir = os.path.expanduser("~/repos/isce2/contrib/stack/stripmapStack") - -# %% -# List available SLC TIFFs -tiffs = sorted(glob.glob(os.path.join(slc_tiff_dir, "CAPELLA*.tif"))) -print(f"Found {len(tiffs)} Capella SLC TIFFs:") -for t in tiffs: - print(f" {os.path.basename(t)}") - -# %% [markdown] -# ## 2. Prepare DEM -# -# Convert the GeoTIFF DEM to ISCE2 format using `dem_to_isce.py` or GDAL. -# The DEM must have `.xml` and `.vrt` sidecar files. - -# %% -dem_isce = os.path.join(work_dir, "dem.wgs84") - -if not os.path.exists(dem_isce + ".xml"): - subprocess.run( - ["gdal_translate", "-of", "ENVI", dem_tiff, dem_isce], - check=True - ) - # Create ISCE XML/VRT metadata - subprocess.run( - ["fixImageXml.py", "-f", "-i", dem_isce], - check=True - ) - print("DEM converted to ISCE format") -else: - print("DEM already exists in ISCE format") - -# %% [markdown] -# ## 3. Organize SLCs by date and unpack -# -# `prepSlcCapella.py` reads each TIFF's metadata to extract the acquisition -# date, moves files into date folders, and generates an unpacking run file. -# -# `unpackFrame_Capella.py` converts each GeoTIFF to ISCE2's binary SLC format -# and extracts orbit/timing metadata into a shelve file. - -# %% -# Step 1: Organize TIFFs into date folders and generate unpack commands -os.chdir(work_dir) -subprocess.run( - [ - "python", os.path.join(script_dir, "prepSlcCapella.py"), - "-i", slc_tiff_dir, - "-o", slc_dir, - ], - check=True -) - -# %% -# Step 2: Run the unpack commands -run_file = os.path.join(work_dir, "run_unPackCapella") -if os.path.exists(run_file): - with open(run_file) as f: - cmds = f.read().strip().split("\n") - - print(f"Unpacking {len(cmds)} SLCs...") - t0 = time.time() - for cmd in cmds: - print(f" {cmd}") - subprocess.run(cmd.split(), check=True) - print(f"Unpacking done in {time.time() - t0:.1f}s") - -# %% -# Verify unpacked SLCs -dates = sorted([d for d in os.listdir(slc_dir) if os.path.isdir(os.path.join(slc_dir, d))]) -print(f"\nUnpacked {len(dates)} dates: {dates}") - -# %% [markdown] -# ## 4. Generate stripmapStack run files -# -# `stackStripMap.py` creates the coregistration processing steps (run files). -# We use `--nofocus` since Capella data is already focused to SLC. - -# %% -os.chdir(work_dir) -subprocess.run( - [ - "python", os.path.join(script_dir, "stackStripMap.py"), - "-s", slc_dir, - "-d", dem_isce, - "-w", work_dir, - "-a", "1", - "-r", "1", - "--nofocus", - "-W", "slc", - ], - check=True -) - -# %% -# List generated run files -run_files = sorted(glob.glob(os.path.join(work_dir, "run_files", "run_*"))) -print("Generated run files:") -for rf in run_files: - print(f" {os.path.basename(rf)}") - -# %% [markdown] -# ## 5. Execute coregistration steps -# -# Run each step sequentially. Steps: -# 1. Reference geometry (topo) -# 2. Secondary focus/split (no-op for pre-focused SLCs) -# 3. Geo2rdr + coarse resampling -# 4. Refine secondary timing (ampcor) - **slowest step** -# 5. Invert misregistration -# 6. Fine resampling -# 7. Grid baseline -# -# **Note on step 4 (refineSecondaryTiming):** The default 60x60 ampcor grid -# can be very slow. For small test scenes, this is fine. For large scenes, -# you can reduce the grid by editing the config files to add `--nwa 10 --nwd 10`. - -# %% -os.chdir(work_dir) -total_t0 = time.time() - -for rf in run_files: - step_name = os.path.basename(rf) - print(f"\n{'='*60}") - print(f"Running: {step_name}") - print(f"{'='*60}") - - with open(rf) as f: - cmds = f.read().strip().split("\n") - - t0 = time.time() - for cmd in cmds: - if cmd.strip(): - subprocess.run(cmd, shell=True, check=True) - elapsed = time.time() - t0 - print(f" {step_name} completed in {elapsed:.1f}s") - -total_elapsed = time.time() - total_t0 -print(f"\nAll steps completed in {total_elapsed:.1f}s") - -# %% [markdown] -# ## 6. Generate interferograms -# -# Create interferograms for consecutive date pairs using `crossmul.py`. -# -# **Multi-looking ratio**: Capella stripmap has range pixel size ~0.617m and -# azimuth spacing ~1.2m. For approximately square pixels, use a 1:2 ratio -# (az:range looks). Example: 5 azimuth x 9 range gives ~5.8m x 5.6m pixels. - -# %% -# Find all coregistered SLC dates -merged_dir = os.path.join(work_dir, "merged") -coreg_dir = os.path.join(merged_dir, "SLC") - -# Reference date -ref_date = dates[0] -ref_slc = os.path.join(coreg_dir, ref_date, ref_date + ".slc") - -all_dates = sorted(os.listdir(coreg_dir)) -print(f"Coregistered dates: {all_dates}") - -# %% -# Generate interferograms for consecutive pairs -igram_dir = os.path.join(work_dir, "interferograms") -os.makedirs(igram_dir, exist_ok=True) - -pairs = list(zip(all_dates[:-1], all_dates[1:])) -print(f"Generating {len(pairs)} consecutive interferograms") - -for date1, date2 in pairs: - pair_dir = os.path.join(igram_dir, f"{date1}_{date2}") - os.makedirs(pair_dir, exist_ok=True) - - slc1 = os.path.join(coreg_dir, date1, date1 + ".slc") - slc2 = os.path.join(coreg_dir, date2, date2 + ".slc") - - out_prefix = os.path.join(pair_dir, f"{date1}_{date2}") - - cmd = [ - "python", os.path.join(script_dir, "crossmul.py"), - "-m", slc1, - "-s", slc2, - "-o", out_prefix, - "-a", "5", - "-r", "9", - ] - print(f" {date1} x {date2}") - subprocess.run(cmd, check=True) - -print("Interferogram generation complete") - -# %% [markdown] -# ## 7. Visualize results -# -# Show phase and coherence for each interferogram pair. - -# %% -def load_igram(pair_dir, pair_name): - """Load interferogram and amplitude from crossmul output.""" - int_file = os.path.join(pair_dir, pair_name + ".int") - amp_file = os.path.join(pair_dir, pair_name + ".amp") - - # Get dimensions from XML - import xml.etree.ElementTree as ET - xml_file = int_file + ".xml" - tree = ET.parse(xml_file) - root = tree.getroot() - - width = None - length = None - for prop in root.iter("property"): - name = prop.get("name") - if name == "width": - width = int(prop.find("value").text) - elif name == "length": - length = int(prop.find("value").text) - assert width and length, f"Could not read dimensions from {xml_file}" - - igram = np.fromfile(int_file, dtype=np.complex64).reshape(length, width) - amp = np.fromfile(amp_file, dtype=np.float32).reshape(length, width * 2) - - return igram, amp, width, length - - -# %% -fig, axes = plt.subplots(len(pairs), 2, figsize=(14, 4 * len(pairs))) -if len(pairs) == 1: - axes = axes[np.newaxis, :] - -for idx, (date1, date2) in enumerate(pairs): - pair_name = f"{date1}_{date2}" - pair_dir = os.path.join(igram_dir, pair_name) - - igram, amp, width, length = load_igram(pair_dir, pair_name) - - # Phase - phase = np.angle(igram) - axes[idx, 0].imshow(phase, cmap="hsv", vmin=-np.pi, vmax=np.pi, aspect="auto") - axes[idx, 0].set_title(f"Phase: {date1} - {date2}") - axes[idx, 0].set_ylabel("Azimuth") - - # Coherence (|E[s1*s2]| / sqrt(E[|s1|^2] * E[|s2|^2])) - # Simple magnitude-based estimate - coherence = np.abs(igram) / (amp[:, :width] * amp[:, width:] + 1e-10) - coherence = np.clip(coherence, 0, 1) - im = axes[idx, 1].imshow(coherence, cmap="gray", vmin=0, vmax=1, aspect="auto") - mean_coh = np.nanmean(coherence[coherence > 0]) - axes[idx, 1].set_title(f"Coherence: {date1} - {date2} (mean={mean_coh:.3f})") - -plt.tight_layout() -plt.savefig(os.path.join(work_dir, "interferograms_overview.png"), dpi=150, bbox_inches="tight") -plt.show() -print("Saved overview figure to interferograms_overview.png") - -# %% [markdown] -# ## 8. Compare with sarlet (optional) -# -# If sarlet/isce3 results exist for the same pairs, load and compare coherence. - -# %% -sarlet_dir = os.path.expanduser("~/repos/sarlet") -sarlet_output = os.path.join(sarlet_dir, "network_output2_mc") - -if os.path.exists(sarlet_output): - print("Found sarlet output, comparing coherence...") - # sarlet interferograms are named by date pair - for date1, date2 in pairs[:3]: # Compare first 3 pairs - sarlet_coh_file = glob.glob( - os.path.join(sarlet_output, f"*{date1}*{date2}*", "coherence.tif") - ) - if sarlet_coh_file: - from osgeo import gdal - ds = gdal.Open(sarlet_coh_file[0]) - sarlet_coh = ds.GetRasterBand(1).ReadAsArray() - ds = None - print(f" {date1}-{date2}: sarlet mean coherence = {np.nanmean(sarlet_coh):.3f}") -else: - print("No sarlet output found, skipping comparison") - -# %% [markdown] -# ## Notes -# -# - **Coherence**: Expect 0.35-0.55 for 3-day Capella stripmap pairs over urban areas -# - **Performance**: For a 1000x1000 crop with 8 dates, total processing is ~2-3 minutes -# - **refineSecondaryTiming**: This step runs ampcor for sub-pixel coregistration. -# The default grid is scaled to image size (max 10x10 windows). For very large scenes, -# processing time scales with the number of ampcor windows. -# - **Multi-looking**: The example uses 5 azimuth x 9 range looks (~5.8m x 5.6m pixels). -# Capella stripmap has azimuth spacing ~1.2m and range ~0.617m, so use ~1:2 (az:range) -# ratio for approximately square pixels. More looks = better coherence estimation -# but lower resolution.