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}") 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/ 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 new file mode 100644 index 000000000..f60c63625 --- /dev/null +++ b/components/isceobj/Sensor/Capella.py @@ -0,0 +1,439 @@ +#!/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 numpy as np + +import isceobj +from isceobj.Scene.Frame import Frame +from isceobj.Planet.Planet import Planet +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', 'spotlight', 'SP'} + +TIFF = Component.Parameter( + 'tiff', + public_name='TIFF', + default='', + type=str, + mandatory=True, + doc='Capella SLC GeoTIFF imagery file' +) + + +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,) + 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 + + def getFrame(self): + return self.frame + + def _loadMetadataFromTiff(self): + """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: + 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') + self._tiff_width = ds.RasterXSize + self._tiff_height = ds.RasterYSize + 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.""" + 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"Supported modes: {SUPPORTED_MODES}. " + f"Sliding Spotlight (SL) mode is not yet supported." + ) + + def populateMetadata(self): + """Create metadata objects from the JSON metadata.""" + + collect = self._metadata.get('collect', {}) + # 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 + 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', 'V') + rxPol = radar.get('receive_polarization', 'V') + polarization = txPol + rxPol + + # 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 + + # 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 + image_geometry = image.get('image_geometry', {}) + 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: 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: + 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' + + # 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) + sensingMid = dataStartTime + datetime.timedelta( + seconds=(lines - 1) * delta_line_time / 2.0 + ) + 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) + + # 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) + + # 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" + # 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): + """ + 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 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 = dc_poly.get('coefficients', [[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 + mid_az_pixel = lines / 2.0 + + # 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: + coeffs_2d = np.array(coeffs, dtype=np.float64) + degree_az = coeffs_2d.shape[0] - 1 + degree_rg = coeffs_2d.shape[1] - 1 + + # 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] * (mid_az_pixel ** 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.""" + 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)) + + # 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) + + 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 + 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/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) 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 100755 index 000000000..5246f535f --- /dev/null +++ b/contrib/stack/stripmapStack/prepSlcCapella.py @@ -0,0 +1,141 @@ +#!/usr/bin/env python3 + +import os +import glob +import argparse +import json + + +def createParser(): + ''' + Create command line parser. + ''' + + 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 TIFF files') + parser.add_argument('-o', '--output', dest='output', type=str, required=False, + 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 + + +def cmdLineParse(iargs=None): + ''' + Command line parser. + ''' + + parser = createParser() + return parser.parse_args(args=iargs) + + +def get_date_from_tiff(tiffFile): + ''' + 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}') + + return False, 'FAIL' + + +def main(iargs=None): + ''' + The main driver. + ''' + + inps = cmdLineParse(iargs) + inputDir = os.path.abspath(inps.input) + + if inps.output: + outputDir = os.path.abspath(inps.output) + else: + outputDir = None + + # filename of the runfile + run_unPack = 'run_unPackCapella' + + # 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')) + + # Organize TIFF files into date folders + for tiffFile in tiffFiles: + successflag, imgDate = get_date_from_tiff(tiffFile) + + 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(f'Failed to get date from: {tiffFile}') + + # Generate the unpacking script for all date dirs + dateDirs = sorted(glob.glob(os.path.join(inputDir, '2*'))) + + if outputDir is not None: + 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__': + + main() diff --git a/contrib/stack/stripmapStack/prepSlcSensors.py b/contrib/stack/stripmapStack/prepSlcSensors.py index 7f81b46c0..45e456e0b 100755 --- a/contrib/stack/stripmapStack/prepSlcSensors.py +++ b/contrib/stack/stripmapStack/prepSlcSensors.py @@ -65,11 +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*.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) - 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) + 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/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 new file mode 100755 index 000000000..c15366ffa --- /dev/null +++ b/contrib/stack/stripmapStack/spotlightPhaseCorrection.py @@ -0,0 +1,222 @@ +#!/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 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}') + + # 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() diff --git a/contrib/stack/stripmapStack/unpackFrame_Capella.py b/contrib/stack/stripmapStack/unpackFrame_Capella.py new file mode 100755 index 000000000..e4ce2ea7b --- /dev/null +++ b/contrib/stack/stripmapStack/unpackFrame_Capella.py @@ -0,0 +1,75 @@ +#!/usr/bin/env python3 + +import isce +from isceobj.Sensor import createSensor +import shelve +import argparse +import glob +from isceobj.Util import Poly1D +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='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') + + return parser.parse_args() + + +@use_api +def unpack(tiffFile, slcname): + ''' + Unpack Capella data to binary SLC file. + ''' + + # 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.tiff = tiffFile + 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 (Capella provides doppler_centroid_polynomial) + coeffs = obj.doppler_coeff + poly = Poly1D.Poly1D() + poly.initPoly(order=len(coeffs) - 1) + poly.setCoeffs(coeffs) + + # Save required metadata for further use + # 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 + + +if __name__ == '__main__': + ''' + Main driver. + ''' + + inps = cmdLineParse() + if inps.slcdir.endswith('/'): + inps.slcdir = inps.slcdir[:-1] + + if inps.tiffFile.endswith('/'): + inps.tiffFile = inps.tiffFile[:-1] + + unpack(inps.tiffFile, inps.slcdir)