From c027e1690bb502137be2da82ea5df69efe1896f5 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Wed, 21 May 2025 07:12:32 -0500 Subject: [PATCH 01/12] cleanup --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index bbca664..97650d9 100644 --- a/.gitignore +++ b/.gitignore @@ -240,3 +240,5 @@ Sessionx.vim tags # Persistent undo [._]*.un~ +# Data +ale From 78eea4af4684f28e465622b4e21bd27ffe770e01 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Wed, 21 May 2025 15:30:22 -0500 Subject: [PATCH 02/12] first kinda-close capella geocode --- src/multirtc/create_rtc.py | 23 ++--- src/multirtc/multirtc.py | 41 ++++++++- src/multirtc/prep_capella.py | 164 +++++++++++++++++++++++++++++++++++ 3 files changed, 211 insertions(+), 17 deletions(-) create mode 100644 src/multirtc/prep_capella.py diff --git a/src/multirtc/create_rtc.py b/src/multirtc/create_rtc.py index 862002e..a2508ed 100755 --- a/src/multirtc/create_rtc.py +++ b/src/multirtc/create_rtc.py @@ -32,7 +32,6 @@ def compute_layover_shadow_mask( radar_grid: isce3.product.RadarGridParameters, orbit: isce3.core.Orbit, geogrid_in: isce3.product.GeoGridParameters, - slc_obj: Sentinel1BurstSlc, dem_raster: isce3.io.Raster, filename_out: str, output_raster_format: str, @@ -59,8 +58,6 @@ def compute_layover_shadow_mask( Orbit defining radar motion on input path geogrid_in: isce3.product.GeoGridParameters Geogrid to geocode the layover/shadow mask in radar grid - burst_in: Sentinel1BurstSlc - Input burst geogrid_in: isce3.product.GeoGridParameters Geogrid to geocode the layover/shadow mask in radar grid dem_raster: isce3.io.Raster @@ -99,9 +96,6 @@ def compute_layover_shadow_mask( if doppler is None: doppler = isce3.core.LUT2d() - # determine the output filename - str_datetime = slc_obj.sensing_start.strftime('%Y%m%d_%H%M%S.%f') - # Run topo to get layover/shadow mask ellipsoid = isce3.core.Ellipsoid() grid_doppler = doppler @@ -122,9 +116,8 @@ def compute_layover_shadow_mask( path_layover_shadow_mask_file, radar_grid.width, radar_grid.length, 1, gdal.GDT_Byte, 'GTiff' ) else: - path_layover_shadow_mask = f'layover_shadow_mask_{str_datetime}' slantrange_layover_shadow_mask_raster = isce3.io.Raster( - path_layover_shadow_mask, radar_grid.width, radar_grid.length, 1, gdal.GDT_Byte, 'MEM' + 'layover_shadow_mask', radar_grid.width, radar_grid.length, 1, gdal.GDT_Byte, 'MEM' ) rdr2geo_obj.topo(dem_raster, layover_shadow_raster=slantrange_layover_shadow_mask_raster) @@ -416,7 +409,6 @@ def run_single_job(product_id: str, burst: Sentinel1BurstSlc, geogrid, opts: Rtc radar_grid_layover_shadow_mask, orbit, geogrid, - burst, dem_raster, layover_shadow_mask_file, raster_format, @@ -558,11 +550,11 @@ def run_single_job(product_id: str, burst: Sentinel1BurstSlc, geogrid, opts: Rtc logger.info(f'elapsed time: {t_end - t_start}') -def umbra_rtc_with_radargrid(umbra_sicd, geogrid, opts): +def capella_rtc(capella_sicd, geogrid, opts): # Common initializations t_start = time.time() output_dir = str(opts.output_dir) - product_id = umbra_sicd.id + product_id = capella_sicd.id os.makedirs(output_dir, exist_ok=True) raster_format = 'GTiff' @@ -575,14 +567,14 @@ def umbra_rtc_with_radargrid(umbra_sicd, geogrid, opts): rtc_anf_gamma0_to_sigma0_file = ( f'{output_dir}/{product_id}_{LAYER_NAME_RTC_ANF_GAMMA0_TO_SIGMA0}.{raster_extension}' ) - radar_grid = umbra_sicd.as_isce3_radargrid() - orbit = umbra_sicd.orbit - wavelength = umbra_sicd.wavelength + radar_grid = capella_sicd.as_isce3_radargrid() + orbit = capella_sicd.orbit + wavelength = capella_sicd.wavelength lookside = radar_grid.lookside dem_raster = isce3.io.Raster(opts.dem_path) ellipsoid = isce3.core.Ellipsoid() - doppler = umbra_sicd.get_doppler_centroid_grid() + doppler = capella_sicd.get_doppler_centroid_grid() exponent = 2 x_snap = geogrid.spacing_x @@ -605,7 +597,6 @@ def umbra_rtc_with_radargrid(umbra_sicd, geogrid, opts): radar_grid_layover_shadow_mask, orbit, geogrid, - umbra_sicd, dem_raster, layover_shadow_mask_file, raster_format, diff --git a/src/multirtc/multirtc.py b/src/multirtc/multirtc.py index 584a442..9ed5d98 100644 --- a/src/multirtc/multirtc.py +++ b/src/multirtc/multirtc.py @@ -2,13 +2,27 @@ from pathlib import Path from typing import Optional +import isce3 + from multirtc.create_rtc import run_single_job, umbra_rtc from multirtc.define_geogrid import generate_geogrids from multirtc.prep_burst import prep_burst +from multirtc.prep_capella import prep_capella from multirtc.prep_umbra import prep_umbra from multirtc.rtc_options import RtcOptions +def print_wkt(sicd): + radar_grid = sicd.as_isce3_radar_grid() + dem = isce3.geometry.DEMInterpolator() + doppler = sicd.get_doppler_centroid_grid() + breakpoint() + wkt = isce3.geometry.get_geo_perimeter_wkt( + grid=radar_grid, orbit=sicd.orbit, doppler=doppler, dem=dem, points_per_edge=3 + ) + print(wkt) + + def opera_rtc_s1_burst(granule: str, resolution: int = 30, work_dir: Optional[Path] = None) -> None: """Create an OPERA RTC for a Sentinel-1 burst @@ -50,6 +64,29 @@ def opera_rtc_umbra_sicd(granule: str, resolution: int = 30, work_dir: Optional[ umbra_rtc(umbra_sicd, geogrid, dem_path, output_dir=output_dir) +def opera_rtc_capella_sicd(granule: str, resolution: int = 30, work_dir: Optional[Path] = None) -> None: + """Create an OPERA RTC for an CAPELLA SICD file + + Args: + granule: Capella SICD file name to create an RTC for + resolution: Resolution of the output RTC (m) + work_dir: Working directory for processing + """ + if work_dir is None: + work_dir = Path.cwd() + input_dir = work_dir / 'input' + output_dir = work_dir / 'output' + granule_path = input_dir / granule + if not granule_path.exists(): + raise FileNotFoundError(f'Capella SICD must be present in input dir {input_dir} for processing.') + [d.mkdir(parents=True, exist_ok=True) for d in [input_dir, output_dir]] + capella_sicd, dem_path = prep_capella(granule_path, work_dir=input_dir) + print_wkt(capella_sicd) + # opts = RtcOptions(dem_path=str(dem_path), output_dir=str(output_dir), resolution=resolution) + # geogrid = generate_geogrids(burst, opts.resolution) + # run_single_job(granule, burst, geogrid, opts) + + def main(): """Create an OPERA RTC for an Umbra SICD SLC granule @@ -57,7 +94,7 @@ def main(): multirtc umbra_image.ntif --resolution 40 """ parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.ArgumentDefaultsHelpFormatter) - parser.add_argument('platform', choices=['S1', 'UMBRA'], help='Platform to create RTC for') + parser.add_argument('platform', choices=['S1', 'UMBRA', 'CAPELLA'], help='Platform to create RTC for') parser.add_argument('granule', help='Data granule to create an RTC for.') parser.add_argument('--resolution', default=30, type=float, help='Resolution of the output RTC (m)') parser.add_argument('--work-dir', type=Path, default=None, help='Working directory for processing') @@ -67,6 +104,8 @@ def main(): opera_rtc_s1_burst(args.granule, args.resolution, args.work_dir) elif args.platform == 'UMBRA': opera_rtc_umbra_sicd(args.granule, args.resolution, args.work_dir) + elif args.platform == 'CAPELLA': + opera_rtc_capella_sicd(args.granule, args.resolution, args.work_dir) else: raise NotImplementedError('Only Sentinel-1 burst and Umbra processing are supported at this time') diff --git a/src/multirtc/prep_capella.py b/src/multirtc/prep_capella.py new file mode 100644 index 0000000..4adfcd0 --- /dev/null +++ b/src/multirtc/prep_capella.py @@ -0,0 +1,164 @@ +from dataclasses import dataclass +from datetime import datetime, timedelta +from pathlib import Path +from typing import Optional + +import isce3 +import numpy as np +from numpy.polynomial.polynomial import polyder, polyval, polyval2d +from sarpy.io.complex.sicd import SICDReader +from shapely.geometry import Polygon + +from multirtc import dem +from multirtc.define_geogrid import get_point_epsg + + +def check_poly_order(poly): + assert len(poly.Coefs) == poly.order1 + 1, 'Polynomial order does not match number of coefficients' + + +def to_isce_datetime(dt): + if isinstance(dt, datetime): + return isce3.core.DateTime(dt) + elif isinstance(dt, np.datetime64): + return isce3.core.DateTime(dt.item()) + else: + raise ValueError(f'Unsupported datetime type: {type(dt)}. Expected datetime or np.datetime64.') + + +@dataclass +class Point: + x: float + y: float + + +@dataclass +class CapellaSICD: + id: str + file_path: Path + footprint: Polygon + wavelength: float + polarization: str + lookside: str # 'right' or 'left' + shape: tuple[int, int] + scp_index: tuple[int, int] + range_pixel_spacing: float + collect_start: datetime + sensing_start: float + sensing_end: float + prf: float + starting_range: float + orbit: isce3.core.Orbit + beta0_coeff: np.ndarray + + @staticmethod + def get_arp_pos_at_time(time, arp_pos_poly): + x = polyval(time, arp_pos_poly.X.get_array()) + y = polyval(time, arp_pos_poly.Y.get_array()) + z = polyval(time, arp_pos_poly.Z.get_array()) + return np.array([x, y, z]) + + @staticmethod + def get_arp_vel_at_time(time, arp_pos_poly): + x = polyval(time, polyder(arp_pos_poly.X.get_array())) + y = polyval(time, polyder(arp_pos_poly.Y.get_array())) + z = polyval(time, polyder(arp_pos_poly.Z.get_array())) + return np.array([x, y, z]) + + @staticmethod + def calculate_orbit(epoch, sensing_start, sensing_end, arp_pos_poly): + svs = [] + orbit_start = np.floor(sensing_start) - 5 + orbit_end = np.ceil(sensing_end) + 5 + for offset_sec in np.arange(orbit_start, orbit_end + 1, 1): + t = sensing_start + offset_sec + pos = CapellaSICD.get_arp_pos_at_time(t, arp_pos_poly) + vel = CapellaSICD.get_arp_vel_at_time(t, arp_pos_poly) + t_isce = to_isce_datetime(epoch + timedelta(seconds=t)) + svs.append(isce3.core.StateVector(t_isce, pos, vel)) + return isce3.core.Orbit(svs, to_isce_datetime(epoch)) + + def as_isce3_radar_grid(self): + radar_grid = isce3.product.RadarGridParameters( + sensing_start=self.sensing_start, + wavelength=self.wavelength, + prf=self.prf, + starting_range=self.starting_range, + range_pixel_spacing=self.range_pixel_spacing, + lookside=isce3.core.LookSide.Right if self.lookside == 'right' else isce3.core.LookSide.Left, + length=self.shape[1], # flipped for "shadows down" convention + width=self.shape[0], # flipped for "shadows down" convention + ref_epoch=to_isce_datetime(self.collect_start), + ) + return radar_grid + + def get_doppler_centroid_grid(self): + return isce3.core.LUT2d() + + @classmethod + def from_sarpy_sicd(cls, sicd, file_path): + center_frequency = sicd.RadarCollection.TxFrequency.Min + sicd.RadarCollection.TxFrequency.Max / 2 + wavelength = isce3.core.speed_of_light / center_frequency + polarization = sicd.RadarCollection.RcvChannels[0].TxRcvPolarization.replace(':', '') + lookside = 'right' if sicd.SCPCOA.SideOfTrack == 'R' else 'left' + footprint = Polygon([(ic.Lon, ic.Lat) for ic in sicd.GeoData.ImageCorners]) + shape = np.array([sicd.ImageData.NumRows, sicd.ImageData.NumCols]) + scp_index = np.array([sicd.ImageData.SCPPixel.Row, sicd.ImageData.SCPPixel.Col]) + start_index = scp_index - shape + r_ss = sicd.Grid.Row.SS + c_ss = sicd.Grid.Col.SS + range_pixel_spacing = r_ss + assert sicd.Grid.Type == 'RGZERO', 'Only range zero doppler grids supported for Capella data' + collect_start = sicd.Timeline.CollectStart + # seconds after collect start + tcoa_coefs = sicd.Grid.TimeCOAPoly.Coefs + first_col_time = polyval2d(start_index[0] * r_ss, start_index[1] * c_ss, tcoa_coefs) + last_col_time = polyval2d(start_index[0] * r_ss, (start_index[1] + shape[1]) * c_ss, tcoa_coefs) + sensing_start = min(first_col_time, last_col_time) + sensing_end = max(first_col_time, last_col_time) + prf = np.mean([ipp.IPPPoly.Coefs[1] for ipp in sicd.Timeline.IPP]) + starting_grid_pos = ( + sicd.Grid.Row.UVectECF.get_array() * r_ss * start_index[0] + + sicd.Grid.Col.UVectECF.get_array() * c_ss * start_index[1] + + sicd.GeoData.SCP.ECF.get_array() + ) + starting_arp_pos = cls.get_arp_pos_at_time(sensing_start, sicd.Position.ARPPoly) + starting_range = np.linalg.norm(starting_grid_pos - starting_arp_pos) + orbit = CapellaSICD.calculate_orbit(collect_start, sensing_start, sensing_end, sicd.Position.ARPPoly) + beta0_coeff = sicd.Radiometric.BetaZeroSFPoly.Coefs + capella_sicd = cls( + id=Path(file_path).with_suffix('').name, + file_path=file_path, + footprint=footprint, + shape=shape, + wavelength=wavelength, + polarization=polarization, + lookside=lookside, + range_pixel_spacing=range_pixel_spacing, + scp_index=scp_index, + collect_start=collect_start, + sensing_start=sensing_start, + sensing_end=sensing_end, + prf=prf, + starting_range=starting_range, + orbit=orbit, + beta0_coeff=beta0_coeff, + ) + return capella_sicd + + +def prep_capella(granule_path: Path, work_dir: Optional[Path] = None) -> Path: + """Prepare data for burst-based processing. + + Args: + granule_path: Path to the UMBRA SICD file + work_dir: Working directory for processing + """ + if work_dir is None: + work_dir = Path.cwd() + reader = SICDReader(str(granule_path)) + sicd = reader.get_sicds_as_tuple()[0] + capella_sicd = CapellaSICD.from_sarpy_sicd(sicd, granule_path) + dem_path = work_dir / 'dem.tif' + dem.download_opera_dem_for_footprint(dem_path, capella_sicd.footprint) + return capella_sicd, dem_path From 4859756a03f8e2585e0048cefe8d6132b03e3cad Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Thu, 22 May 2025 08:12:04 -0500 Subject: [PATCH 03/12] corrected range --- src/multirtc/multirtc.py | 3 +-- src/multirtc/prep_capella.py | 49 ++++++++++++------------------------ 2 files changed, 17 insertions(+), 35 deletions(-) diff --git a/src/multirtc/multirtc.py b/src/multirtc/multirtc.py index 9ed5d98..d478315 100644 --- a/src/multirtc/multirtc.py +++ b/src/multirtc/multirtc.py @@ -14,9 +14,8 @@ def print_wkt(sicd): radar_grid = sicd.as_isce3_radar_grid() - dem = isce3.geometry.DEMInterpolator() + dem = isce3.geometry.DEMInterpolator(sicd.hae) doppler = sicd.get_doppler_centroid_grid() - breakpoint() wkt = isce3.geometry.get_geo_perimeter_wkt( grid=radar_grid, orbit=sicd.orbit, doppler=doppler, dem=dem, points_per_edge=3 ) diff --git a/src/multirtc/prep_capella.py b/src/multirtc/prep_capella.py index 4adfcd0..dae5626 100644 --- a/src/multirtc/prep_capella.py +++ b/src/multirtc/prep_capella.py @@ -5,12 +5,10 @@ import isce3 import numpy as np -from numpy.polynomial.polynomial import polyder, polyval, polyval2d from sarpy.io.complex.sicd import SICDReader from shapely.geometry import Polygon from multirtc import dem -from multirtc.define_geogrid import get_point_epsg def check_poly_order(poly): @@ -40,6 +38,7 @@ class CapellaSICD: wavelength: float polarization: str lookside: str # 'right' or 'left' + hae: float shape: tuple[int, int] scp_index: tuple[int, int] range_pixel_spacing: float @@ -51,20 +50,6 @@ class CapellaSICD: orbit: isce3.core.Orbit beta0_coeff: np.ndarray - @staticmethod - def get_arp_pos_at_time(time, arp_pos_poly): - x = polyval(time, arp_pos_poly.X.get_array()) - y = polyval(time, arp_pos_poly.Y.get_array()) - z = polyval(time, arp_pos_poly.Z.get_array()) - return np.array([x, y, z]) - - @staticmethod - def get_arp_vel_at_time(time, arp_pos_poly): - x = polyval(time, polyder(arp_pos_poly.X.get_array())) - y = polyval(time, polyder(arp_pos_poly.Y.get_array())) - z = polyval(time, polyder(arp_pos_poly.Z.get_array())) - return np.array([x, y, z]) - @staticmethod def calculate_orbit(epoch, sensing_start, sensing_end, arp_pos_poly): svs = [] @@ -72,8 +57,8 @@ def calculate_orbit(epoch, sensing_start, sensing_end, arp_pos_poly): orbit_end = np.ceil(sensing_end) + 5 for offset_sec in np.arange(orbit_start, orbit_end + 1, 1): t = sensing_start + offset_sec - pos = CapellaSICD.get_arp_pos_at_time(t, arp_pos_poly) - vel = CapellaSICD.get_arp_vel_at_time(t, arp_pos_poly) + pos = arp_pos_poly(t) + vel = arp_pos_poly.derivative_eval(t) t_isce = to_isce_datetime(epoch + timedelta(seconds=t)) svs.append(isce3.core.StateVector(t_isce, pos, vel)) return isce3.core.Orbit(svs, to_isce_datetime(epoch)) @@ -104,26 +89,23 @@ def from_sarpy_sicd(cls, sicd, file_path): footprint = Polygon([(ic.Lon, ic.Lat) for ic in sicd.GeoData.ImageCorners]) shape = np.array([sicd.ImageData.NumRows, sicd.ImageData.NumCols]) scp_index = np.array([sicd.ImageData.SCPPixel.Row, sicd.ImageData.SCPPixel.Col]) - start_index = scp_index - shape - r_ss = sicd.Grid.Row.SS - c_ss = sicd.Grid.Col.SS - range_pixel_spacing = r_ss + row_shift = sicd.ImageData.SCPPixel.Row - sicd.ImageData.FirstRow + col_shift = sicd.ImageData.SCPPixel.Col - sicd.ImageData.FirstCol + row_mult = sicd.Grid.Row.SS + col_mult = sicd.Grid.Col.SS + range_pixel_spacing = row_mult assert sicd.Grid.Type == 'RGZERO', 'Only range zero doppler grids supported for Capella data' collect_start = sicd.Timeline.CollectStart # seconds after collect start - tcoa_coefs = sicd.Grid.TimeCOAPoly.Coefs - first_col_time = polyval2d(start_index[0] * r_ss, start_index[1] * c_ss, tcoa_coefs) - last_col_time = polyval2d(start_index[0] * r_ss, (start_index[1] + shape[1]) * c_ss, tcoa_coefs) - sensing_start = min(first_col_time, last_col_time) + first_col_time = sicd.Grid.TimeCOAPoly((0 - row_shift) * row_mult, (0 - col_shift) * col_mult) + last_col_time = sicd.Grid.TimeCOAPoly((shape[0] - row_shift) * row_mult, (shape[1] - col_shift) * col_mult) + sensing_start = min(first_col_time, last_col_time) # + 2 # fudged sensing_end = max(first_col_time, last_col_time) - prf = np.mean([ipp.IPPPoly.Coefs[1] for ipp in sicd.Timeline.IPP]) - starting_grid_pos = ( - sicd.Grid.Row.UVectECF.get_array() * r_ss * start_index[0] - + sicd.Grid.Col.UVectECF.get_array() * c_ss * start_index[1] - + sicd.GeoData.SCP.ECF.get_array() + prf = np.mean([ipp.IPPPoly.derivative_eval((ipp.TStart + ipp.TEnd) / 2) for ipp in sicd.Timeline.IPP]) + starting_row_pos = ( + sicd.GeoData.SCP.ECF.get_array() + sicd.Grid.Row.UVectECF.get_array() * (0 - row_shift) * row_mult ) - starting_arp_pos = cls.get_arp_pos_at_time(sensing_start, sicd.Position.ARPPoly) - starting_range = np.linalg.norm(starting_grid_pos - starting_arp_pos) + starting_range = np.linalg.norm(sicd.SCPCOA.ARPPos.get_array() - starting_row_pos) orbit = CapellaSICD.calculate_orbit(collect_start, sensing_start, sensing_end, sicd.Position.ARPPoly) beta0_coeff = sicd.Radiometric.BetaZeroSFPoly.Coefs capella_sicd = cls( @@ -134,6 +116,7 @@ def from_sarpy_sicd(cls, sicd, file_path): wavelength=wavelength, polarization=polarization, lookside=lookside, + hae=float(sicd.GeoData.SCP.LLH.HAE), range_pixel_spacing=range_pixel_spacing, scp_index=scp_index, collect_start=collect_start, From c8c7aaa3ce9446897256f72284dbc619275b892b Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Thu, 22 May 2025 08:22:30 -0500 Subject: [PATCH 04/12] fix row for tcoa --- src/multirtc/prep_capella.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/multirtc/prep_capella.py b/src/multirtc/prep_capella.py index dae5626..8e75fa4 100644 --- a/src/multirtc/prep_capella.py +++ b/src/multirtc/prep_capella.py @@ -98,7 +98,7 @@ def from_sarpy_sicd(cls, sicd, file_path): collect_start = sicd.Timeline.CollectStart # seconds after collect start first_col_time = sicd.Grid.TimeCOAPoly((0 - row_shift) * row_mult, (0 - col_shift) * col_mult) - last_col_time = sicd.Grid.TimeCOAPoly((shape[0] - row_shift) * row_mult, (shape[1] - col_shift) * col_mult) + last_col_time = sicd.Grid.TimeCOAPoly((0 - row_shift) * row_mult, (shape[1] - col_shift) * col_mult) sensing_start = min(first_col_time, last_col_time) # + 2 # fudged sensing_end = max(first_col_time, last_col_time) prf = np.mean([ipp.IPPPoly.derivative_eval((ipp.TStart + ipp.TEnd) / 2) for ipp in sicd.Timeline.IPP]) From aa950d9692bc703857caeeefac708b986d444cdf Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Thu, 22 May 2025 08:37:31 -0500 Subject: [PATCH 05/12] use INCA azimuth times --- src/multirtc/prep_capella.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/multirtc/prep_capella.py b/src/multirtc/prep_capella.py index 8e75fa4..bb3fb44 100644 --- a/src/multirtc/prep_capella.py +++ b/src/multirtc/prep_capella.py @@ -53,8 +53,8 @@ class CapellaSICD: @staticmethod def calculate_orbit(epoch, sensing_start, sensing_end, arp_pos_poly): svs = [] - orbit_start = np.floor(sensing_start) - 5 - orbit_end = np.ceil(sensing_end) + 5 + orbit_start = np.floor(sensing_start) - 10 + orbit_end = np.ceil(sensing_end) + 10 for offset_sec in np.arange(orbit_start, orbit_end + 1, 1): t = sensing_start + offset_sec pos = arp_pos_poly(t) @@ -91,19 +91,18 @@ def from_sarpy_sicd(cls, sicd, file_path): scp_index = np.array([sicd.ImageData.SCPPixel.Row, sicd.ImageData.SCPPixel.Col]) row_shift = sicd.ImageData.SCPPixel.Row - sicd.ImageData.FirstRow col_shift = sicd.ImageData.SCPPixel.Col - sicd.ImageData.FirstCol - row_mult = sicd.Grid.Row.SS - col_mult = sicd.Grid.Col.SS - range_pixel_spacing = row_mult + range_pixel_spacing = sicd.Grid.Row.SS assert sicd.Grid.Type == 'RGZERO', 'Only range zero doppler grids supported for Capella data' collect_start = sicd.Timeline.CollectStart # seconds after collect start - first_col_time = sicd.Grid.TimeCOAPoly((0 - row_shift) * row_mult, (0 - col_shift) * col_mult) - last_col_time = sicd.Grid.TimeCOAPoly((0 - row_shift) * row_mult, (shape[1] - col_shift) * col_mult) - sensing_start = min(first_col_time, last_col_time) # + 2 # fudged + first_col_time = sicd.RMA.INCA.TimeCAPoly(0 - col_shift) + last_col_time = sicd.RMA.INCA.TimeCAPoly(shape[1] - col_shift) + sensing_start = min(first_col_time, last_col_time) sensing_end = max(first_col_time, last_col_time) prf = np.mean([ipp.IPPPoly.derivative_eval((ipp.TStart + ipp.TEnd) / 2) for ipp in sicd.Timeline.IPP]) starting_row_pos = ( - sicd.GeoData.SCP.ECF.get_array() + sicd.Grid.Row.UVectECF.get_array() * (0 - row_shift) * row_mult + sicd.GeoData.SCP.ECF.get_array() + + sicd.Grid.Row.UVectECF.get_array() * (0 - row_shift) * range_pixel_spacing ) starting_range = np.linalg.norm(sicd.SCPCOA.ARPPos.get_array() - starting_row_pos) orbit = CapellaSICD.calculate_orbit(collect_start, sensing_start, sensing_end, sicd.Position.ARPPoly) From c1a5f9d27e2310baae51c74c2f44064bb520ddae Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Thu, 22 May 2025 08:44:30 -0500 Subject: [PATCH 06/12] alternative PRF --- src/multirtc/prep_capella.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/multirtc/prep_capella.py b/src/multirtc/prep_capella.py index bb3fb44..f767226 100644 --- a/src/multirtc/prep_capella.py +++ b/src/multirtc/prep_capella.py @@ -99,7 +99,8 @@ def from_sarpy_sicd(cls, sicd, file_path): last_col_time = sicd.RMA.INCA.TimeCAPoly(shape[1] - col_shift) sensing_start = min(first_col_time, last_col_time) sensing_end = max(first_col_time, last_col_time) - prf = np.mean([ipp.IPPPoly.derivative_eval((ipp.TStart + ipp.TEnd) / 2) for ipp in sicd.Timeline.IPP]) + # prf = np.mean([ipp.IPPPoly.derivative_eval((ipp.TStart + ipp.TEnd) / 2) for ipp in sicd.Timeline.IPP]) + prf = shape[1] / (sensing_end - sensing_start) starting_row_pos = ( sicd.GeoData.SCP.ECF.get_array() + sicd.Grid.Row.UVectECF.get_array() * (0 - row_shift) * range_pixel_spacing From 989f6807088a51bd99e16af0e9d78b8576fc1245 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Fri, 23 May 2025 14:17:31 -0500 Subject: [PATCH 07/12] working capella rtc --- src/multirtc/create_rtc.py | 17 +++++------ src/multirtc/multirtc.py | 11 ++++--- src/multirtc/prep_capella.py | 56 +++++++++++++++++++++++++++++++++--- 3 files changed, 66 insertions(+), 18 deletions(-) diff --git a/src/multirtc/create_rtc.py b/src/multirtc/create_rtc.py index a2508ed..fe88861 100755 --- a/src/multirtc/create_rtc.py +++ b/src/multirtc/create_rtc.py @@ -550,11 +550,11 @@ def run_single_job(product_id: str, burst: Sentinel1BurstSlc, geogrid, opts: Rtc logger.info(f'elapsed time: {t_end - t_start}') -def capella_rtc(capella_sicd, geogrid, opts): +def capella_rtc(sicd, geogrid, opts): # Common initializations t_start = time.time() output_dir = str(opts.output_dir) - product_id = capella_sicd.id + product_id = sicd.id os.makedirs(output_dir, exist_ok=True) raster_format = 'GTiff' @@ -567,14 +567,14 @@ def capella_rtc(capella_sicd, geogrid, opts): rtc_anf_gamma0_to_sigma0_file = ( f'{output_dir}/{product_id}_{LAYER_NAME_RTC_ANF_GAMMA0_TO_SIGMA0}.{raster_extension}' ) - radar_grid = capella_sicd.as_isce3_radargrid() - orbit = capella_sicd.orbit - wavelength = capella_sicd.wavelength + radar_grid = sicd.as_isce3_radargrid() + orbit = sicd.orbit + wavelength = sicd.wavelength lookside = radar_grid.lookside dem_raster = isce3.io.Raster(opts.dem_path) ellipsoid = isce3.core.Ellipsoid() - doppler = capella_sicd.get_doppler_centroid_grid() + doppler = sicd.get_doppler_centroid_grid() exponent = 2 x_snap = geogrid.spacing_x @@ -582,8 +582,9 @@ def capella_rtc(capella_sicd, geogrid, opts): geogrid.start_x = np.floor(float(geogrid.start_x) / x_snap) * x_snap geogrid.start_y = np.ceil(float(geogrid.start_y) / y_snap) * y_snap - # input_filename = save_as_beta0(umbra_sicd, Path(opts.output_dir)) - input_filename = 'beta0.tif' + input_filename = sicd.file_path.parent / (sicd.file_path.stem + '_beta0.tif') + if not input_filename.exists(): + sicd.create_complex_beta0(input_filename) input_filename = str(input_filename) # geocoding optional arguments diff --git a/src/multirtc/multirtc.py b/src/multirtc/multirtc.py index d478315..f1ce95a 100644 --- a/src/multirtc/multirtc.py +++ b/src/multirtc/multirtc.py @@ -4,7 +4,7 @@ import isce3 -from multirtc.create_rtc import run_single_job, umbra_rtc +from multirtc.create_rtc import capella_rtc, run_single_job, umbra_rtc from multirtc.define_geogrid import generate_geogrids from multirtc.prep_burst import prep_burst from multirtc.prep_capella import prep_capella @@ -13,7 +13,7 @@ def print_wkt(sicd): - radar_grid = sicd.as_isce3_radar_grid() + radar_grid = sicd.as_isce3_radargrid() dem = isce3.geometry.DEMInterpolator(sicd.hae) doppler = sicd.get_doppler_centroid_grid() wkt = isce3.geometry.get_geo_perimeter_wkt( @@ -80,10 +80,9 @@ def opera_rtc_capella_sicd(granule: str, resolution: int = 30, work_dir: Optiona raise FileNotFoundError(f'Capella SICD must be present in input dir {input_dir} for processing.') [d.mkdir(parents=True, exist_ok=True) for d in [input_dir, output_dir]] capella_sicd, dem_path = prep_capella(granule_path, work_dir=input_dir) - print_wkt(capella_sicd) - # opts = RtcOptions(dem_path=str(dem_path), output_dir=str(output_dir), resolution=resolution) - # geogrid = generate_geogrids(burst, opts.resolution) - # run_single_job(granule, burst, geogrid, opts) + opts = RtcOptions(dem_path=str(dem_path), output_dir=str(output_dir), resolution=resolution) + geogrid = generate_geogrids(capella_sicd, opts.resolution) + capella_rtc(capella_sicd, geogrid, opts) def main(): diff --git a/src/multirtc/prep_capella.py b/src/multirtc/prep_capella.py index f767226..f4b3510 100644 --- a/src/multirtc/prep_capella.py +++ b/src/multirtc/prep_capella.py @@ -5,12 +5,17 @@ import isce3 import numpy as np +from numpy.polynomial.polynomial import polyval2d +from osgeo import gdal from sarpy.io.complex.sicd import SICDReader from shapely.geometry import Polygon from multirtc import dem +gdal.UseExceptions() + + def check_poly_order(poly): assert len(poly.Coefs) == poly.order1 + 1, 'Polynomial order does not match number of coefficients' @@ -35,16 +40,20 @@ class CapellaSICD: id: str file_path: Path footprint: Polygon + center: Point wavelength: float polarization: str lookside: str # 'right' or 'left' hae: float shape: tuple[int, int] scp_index: tuple[int, int] + row_mult: float + col_mult: float range_pixel_spacing: float collect_start: datetime sensing_start: float sensing_end: float + az_reversed: bool prf: float starting_range: float orbit: isce3.core.Orbit @@ -63,7 +72,7 @@ def calculate_orbit(epoch, sensing_start, sensing_end, arp_pos_poly): svs.append(isce3.core.StateVector(t_isce, pos, vel)) return isce3.core.Orbit(svs, to_isce_datetime(epoch)) - def as_isce3_radar_grid(self): + def as_isce3_radargrid(self): radar_grid = isce3.product.RadarGridParameters( sensing_start=self.sensing_start, wavelength=self.wavelength, @@ -80,6 +89,42 @@ def as_isce3_radar_grid(self): def get_doppler_centroid_grid(self): return isce3.core.LUT2d() + def get_xrow_ycol(self) -> np.ndarray: + """Calculate xrow and ycol SICD.""" + irow = np.tile(np.arange(self.shape[0]), (self.shape[1], 1)).T + irow -= self.scp_index[0] + xrow = irow * self.row_mult + + icol = np.tile(np.arange(self.shape[1]), (self.shape[0], 1)) + icol -= self.scp_index[1] + ycol = icol * self.col_mult + return xrow, ycol + + def load_data(self): + """Load data from the UMBRA SICD file.""" + reader = SICDReader(str(self.file_path)) + data = reader[:, :] + return data + + def create_complex_beta0(self, outpath, isce_format=True): + xrow, ycol = self.get_xrow_ycol() + scale_factor = np.sqrt(polyval2d(xrow, ycol, self.beta0_coeff)) + data = self.load_data() + scaled_data = data * scale_factor + + if isce_format: + if self.az_reversed: + scaled_data = scaled_data[:, ::-1].T + else: + scaled_data = scaled_data.T + + driver = gdal.GetDriverByName('GTiff') + ds = driver.Create(str(outpath), scaled_data.shape[1], scaled_data.shape[0], 1, gdal.GDT_CFloat32) + band = ds.GetRasterBand(1) + band.WriteArray(scaled_data) + band.FlushCache() + ds = None + @classmethod def from_sarpy_sicd(cls, sicd, file_path): center_frequency = sicd.RadarCollection.TxFrequency.Min + sicd.RadarCollection.TxFrequency.Max / 2 @@ -99,7 +144,7 @@ def from_sarpy_sicd(cls, sicd, file_path): last_col_time = sicd.RMA.INCA.TimeCAPoly(shape[1] - col_shift) sensing_start = min(first_col_time, last_col_time) sensing_end = max(first_col_time, last_col_time) - # prf = np.mean([ipp.IPPPoly.derivative_eval((ipp.TStart + ipp.TEnd) / 2) for ipp in sicd.Timeline.IPP]) + az_reversed = last_col_time < first_col_time prf = shape[1] / (sensing_end - sensing_start) starting_row_pos = ( sicd.GeoData.SCP.ECF.get_array() @@ -107,25 +152,28 @@ def from_sarpy_sicd(cls, sicd, file_path): ) starting_range = np.linalg.norm(sicd.SCPCOA.ARPPos.get_array() - starting_row_pos) orbit = CapellaSICD.calculate_orbit(collect_start, sensing_start, sensing_end, sicd.Position.ARPPoly) - beta0_coeff = sicd.Radiometric.BetaZeroSFPoly.Coefs capella_sicd = cls( id=Path(file_path).with_suffix('').name, file_path=file_path, footprint=footprint, shape=shape, + center=Point(sicd.GeoData.SCP.LLH.Lon, sicd.GeoData.SCP.LLH.Lat), wavelength=wavelength, polarization=polarization, lookside=lookside, hae=float(sicd.GeoData.SCP.LLH.HAE), range_pixel_spacing=range_pixel_spacing, scp_index=scp_index, + col_mult=sicd.Grid.Col.SS, + row_mult=sicd.Grid.Row.SS, collect_start=collect_start, sensing_start=sensing_start, sensing_end=sensing_end, + az_reversed=az_reversed, prf=prf, starting_range=starting_range, orbit=orbit, - beta0_coeff=beta0_coeff, + beta0_coeff=sicd.Radiometric.BetaZeroSFPoly.Coefs, ) return capella_sicd From 696a3d88c54e830ae1112b652d2a3537c8b1b38a Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 26 May 2025 19:33:29 +0000 Subject: [PATCH 08/12] Bump ASFHyP3/actions from 0.18.1 to 0.20.0 Bumps [ASFHyP3/actions](https://github.com/asfhyp3/actions) from 0.18.1 to 0.20.0. - [Release notes](https://github.com/asfhyp3/actions/releases) - [Changelog](https://github.com/ASFHyP3/actions/blob/develop/CHANGELOG.md) - [Commits](https://github.com/asfhyp3/actions/compare/v0.18.1...v0.20.0) --- updated-dependencies: - dependency-name: ASFHyP3/actions dependency-version: 0.20.0 dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] --- .github/workflows/changelog.yml | 2 +- .github/workflows/labeled-pr.yml | 2 +- .github/workflows/release-checklist-comment.yml | 2 +- .github/workflows/release.yml | 2 +- .github/workflows/static-analysis.yml | 4 ++-- .github/workflows/tag-version.yml | 2 +- .github/workflows/test-and-build.yml | 2 +- 7 files changed, 8 insertions(+), 8 deletions(-) diff --git a/.github/workflows/changelog.yml b/.github/workflows/changelog.yml index 3e4a03f..5360e08 100644 --- a/.github/workflows/changelog.yml +++ b/.github/workflows/changelog.yml @@ -14,4 +14,4 @@ on: jobs: call-changelog-check-workflow: # Docs: https://github.com/ASFHyP3/actions - uses: ASFHyP3/actions/.github/workflows/reusable-changelog-check.yml@v0.18.1 + uses: ASFHyP3/actions/.github/workflows/reusable-changelog-check.yml@v0.20.0 diff --git a/.github/workflows/labeled-pr.yml b/.github/workflows/labeled-pr.yml index 4f8dcd9..d729153 100644 --- a/.github/workflows/labeled-pr.yml +++ b/.github/workflows/labeled-pr.yml @@ -13,4 +13,4 @@ on: jobs: call-labeled-pr-check-workflow: # Docs: https://github.com/ASFHyP3/actions - uses: ASFHyP3/actions/.github/workflows/reusable-labeled-pr-check.yml@v0.18.1 + uses: ASFHyP3/actions/.github/workflows/reusable-labeled-pr-check.yml@v0.20.0 diff --git a/.github/workflows/release-checklist-comment.yml b/.github/workflows/release-checklist-comment.yml index 5170ee3..a2bb05d 100644 --- a/.github/workflows/release-checklist-comment.yml +++ b/.github/workflows/release-checklist-comment.yml @@ -10,7 +10,7 @@ on: jobs: call-release-workflow: # Docs: https://github.com/ASFHyP3/actions - uses: ASFHyP3/actions/.github/workflows/reusable-release-checklist-comment.yml@v0.18.1 + uses: ASFHyP3/actions/.github/workflows/reusable-release-checklist-comment.yml@v0.20.0 permissions: pull-requests: write secrets: diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index d43f1f3..f49fd13 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -8,7 +8,7 @@ on: jobs: call-release-workflow: # Docs: https://github.com/ASFHyP3/actions - uses: ASFHyP3/actions/.github/workflows/reusable-release.yml@v0.18.1 + uses: ASFHyP3/actions/.github/workflows/reusable-release.yml@v0.20.0 with: release_prefix: MultiRTC release_branch: main diff --git a/.github/workflows/static-analysis.yml b/.github/workflows/static-analysis.yml index 15d26d0..e72b3ae 100644 --- a/.github/workflows/static-analysis.yml +++ b/.github/workflows/static-analysis.yml @@ -5,8 +5,8 @@ on: push jobs: call-secrets-analysis-workflow: # Docs: https://github.com/ASFHyP3/actions - uses: ASFHyP3/actions/.github/workflows/reusable-secrets-analysis.yml@v0.18.1 + uses: ASFHyP3/actions/.github/workflows/reusable-secrets-analysis.yml@v0.20.0 call-ruff-workflow: # Docs: https://github.com/ASFHyP3/actions - uses: ASFHyP3/actions/.github/workflows/reusable-ruff.yml@v0.18.1 + uses: ASFHyP3/actions/.github/workflows/reusable-ruff.yml@v0.20.0 diff --git a/.github/workflows/tag-version.yml b/.github/workflows/tag-version.yml index 4679b3c..1bfc245 100644 --- a/.github/workflows/tag-version.yml +++ b/.github/workflows/tag-version.yml @@ -8,7 +8,7 @@ on: jobs: call-bump-version-workflow: # Docs: https://github.com/ASFHyP3/actions - uses: ASFHyP3/actions/.github/workflows/reusable-bump-version.yml@v0.18.1 + uses: ASFHyP3/actions/.github/workflows/reusable-bump-version.yml@v0.20.0 with: user: forrest-bot email: ffwilliams2@alaska.edu diff --git a/.github/workflows/test-and-build.yml b/.github/workflows/test-and-build.yml index 1db4d2c..193470f 100644 --- a/.github/workflows/test-and-build.yml +++ b/.github/workflows/test-and-build.yml @@ -13,7 +13,7 @@ on: jobs: call-pytest-workflow: # Docs: https://github.com/ASFHyP3/actions - uses: ASFHyP3/actions/.github/workflows/reusable-pytest.yml@v0.18.1 + uses: ASFHyP3/actions/.github/workflows/reusable-pytest.yml@v0.20.0 with: local_package_name: mulitrtc python_versions: >- From 2083e624ed94bb25af40fc98bd3075e22d4c33e6 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Tue, 27 May 2025 09:00:52 -0500 Subject: [PATCH 09/12] switch capella class to new framework --- src/multirtc/base.py | 43 ++++++++ src/multirtc/create_rtc.py | 6 +- src/multirtc/define_geogrid.py | 15 ++- src/multirtc/prep_capella.py | 174 +-------------------------------- src/multirtc/sicd.py | 170 ++++++++++++++++++++++++++++++++ 5 files changed, 225 insertions(+), 183 deletions(-) create mode 100644 src/multirtc/base.py create mode 100644 src/multirtc/sicd.py diff --git a/src/multirtc/base.py b/src/multirtc/base.py new file mode 100644 index 0000000..476eba6 --- /dev/null +++ b/src/multirtc/base.py @@ -0,0 +1,43 @@ +from abc import ABC +from datetime import datetime +from pathlib import Path + +from shapely.geometry import Point, Polygon + + +class SlcTemplate(ABC): + required_attributes = { + 'id': str, + 'filepath': Path, + 'footprint': Polygon, + 'center': Point, + 'lookside': str, # 'right' or 'left' + 'wavelength': float, + 'polarization': str, + 'shape': tuple, + 'range_pixel_spacing': float, + 'reference_time': datetime, + 'sensing_start': float, + 'prf': float, + 'orbit': object, # Replace with actual orbit type + 'radar_grid': object, # Replace with actual radar grid type + 'doppler_centroid_grid': object, # Replace with actual doppler centroid grid type + } + + # I prefer this setup to enforce properties over forcing subclasses to have a bunch of @property statements + def __init_subclass__(cls): + super().__init_subclass__() + original_init = cls.__init__ + + def wrapped_init(self, *args, **kwargs): + original_init(self, *args, **kwargs) + for attr, expected_type in cls.required_attributes.items(): + if not hasattr(self, attr): + raise NotImplementedError(f'{cls.__name__} must define self.{attr}') + if not isinstance(getattr(self, attr), expected_type): + raise TypeError( + f'{cls.__name__}.{attr} must be of type {expected_type.__name__},' + f'got {type(getattr(self, attr)).__name__}' + ) + + cls.__init__ = wrapped_init diff --git a/src/multirtc/create_rtc.py b/src/multirtc/create_rtc.py index fe88861..206a5c9 100755 --- a/src/multirtc/create_rtc.py +++ b/src/multirtc/create_rtc.py @@ -567,14 +567,14 @@ def capella_rtc(sicd, geogrid, opts): rtc_anf_gamma0_to_sigma0_file = ( f'{output_dir}/{product_id}_{LAYER_NAME_RTC_ANF_GAMMA0_TO_SIGMA0}.{raster_extension}' ) - radar_grid = sicd.as_isce3_radargrid() + radar_grid = sicd.radar_grid orbit = sicd.orbit wavelength = sicd.wavelength lookside = radar_grid.lookside dem_raster = isce3.io.Raster(opts.dem_path) ellipsoid = isce3.core.Ellipsoid() - doppler = sicd.get_doppler_centroid_grid() + doppler = sicd.doppler_centroid_grid exponent = 2 x_snap = geogrid.spacing_x @@ -582,7 +582,7 @@ def capella_rtc(sicd, geogrid, opts): geogrid.start_x = np.floor(float(geogrid.start_x) / x_snap) * x_snap geogrid.start_y = np.ceil(float(geogrid.start_y) / y_snap) * y_snap - input_filename = sicd.file_path.parent / (sicd.file_path.stem + '_beta0.tif') + input_filename = sicd.filepath.parent / (sicd.filepath.stem + '_beta0.tif') if not input_filename.exists(): sicd.create_complex_beta0(input_filename) input_filename = str(input_filename) diff --git a/src/multirtc/define_geogrid.py b/src/multirtc/define_geogrid.py index d3d0aca..87c85c5 100644 --- a/src/multirtc/define_geogrid.py +++ b/src/multirtc/define_geogrid.py @@ -128,13 +128,12 @@ def generate_geogrids(slc_obj, resolution: int, epsg: int = None, rda: bool = Tr if epsg is None: epsg = get_point_epsg(slc_obj.center.y, slc_obj.center.x) - if rda: - radar_grid = slc_obj.as_isce3_radargrid() - geogrid = isce3.product.bbox_to_geogrid( - radar_grid, slc_obj.orbit, isce3.core.LUT2d(), x_spacing, y_spacing, epsg - ) - else: - geogrid = slc_obj.get_geogrid(x_spacing, y_spacing) - + # if rda: + # radar_grid = slc_obj.as_isce3_radargrid() + # else: + # geogrid = slc_obj.get_geogrid(x_spacing, y_spacing) + geogrid = isce3.product.bbox_to_geogrid( + slc_obj.radar_grid, slc_obj.orbit, slc_obj.doppler_centroid_grid, x_spacing, y_spacing, epsg + ) geogrid_snapped = snap_geogrid(geogrid, geogrid.spacing_x, geogrid.spacing_y) return geogrid_snapped diff --git a/src/multirtc/prep_capella.py b/src/multirtc/prep_capella.py index f4b3510..ad8bcfd 100644 --- a/src/multirtc/prep_capella.py +++ b/src/multirtc/prep_capella.py @@ -1,183 +1,15 @@ -from dataclasses import dataclass -from datetime import datetime, timedelta from pathlib import Path from typing import Optional -import isce3 -import numpy as np -from numpy.polynomial.polynomial import polyval2d from osgeo import gdal -from sarpy.io.complex.sicd import SICDReader -from shapely.geometry import Polygon from multirtc import dem +from multirtc.sicd import SicdRzdSlc gdal.UseExceptions() -def check_poly_order(poly): - assert len(poly.Coefs) == poly.order1 + 1, 'Polynomial order does not match number of coefficients' - - -def to_isce_datetime(dt): - if isinstance(dt, datetime): - return isce3.core.DateTime(dt) - elif isinstance(dt, np.datetime64): - return isce3.core.DateTime(dt.item()) - else: - raise ValueError(f'Unsupported datetime type: {type(dt)}. Expected datetime or np.datetime64.') - - -@dataclass -class Point: - x: float - y: float - - -@dataclass -class CapellaSICD: - id: str - file_path: Path - footprint: Polygon - center: Point - wavelength: float - polarization: str - lookside: str # 'right' or 'left' - hae: float - shape: tuple[int, int] - scp_index: tuple[int, int] - row_mult: float - col_mult: float - range_pixel_spacing: float - collect_start: datetime - sensing_start: float - sensing_end: float - az_reversed: bool - prf: float - starting_range: float - orbit: isce3.core.Orbit - beta0_coeff: np.ndarray - - @staticmethod - def calculate_orbit(epoch, sensing_start, sensing_end, arp_pos_poly): - svs = [] - orbit_start = np.floor(sensing_start) - 10 - orbit_end = np.ceil(sensing_end) + 10 - for offset_sec in np.arange(orbit_start, orbit_end + 1, 1): - t = sensing_start + offset_sec - pos = arp_pos_poly(t) - vel = arp_pos_poly.derivative_eval(t) - t_isce = to_isce_datetime(epoch + timedelta(seconds=t)) - svs.append(isce3.core.StateVector(t_isce, pos, vel)) - return isce3.core.Orbit(svs, to_isce_datetime(epoch)) - - def as_isce3_radargrid(self): - radar_grid = isce3.product.RadarGridParameters( - sensing_start=self.sensing_start, - wavelength=self.wavelength, - prf=self.prf, - starting_range=self.starting_range, - range_pixel_spacing=self.range_pixel_spacing, - lookside=isce3.core.LookSide.Right if self.lookside == 'right' else isce3.core.LookSide.Left, - length=self.shape[1], # flipped for "shadows down" convention - width=self.shape[0], # flipped for "shadows down" convention - ref_epoch=to_isce_datetime(self.collect_start), - ) - return radar_grid - - def get_doppler_centroid_grid(self): - return isce3.core.LUT2d() - - def get_xrow_ycol(self) -> np.ndarray: - """Calculate xrow and ycol SICD.""" - irow = np.tile(np.arange(self.shape[0]), (self.shape[1], 1)).T - irow -= self.scp_index[0] - xrow = irow * self.row_mult - - icol = np.tile(np.arange(self.shape[1]), (self.shape[0], 1)) - icol -= self.scp_index[1] - ycol = icol * self.col_mult - return xrow, ycol - - def load_data(self): - """Load data from the UMBRA SICD file.""" - reader = SICDReader(str(self.file_path)) - data = reader[:, :] - return data - - def create_complex_beta0(self, outpath, isce_format=True): - xrow, ycol = self.get_xrow_ycol() - scale_factor = np.sqrt(polyval2d(xrow, ycol, self.beta0_coeff)) - data = self.load_data() - scaled_data = data * scale_factor - - if isce_format: - if self.az_reversed: - scaled_data = scaled_data[:, ::-1].T - else: - scaled_data = scaled_data.T - - driver = gdal.GetDriverByName('GTiff') - ds = driver.Create(str(outpath), scaled_data.shape[1], scaled_data.shape[0], 1, gdal.GDT_CFloat32) - band = ds.GetRasterBand(1) - band.WriteArray(scaled_data) - band.FlushCache() - ds = None - - @classmethod - def from_sarpy_sicd(cls, sicd, file_path): - center_frequency = sicd.RadarCollection.TxFrequency.Min + sicd.RadarCollection.TxFrequency.Max / 2 - wavelength = isce3.core.speed_of_light / center_frequency - polarization = sicd.RadarCollection.RcvChannels[0].TxRcvPolarization.replace(':', '') - lookside = 'right' if sicd.SCPCOA.SideOfTrack == 'R' else 'left' - footprint = Polygon([(ic.Lon, ic.Lat) for ic in sicd.GeoData.ImageCorners]) - shape = np.array([sicd.ImageData.NumRows, sicd.ImageData.NumCols]) - scp_index = np.array([sicd.ImageData.SCPPixel.Row, sicd.ImageData.SCPPixel.Col]) - row_shift = sicd.ImageData.SCPPixel.Row - sicd.ImageData.FirstRow - col_shift = sicd.ImageData.SCPPixel.Col - sicd.ImageData.FirstCol - range_pixel_spacing = sicd.Grid.Row.SS - assert sicd.Grid.Type == 'RGZERO', 'Only range zero doppler grids supported for Capella data' - collect_start = sicd.Timeline.CollectStart - # seconds after collect start - first_col_time = sicd.RMA.INCA.TimeCAPoly(0 - col_shift) - last_col_time = sicd.RMA.INCA.TimeCAPoly(shape[1] - col_shift) - sensing_start = min(first_col_time, last_col_time) - sensing_end = max(first_col_time, last_col_time) - az_reversed = last_col_time < first_col_time - prf = shape[1] / (sensing_end - sensing_start) - starting_row_pos = ( - sicd.GeoData.SCP.ECF.get_array() - + sicd.Grid.Row.UVectECF.get_array() * (0 - row_shift) * range_pixel_spacing - ) - starting_range = np.linalg.norm(sicd.SCPCOA.ARPPos.get_array() - starting_row_pos) - orbit = CapellaSICD.calculate_orbit(collect_start, sensing_start, sensing_end, sicd.Position.ARPPoly) - capella_sicd = cls( - id=Path(file_path).with_suffix('').name, - file_path=file_path, - footprint=footprint, - shape=shape, - center=Point(sicd.GeoData.SCP.LLH.Lon, sicd.GeoData.SCP.LLH.Lat), - wavelength=wavelength, - polarization=polarization, - lookside=lookside, - hae=float(sicd.GeoData.SCP.LLH.HAE), - range_pixel_spacing=range_pixel_spacing, - scp_index=scp_index, - col_mult=sicd.Grid.Col.SS, - row_mult=sicd.Grid.Row.SS, - collect_start=collect_start, - sensing_start=sensing_start, - sensing_end=sensing_end, - az_reversed=az_reversed, - prf=prf, - starting_range=starting_range, - orbit=orbit, - beta0_coeff=sicd.Radiometric.BetaZeroSFPoly.Coefs, - ) - return capella_sicd - - def prep_capella(granule_path: Path, work_dir: Optional[Path] = None) -> Path: """Prepare data for burst-based processing. @@ -187,9 +19,7 @@ def prep_capella(granule_path: Path, work_dir: Optional[Path] = None) -> Path: """ if work_dir is None: work_dir = Path.cwd() - reader = SICDReader(str(granule_path)) - sicd = reader.get_sicds_as_tuple()[0] - capella_sicd = CapellaSICD.from_sarpy_sicd(sicd, granule_path) + capella_sicd = SicdRzdSlc(granule_path) dem_path = work_dir / 'dem.tif' dem.download_opera_dem_for_footprint(dem_path, capella_sicd.footprint) return capella_sicd, dem_path diff --git a/src/multirtc/sicd.py b/src/multirtc/sicd.py new file mode 100644 index 0000000..a2ffe91 --- /dev/null +++ b/src/multirtc/sicd.py @@ -0,0 +1,170 @@ +from datetime import datetime, timedelta +from pathlib import Path + +import isce3 +import numpy as np +from numpy.polynomial.polynomial import polyval2d +from osgeo import gdal +from sarpy.io.complex.sicd import SICDReader +from shapely.geometry import Point, Polygon + +from multirtc.base import SlcTemplate + + +def to_isce_datetime(dt): + if isinstance(dt, datetime): + return isce3.core.DateTime(dt) + elif isinstance(dt, np.datetime64): + return isce3.core.DateTime(dt.item()) + else: + raise ValueError(f'Unsupported datetime type: {type(dt)}. Expected datetime or np.datetime64.') + + +def check_poly_order(poly): + assert len(poly.Coefs) == poly.order1 + 1, 'Polynomial order does not match number of coefficients' + + +class Sicd: + def __init__(self, sicd_path): + reader = SICDReader(str(sicd_path.expanduser().resolve())) + sicd = reader.get_sicds_as_tuple()[0] + self.source = sicd + self.id = Path(sicd_path).with_suffix('').name + self.filepath = Path(sicd_path) + self.footprint = Polygon([(ic.Lon, ic.Lat) for ic in sicd.GeoData.ImageCorners]) + self.center = Point(sicd.GeoData.SCP.LLH.Lon, sicd.GeoData.SCP.LLH.Lat) + self.lookside = 'right' if sicd.SCPCOA.SideOfTrack == 'R' else 'left' + + center_frequency = sicd.RadarCollection.TxFrequency.Min + sicd.RadarCollection.TxFrequency.Max / 2 + self.wavelength = isce3.core.speed_of_light / center_frequency + self.polarization = sicd.RadarCollection.RcvChannels[0].TxRcvPolarization.replace(':', '') + self.shape = (sicd.ImageData.NumRows, sicd.ImageData.NumCols) + self.spacing = (sicd.Grid.Row.SS, sicd.Grid.Col.SS) + self.scp_index = (sicd.ImageData.SCPPixel.Row, sicd.ImageData.SCPPixel.Col) + self.range_pixel_spacing = sicd.Grid.Row.SS + self.reference_time = sicd.Timeline.CollectStart.item() + self.shift = ( + sicd.ImageData.SCPPixel.Row - sicd.ImageData.FirstRow, + sicd.ImageData.SCPPixel.Col - sicd.ImageData.FirstCol, + ) + self.arp_pos_poly = sicd.Position.ARPPoly + starting_row_pos = ( + sicd.GeoData.SCP.ECF.get_array() + + sicd.Grid.Row.UVectECF.get_array() * (0 - self.shift[0]) * self.spacing[0] + ) + self.starting_range = np.linalg.norm(sicd.SCPCOA.ARPPos.get_array() - starting_row_pos) + last_line_time = sicd.Grid.TimeCOAPoly(0, self.shape[1] - self.shift[1]) + first_line_time = sicd.Grid.TimeCOAPoly(0, -self.shift[1]) + self.az_reversed = last_line_time >= first_line_time + self.beta0 = sicd.Radiometric.BetaZeroSFPoly + self.sigma0 = sicd.Radiometric.SigmaZeroSFPoly + + def get_xrow_ycol(self) -> np.ndarray: + """Calculate xrow and ycol SICD.""" + irow = np.tile(np.arange(self.shape[0]), (self.shape[1], 1)).T + irow -= self.scp_index[0] + xrow = irow * self.row_mult + + icol = np.tile(np.arange(self.shape[1]), (self.shape[0], 1)) + icol -= self.scp_index[1] + ycol = icol * self.col_mult + return xrow, ycol + + def load_data(self): + return self.source[:, :] + + def load_scaled_data(self, scale, power=False): + if scale == 'beta0': + coeff = self.beta0.Coefs + elif scale == 'sigma0': + coeff = self.sigma0.Coefs + else: + raise ValueError(f'Scale must be either "beta0" or "sigma0", got {scale}') + + xrow, ycol = self.get_xrow_ycol() + scale_factor = polyval2d(xrow, ycol, coeff) + data = self.load_data() + if power: + scaled_data = (data.real**2 + data.imag**2) * scale_factor + else: + scaled_data = data * np.sqrt(scale_factor) + return scaled_data + + def write_complex_beta0(self, outpath, isce_format=True): + scaled_data = self.load_scaled_data('beta0', power=False) + if isce_format: + if self.az_reversed: + scaled_data = scaled_data[:, ::-1].T + else: + scaled_data = scaled_data.T + + driver = gdal.GetDriverByName('GTiff') + ds = driver.Create(str(outpath), scaled_data.shape[1], scaled_data.shape[0], 1, gdal.GDT_CFloat32) + band = ds.GetRasterBand(1) + band.WriteArray(scaled_data) + band.FlushCache() + ds = None + + def create_complex_beta0(self, outpath, isce_format=True): + xrow, ycol = self.get_xrow_ycol() + scale_factor = np.sqrt(polyval2d(xrow, ycol, self.beta0_coeff)) + data = self.load_data() + scaled_data = data * scale_factor + + if isce_format: + if self.az_reversed: + scaled_data = scaled_data[:, ::-1].T + else: + scaled_data = scaled_data.T + + driver = gdal.GetDriverByName('GTiff') + ds = driver.Create(str(outpath), scaled_data.shape[1], scaled_data.shape[0], 1, gdal.GDT_CFloat32) + band = ds.GetRasterBand(1) + band.WriteArray(scaled_data) + band.FlushCache() + ds = None + + +class SicdRzdSlc(SlcTemplate, Sicd): + def __init__(self, sicd_path): + super().__init__(sicd_path) + assert self.source.Grid.Type == 'RGZERO', 'Only range zero doppler grids supported for Capella data' + first_col_time = self.source.RMA.INCA.TimeCAPoly(0 - self.shift[1]) + last_col_time = self.source.RMA.INCA.TimeCAPoly(self.shape[1] - self.shift[1]) + self.sensing_start = min(first_col_time, last_col_time) + self.sensing_end = max(first_col_time, last_col_time) + self.prf = self.shape[1] / (self.sensing_end - self.sensing_start) + self.orbit = self.get_orbit() + self.radar_grid = self.get_radar_grid() + self.doppler_centroid_grid = isce3.core.LUT2d() + + def get_orbit(self): + svs = [] + orbit_start = np.floor(self.sensing_start) - 10 + orbit_end = np.ceil(self.sensing_end) + 10 + for offset_sec in np.arange(orbit_start, orbit_end + 1, 1): + t = self.sensing_start + offset_sec + pos = self.arp_pos_poly(t) + vel = self.arp_pos_poly.derivative_eval(t) + t_isce = to_isce_datetime(self.reference_time + timedelta(seconds=t)) + svs.append(isce3.core.StateVector(t_isce, pos, vel)) + return isce3.core.Orbit(svs, to_isce_datetime(self.reference_time)) + + def get_radar_grid(self): + radar_grid = isce3.product.RadarGridParameters( + sensing_start=self.sensing_start, + wavelength=self.wavelength, + prf=self.prf, + starting_range=self.starting_range, + range_pixel_spacing=self.range_pixel_spacing, + lookside=isce3.core.LookSide.Right if self.lookside == 'right' else isce3.core.LookSide.Left, + length=self.shape[1], # flipped for "shadows down" convention + width=self.shape[0], # flipped for "shadows down" convention + ref_epoch=to_isce_datetime(self.reference_time), + ) + return radar_grid + + +if __name__ == '__main__': + slc = SicdRzdSlc(Path('~/Data/capella/input/CAPELLA_C15_SS_SICD_VV_20250225215110_20250225215124.ntf')) + breakpoint() From 3cb0f88ded782f539e3572d1aa06b0d947417939 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Wed, 28 May 2025 08:42:22 -0500 Subject: [PATCH 10/12] unified s1 and capella --- src/multirtc/base.py | 15 +++++++ src/multirtc/create_rtc.py | 56 ++++++++++++++++++++++++-- src/multirtc/multirtc.py | 19 +++++++-- src/multirtc/prep_burst.py | 62 ++++++++++++++++++++++++----- src/multirtc/s1burst_corrections.py | 7 ++-- src/multirtc/sicd.py | 22 ++-------- 6 files changed, 142 insertions(+), 39 deletions(-) diff --git a/src/multirtc/base.py b/src/multirtc/base.py index 476eba6..e1257e4 100644 --- a/src/multirtc/base.py +++ b/src/multirtc/base.py @@ -2,9 +2,24 @@ from datetime import datetime from pathlib import Path +import isce3 +import numpy as np from shapely.geometry import Point, Polygon +def to_isce_datetime(dt): + if isinstance(dt, datetime): + return isce3.core.DateTime(dt) + elif isinstance(dt, np.datetime64): + return isce3.core.DateTime(dt.item()) + else: + raise ValueError(f'Unsupported datetime type: {type(dt)}. Expected datetime or np.datetime64.') + + +def from_isce_datetime(dt): + return datetime.fromisoformat(dt.isoformat()) + + class SlcTemplate(ABC): required_attributes = { 'id': str, diff --git a/src/multirtc/create_rtc.py b/src/multirtc/create_rtc.py index 206a5c9..25f770d 100755 --- a/src/multirtc/create_rtc.py +++ b/src/multirtc/create_rtc.py @@ -11,8 +11,10 @@ from scipy import ndimage from tqdm import tqdm +from multirtc.prep_burst import S1BurstSlc from multirtc.rtc_options import RtcOptions from multirtc.s1burst_corrections import apply_slc_corrections, compute_correction_lut +from multirtc.sicd import SicdSlc logger = logging.getLogger('rtc_s1') @@ -582,15 +584,48 @@ def capella_rtc(sicd, geogrid, opts): geogrid.start_x = np.floor(float(geogrid.start_x) / x_snap) * x_snap geogrid.start_y = np.ceil(float(geogrid.start_y) / y_snap) * y_snap - input_filename = sicd.filepath.parent / (sicd.filepath.stem + '_beta0.tif') - if not input_filename.exists(): - sicd.create_complex_beta0(input_filename) - input_filename = str(input_filename) + if isinstance(sicd, SicdSlc): + input_filename = sicd.filepath.parent / (sicd.filepath.stem + '_beta0.tif') + if not input_filename.exists(): + sicd.create_complex_beta0(input_filename) + input_filename = str(input_filename) + elif isinstance(sicd, S1BurstSlc) and opts.apply_thermal_noise or opts.apply_abs_rad: + # Convert to beta0 and apply thermal noise correction + input_filename = str(sicd.filepath.parent / (sicd.filepath.stem + '_beta0.tif')) + apply_slc_corrections( + sicd.source, + sicd.filepath, + input_filename, + flag_output_complex=True, + flag_thermal_correction=opts.apply_thermal_noise, + flag_apply_abs_rad_correction=opts.apply_abs_rad, + ) + else: + input_filename = str(sicd.filepath) # geocoding optional arguments geocode_kwargs = {} layover_shadow_mask_geocode_kwargs = {} + # get sub_swaths metadata + if isinstance(sicd, S1BurstSlc) and opts.apply_valid_samples_sub_swath_masking: + # Extract burst boundaries and create sub_swaths object to mask + # invalid radar samples + n_subswaths = 1 + sub_swaths = isce3.product.SubSwaths(radar_grid.length, radar_grid.width, n_subswaths) + last_range_sample = min([sicd.last_valid_sample, radar_grid.width]) + valid_samples_sub_swath = np.repeat( + [[sicd.first_valid_sample, last_range_sample + 1]], radar_grid.length, axis=0 + ) + for i in range(sicd.first_valid_line): + valid_samples_sub_swath[i, :] = 0 + for i in range(sicd.last_valid_line, radar_grid.length): + valid_samples_sub_swath[i, :] = 0 + + sub_swaths.set_valid_samples_array(1, valid_samples_sub_swath) + geocode_kwargs['sub_swaths'] = sub_swaths + layover_shadow_mask_geocode_kwargs['sub_swaths'] = sub_swaths + layover_shadow_mask_file = f'{output_dir}/{product_id}_{LAYER_NAME_LAYOVER_SHADOW_MASK}.{raster_extension}' logger.info(f' computing layover shadow mask for {product_id}') radar_grid_layover_shadow_mask = radar_grid @@ -621,6 +656,19 @@ def capella_rtc(sicd, geogrid, opts): rtc_anf_gamma0_to_sigma0_file, geogrid.width, geogrid.length, 1, gdal.GDT_Float32, raster_format ) geocode_kwargs['out_geo_rtc_gamma0_to_sigma0'] = out_geo_rtc_gamma0_to_sigma0_obj + if opts.apply_bistatic_delay or opts.apply_static_tropo: + rg_lut, az_lut = compute_correction_lut( + sicd.source, + dem_raster, + output_dir, + opts.correction_lut_range_spacing_in_meters, + opts.correction_lut_azimuth_spacing_in_meters, + opts.apply_bistatic_delay, + opts.apply_static_tropo, + ) + geocode_kwargs['az_time_correction'] = az_lut + if rg_lut is not None: + geocode_kwargs['slant_range_correction'] = rg_lut rdr_raster = isce3.io.Raster(input_filename) # Generate output geocoded burst raster diff --git a/src/multirtc/multirtc.py b/src/multirtc/multirtc.py index f1ce95a..3fe8f29 100644 --- a/src/multirtc/multirtc.py +++ b/src/multirtc/multirtc.py @@ -37,9 +37,16 @@ def opera_rtc_s1_burst(granule: str, resolution: int = 30, work_dir: Optional[Pa [d.mkdir(parents=True, exist_ok=True) for d in [input_dir, output_dir]] burst, dem_path = prep_burst(granule, work_dir=input_dir) - opts = RtcOptions(dem_path=str(dem_path), output_dir=str(output_dir), resolution=resolution) + opts = RtcOptions( + dem_path=str(dem_path), + output_dir=str(output_dir), + resolution=resolution, + apply_bistatic_delay=True, + apply_static_tropo=True, + ) geogrid = generate_geogrids(burst, opts.resolution) - run_single_job(granule, burst, geogrid, opts) + capella_rtc(burst, geogrid, opts) + # run_single_job(granule, burst, geogrid, opts) def opera_rtc_umbra_sicd(granule: str, resolution: int = 30, work_dir: Optional[Path] = None) -> None: @@ -80,7 +87,13 @@ def opera_rtc_capella_sicd(granule: str, resolution: int = 30, work_dir: Optiona raise FileNotFoundError(f'Capella SICD must be present in input dir {input_dir} for processing.') [d.mkdir(parents=True, exist_ok=True) for d in [input_dir, output_dir]] capella_sicd, dem_path = prep_capella(granule_path, work_dir=input_dir) - opts = RtcOptions(dem_path=str(dem_path), output_dir=str(output_dir), resolution=resolution) + opts = RtcOptions( + dem_path=str(dem_path), + output_dir=str(output_dir), + resolution=resolution, + apply_bistatic_delay=False, + apply_static_tropo=False, + ) geogrid = generate_geogrids(capella_sicd, opts.resolution) capella_rtc(capella_sicd, geogrid, opts) diff --git a/src/multirtc/prep_burst.py b/src/multirtc/prep_burst.py index dd01e99..b236753 100644 --- a/src/multirtc/prep_burst.py +++ b/src/multirtc/prep_burst.py @@ -3,15 +3,60 @@ from typing import Optional from zipfile import ZipFile +import isce3 import lxml.etree as ET import s1reader from burst2safe.burst2safe import burst2safe from shapely.geometry import Polygon, box from multirtc import dem, orbit +from multirtc.base import SlcTemplate, from_isce_datetime, to_isce_datetime -def get_s1_granule_bbox(granule_path: Path, buffer: float = 0.025) -> box: +class S1BurstSlc(SlcTemplate): + def __init__(self, safe_path, orbit_path, burst_name): + _, burst_id, swath, _, polarization, _ = burst_name.split('_') + burst_id = int(burst_id) + swath_num = int(swath[2]) + bursts = s1reader.load_bursts(str(safe_path), str(orbit_path), swath_num, polarization) + burst = [b for b in bursts if str(b.burst_id).endswith(f'{burst_id}_{swath.lower()}')][0] + del bursts + vrt_path = safe_path.parent / f'{burst_name}.vrt' + burst.slc_to_vrt_file(vrt_path) + self.id = burst_name + self.filepath = vrt_path + self.footprint = burst.border[0] + self.center = burst.center + self.lookside = 'right' + self.wavelength = burst.wavelength + self.polarization = burst.polarization + self.shape = burst.shape + self.range_pixel_spacing = burst.range_pixel_spacing + self.reference_time = from_isce_datetime(burst.orbit.reference_epoch) + self.sensing_start = (burst.sensing_start - self.reference_time).total_seconds() + self.starting_range = burst.starting_range + self.prf = 1 / burst.azimuth_time_interval + self.orbit = burst.orbit + self.doppler_centroid_grid = isce3.core.LUT2d() + self.radar_grid = isce3.product.RadarGridParameters( + sensing_start=self.sensing_start, + wavelength=self.wavelength, + prf=self.prf, + starting_range=self.starting_range, + range_pixel_spacing=self.range_pixel_spacing, + lookside=isce3.core.LookSide.Right, + length=self.shape[0], + width=self.shape[1], + ref_epoch=to_isce_datetime(self.reference_time), + ) + self.first_valid_line = burst.first_valid_line + self.last_valid_line = burst.last_valid_line + self.first_valid_sample = burst.first_valid_sample + self.last_valid_sample = burst.last_valid_sample + self.source = burst + + +def get_s1_granule_bbox(granule_path: Path) -> box: if granule_path.suffix == '.zip': with ZipFile(granule_path, 'r') as z: manifest_path = [x for x in z.namelist() if x.endswith('manifest.safe')][0] @@ -25,11 +70,10 @@ def get_s1_granule_bbox(granule_path: Path, buffer: float = 0.025) -> box: frame_string = frame_element.find('.//{http://www.opengis.net/gml}coordinates').text coord_strings = [pair.split(',') for pair in frame_string.split(' ')] coords = [(float(lon), float(lat)) for lat, lon in coord_strings] - footprint = Polygon(coords).buffer(buffer) - return box(*footprint.bounds) + return Polygon(coords) -def prep_burst(granule: str, work_dir: Optional[Path] = None) -> Path: +def prep_burst(burst_granule: str, work_dir: Optional[Path] = None) -> Path: """Prepare data for burst-based processing. Args: @@ -43,17 +87,15 @@ def prep_burst(granule: str, work_dir: Optional[Path] = None) -> Path: print('Downloading data...') if len(list(work_dir.glob('S1*.zip'))) == 0: - granule_path = burst2safe(granules=[granule], all_anns=True, work_dir=work_dir) + granule_path = burst2safe(granules=[burst_granule], all_anns=True, work_dir=work_dir) make_archive(base_name=str(granule_path.with_suffix('')), format='zip', base_dir=str(granule_path)) - granule = granule_path.with_suffix('').name granule_path = granule_path.with_suffix('.zip') else: granule_path = work_dir / list(work_dir.glob('S1*.zip'))[0].name orbit_path = orbit.get_orbit(granule_path.with_suffix('').name, save_dir=work_dir) + burst_slc = S1BurstSlc(granule_path, orbit_path, burst_granule) dem_path = work_dir / 'dem.tif' - granule_bbox = get_s1_granule_bbox(granule_path) - dem.download_opera_dem_for_footprint(dem_path, granule_bbox) - burst = s1reader.load_bursts(str(granule_path), str(orbit_path), 1, 'VV')[0] - return burst, dem_path + dem.download_opera_dem_for_footprint(dem_path, burst_slc.footprint) + return burst_slc, dem_path diff --git a/src/multirtc/s1burst_corrections.py b/src/multirtc/s1burst_corrections.py index 6e9a51c..966db9f 100644 --- a/src/multirtc/s1burst_corrections.py +++ b/src/multirtc/s1burst_corrections.py @@ -130,7 +130,7 @@ def compute_correction_lut( def apply_slc_corrections( burst: Sentinel1BurstSlc, - path_slc_vrt: str, + path_input_slc_vrt: str, path_slc_out: str, flag_output_complex: bool = False, flag_thermal_correction: bool = True, @@ -143,7 +143,7 @@ def apply_slc_corrections( ---------- burst_in: Sentinel1BurstSlc Input burst to apply the correction - path_slc_vrt: str + path_input_slc_vrt: str Path to the input burst to apply correction path_slc_out: str Path to the output SLC which the corrections are applied @@ -157,8 +157,7 @@ def apply_slc_corrections( """ # Load the SLC of the burst - burst.slc_to_vrt_file(path_slc_vrt) - slc_gdal_ds = gdal.Open(path_slc_vrt) + slc_gdal_ds = gdal.Open(path_input_slc_vrt) arr_slc_from = slc_gdal_ds.ReadAsArray() # Apply thermal noise correction diff --git a/src/multirtc/sicd.py b/src/multirtc/sicd.py index a2ffe91..685e45e 100644 --- a/src/multirtc/sicd.py +++ b/src/multirtc/sicd.py @@ -1,4 +1,4 @@ -from datetime import datetime, timedelta +from datetime import timedelta from pathlib import Path import isce3 @@ -8,23 +8,14 @@ from sarpy.io.complex.sicd import SICDReader from shapely.geometry import Point, Polygon -from multirtc.base import SlcTemplate - - -def to_isce_datetime(dt): - if isinstance(dt, datetime): - return isce3.core.DateTime(dt) - elif isinstance(dt, np.datetime64): - return isce3.core.DateTime(dt.item()) - else: - raise ValueError(f'Unsupported datetime type: {type(dt)}. Expected datetime or np.datetime64.') +from multirtc.base import SlcTemplate, to_isce_datetime def check_poly_order(poly): assert len(poly.Coefs) == poly.order1 + 1, 'Polynomial order does not match number of coefficients' -class Sicd: +class SicdSlc: def __init__(self, sicd_path): reader = SICDReader(str(sicd_path.expanduser().resolve())) sicd = reader.get_sicds_as_tuple()[0] @@ -125,7 +116,7 @@ def create_complex_beta0(self, outpath, isce_format=True): ds = None -class SicdRzdSlc(SlcTemplate, Sicd): +class SicdRzdSlc(SlcTemplate, SicdSlc): def __init__(self, sicd_path): super().__init__(sicd_path) assert self.source.Grid.Type == 'RGZERO', 'Only range zero doppler grids supported for Capella data' @@ -163,8 +154,3 @@ def get_radar_grid(self): ref_epoch=to_isce_datetime(self.reference_time), ) return radar_grid - - -if __name__ == '__main__': - slc = SicdRzdSlc(Path('~/Data/capella/input/CAPELLA_C15_SS_SICD_VV_20250225215110_20250225215124.ntf')) - breakpoint() From ef2350d574f703541d34b09a08989f371b22923c Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Wed, 28 May 2025 09:15:49 -0500 Subject: [PATCH 11/12] refactor and simplify --- src/multirtc/create_rtc.py | 294 ++-------------------------- src/multirtc/multirtc.py | 47 +++-- src/multirtc/prep_burst.py | 68 +++++-- src/multirtc/s1burst_corrections.py | 66 ------- 4 files changed, 94 insertions(+), 381 deletions(-) diff --git a/src/multirtc/create_rtc.py b/src/multirtc/create_rtc.py index 25f770d..b165c2e 100755 --- a/src/multirtc/create_rtc.py +++ b/src/multirtc/create_rtc.py @@ -7,13 +7,11 @@ import numpy as np import pyproj from osgeo import gdal -from s1reader.s1_burst_slc import Sentinel1BurstSlc from scipy import ndimage from tqdm import tqdm from multirtc.prep_burst import S1BurstSlc -from multirtc.rtc_options import RtcOptions -from multirtc.s1burst_corrections import apply_slc_corrections, compute_correction_lut +from multirtc.s1burst_corrections import compute_correction_lut from multirtc.sicd import SicdSlc @@ -324,239 +322,11 @@ def save_intermediate_geocode_files( del obj -def run_single_job(product_id: str, burst: Sentinel1BurstSlc, geogrid, opts: RtcOptions): - """ - Run geocode burst workflow with user-defined - args stored in dictionary runconfig `cfg` - - Parameters - --------- - cfg: RunConfig - RunConfig object with user runconfig options - """ - # Common initializations - t_start = time.time() - output_dir = str(opts.output_dir) - os.makedirs(output_dir, exist_ok=True) - - raster_format = 'GTiff' - raster_extension = 'tif' - - # Filenames - temp_slc_path = os.path.join(output_dir, 'slc.vrt') - temp_slc_corrected_path = os.path.join(output_dir, f'slc_corrected.{raster_extension}') - geo_burst_filename = f'{output_dir}/{product_id}.{raster_extension}' - nlooks_file = f'{output_dir}/{product_id}_{LAYER_NAME_NUMBER_OF_LOOKS}.{raster_extension}' - rtc_anf_file = f'{output_dir}/{product_id}_{opts.layer_name_rtc_anf}.{raster_extension}' - rtc_anf_gamma0_to_sigma0_file = ( - f'{output_dir}/{product_id}_{LAYER_NAME_RTC_ANF_GAMMA0_TO_SIGMA0}.{raster_extension}' - ) - radar_grid = burst.as_isce3_radargrid() - orbit = burst.orbit - wavelength = burst.wavelength - lookside = radar_grid.lookside - - dem_raster = isce3.io.Raster(opts.dem_path) - ellipsoid = isce3.core.Ellipsoid() - zero_doppler = isce3.core.LUT2d() - exponent = 1 if (opts.apply_thermal_noise or opts.apply_ads_rad) else 2 - - x_snap = geogrid.spacing_x - y_snap = geogrid.spacing_y - geogrid.start_x = np.floor(float(geogrid.start_x) / x_snap) * x_snap - geogrid.start_y = np.ceil(float(geogrid.start_y) / y_snap) * y_snap - - # Convert to beta0 and apply thermal noise correction - if opts.apply_thermal_noise or opts.apply_abs_rad: - apply_slc_corrections( - burst, - temp_slc_path, - temp_slc_corrected_path, - flag_output_complex=False, - flag_thermal_correction=opts.apply_thermal_noise, - flag_apply_abs_rad_correction=opts.apply_abs_rad, - ) - input_burst_filename = temp_slc_corrected_path - else: - input_burst_filename = temp_slc_path - - # geocoding optional arguments - geocode_kwargs = {} - layover_shadow_mask_geocode_kwargs = {} - - # get sub_swaths metadata - if opts.apply_valid_samples_sub_swath_masking: - # Extract burst boundaries and create sub_swaths object to mask - # invalid radar samples - n_subswaths = 1 - sub_swaths = isce3.product.SubSwaths(radar_grid.length, radar_grid.width, n_subswaths) - last_range_sample = min([burst.last_valid_sample, radar_grid.width]) - valid_samples_sub_swath = np.repeat( - [[burst.first_valid_sample, last_range_sample + 1]], radar_grid.length, axis=0 - ) - for i in range(burst.first_valid_line): - valid_samples_sub_swath[i, :] = 0 - for i in range(burst.last_valid_line, radar_grid.length): - valid_samples_sub_swath[i, :] = 0 - - sub_swaths.set_valid_samples_array(1, valid_samples_sub_swath) - geocode_kwargs['sub_swaths'] = sub_swaths - layover_shadow_mask_geocode_kwargs['sub_swaths'] = sub_swaths - - # Calculate layover/shadow mask - layover_shadow_mask_file = f'{output_dir}/{product_id}_{LAYER_NAME_LAYOVER_SHADOW_MASK}.{raster_extension}' - logger.info(f' computing layover shadow mask for {product_id}') - radar_grid_layover_shadow_mask = radar_grid - slantrange_layover_shadow_mask_raster = compute_layover_shadow_mask( - radar_grid_layover_shadow_mask, - orbit, - geogrid, - dem_raster, - layover_shadow_mask_file, - raster_format, - output_dir, - shadow_dilation_size=opts.shadow_dilation_size, - threshold_rdr2geo=opts.rdr2geo_threshold, - numiter_rdr2geo=opts.rdr2geo_numiter, - threshold_geo2rdr=opts.geo2rdr_threshold, - numiter_geo2rdr=opts.geo2rdr_numiter, - memory_mode=opts.memory_mode_isce3, - geocode_options=layover_shadow_mask_geocode_kwargs, - ) - logger.info(f'file saved: {layover_shadow_mask_file}') - if opts.apply_shadow_masking: - geocode_kwargs['input_layover_shadow_mask_raster'] = slantrange_layover_shadow_mask_raster - - out_geo_nlooks_obj = isce3.io.Raster(nlooks_file, geogrid.width, geogrid.length, 1, gdal.GDT_Float32, raster_format) - out_geo_rtc_obj = isce3.io.Raster(rtc_anf_file, geogrid.width, geogrid.length, 1, gdal.GDT_Float32, raster_format) - out_geo_rtc_gamma0_to_sigma0_obj = isce3.io.Raster( - rtc_anf_gamma0_to_sigma0_file, - geogrid.width, - geogrid.length, - 1, - gdal.GDT_Float32, - raster_format, - ) - geocode_kwargs['out_geo_rtc_gamma0_to_sigma0'] = out_geo_rtc_gamma0_to_sigma0_obj - - # Calculate geolocation correction LUT - if opts.apply_bistatic_delay or opts.apply_static_tropo: - rg_lut, az_lut = compute_correction_lut( - burst, - dem_raster, - output_dir, - opts.correction_lut_range_spacing_in_meters, - opts.correction_lut_azimuth_spacing_in_meters, - opts.apply_bistatic_delay, - opts.apply_static_tropo, - ) - geocode_kwargs['az_time_correction'] = az_lut - if rg_lut is not None: - geocode_kwargs['slant_range_correction'] = rg_lut - - rdr_burst_raster = isce3.io.Raster(input_burst_filename) - # Generate output geocoded burst raster - geo_burst_raster = isce3.io.Raster( - geo_burst_filename, geogrid.width, geogrid.length, rdr_burst_raster.num_bands, gdal.GDT_Float32, raster_format - ) - - # init Geocode object depending on raster type - if rdr_burst_raster.datatype() == gdal.GDT_Float32: - geo_obj = isce3.geocode.GeocodeFloat32() - elif rdr_burst_raster.datatype() == gdal.GDT_Float64: - geo_obj = isce3.geocode.GeocodeFloat64() - elif rdr_burst_raster.datatype() == gdal.GDT_CFloat32: - geo_obj = isce3.geocode.GeocodeCFloat32() - elif rdr_burst_raster.datatype() == gdal.GDT_CFloat64: - geo_obj = isce3.geocode.GeocodeCFloat64() - else: - err_str = 'Unsupported raster type for geocoding' - raise NotImplementedError(err_str) - - # init geocode members - geo_obj.orbit = orbit - geo_obj.ellipsoid = ellipsoid - geo_obj.doppler = zero_doppler - geo_obj.threshold_geo2rdr = opts.geo2rdr_threshold - geo_obj.numiter_geo2rdr = opts.geo2rdr_numiter - - # set data interpolator based on the geocode algorithm - if opts.geocode_algorithm_isce3 == isce3.geocode.GeocodeOutputMode.INTERP: - geo_obj.data_interpolator = opts.geocode_algorithm_isce3 - - geo_obj.geogrid( - geogrid.start_x, - geogrid.start_y, - geogrid.spacing_x, - geogrid.spacing_y, - geogrid.width, - geogrid.length, - geogrid.epsg, - ) - - geo_obj.geocode( - radar_grid=radar_grid, - input_raster=rdr_burst_raster, - output_raster=geo_burst_raster, - dem_raster=dem_raster, - output_mode=opts.geocode_algorithm_isce3, - geogrid_upsampling=opts.geogrid_upsampling, - flag_apply_rtc=opts.apply_rtc, - input_terrain_radiometry=opts.input_terrain_radiometry_isce3, - output_terrain_radiometry=opts.terrain_radiometry_isce3, - exponent=exponent, - rtc_min_value_db=opts.rtc_min_value_db, - rtc_upsampling=opts.rtc_upsampling, - rtc_algorithm=opts.rtc_algorithm_isce3, - abs_cal_factor=opts.abs_cal_factor, - flag_upsample_radar_grid=opts.upsample_radar_grid, - clip_min=opts.clip_min, - clip_max=opts.clip_max, - out_geo_nlooks=out_geo_nlooks_obj, - out_geo_rtc=out_geo_rtc_obj, - rtc_area_beta_mode=opts.rtc_area_beta_mode_isce3, - # out_geo_rtc_gamma0_to_sigma0=out_geo_rtc_gamma0_to_sigma0_obj, - input_rtc=None, - output_rtc=None, - dem_interp_method=opts.dem_interpolation_method_isce3, - memory_mode=opts.memory_mode_isce3, - **geocode_kwargs, - ) - - del geo_burst_raster - - out_geo_nlooks_obj.close_dataset() - del out_geo_nlooks_obj - - out_geo_rtc_obj.close_dataset() - del out_geo_rtc_obj - - out_geo_rtc_gamma0_to_sigma0_obj.close_dataset() - del out_geo_rtc_gamma0_to_sigma0_obj - - radar_grid_file_dict = {} - save_intermediate_geocode_files( - geogrid, - opts.dem_interpolation_method_isce3, - product_id, - output_dir, - raster_extension, - dem_raster, - radar_grid_file_dict, - lookside, - wavelength, - orbit, - ) - - t_end = time.time() - logger.info(f'elapsed time: {t_end - t_start}') - - -def capella_rtc(sicd, geogrid, opts): +def rtc(slc, geogrid, opts): # Common initializations t_start = time.time() output_dir = str(opts.output_dir) - product_id = sicd.id + product_id = slc.id os.makedirs(output_dir, exist_ok=True) raster_format = 'GTiff' @@ -569,14 +339,14 @@ def capella_rtc(sicd, geogrid, opts): rtc_anf_gamma0_to_sigma0_file = ( f'{output_dir}/{product_id}_{LAYER_NAME_RTC_ANF_GAMMA0_TO_SIGMA0}.{raster_extension}' ) - radar_grid = sicd.radar_grid - orbit = sicd.orbit - wavelength = sicd.wavelength + radar_grid = slc.radar_grid + orbit = slc.orbit + wavelength = slc.wavelength lookside = radar_grid.lookside dem_raster = isce3.io.Raster(opts.dem_path) ellipsoid = isce3.core.Ellipsoid() - doppler = sicd.doppler_centroid_grid + doppler = slc.doppler_centroid_grid exponent = 2 x_snap = geogrid.spacing_x @@ -584,47 +354,25 @@ def capella_rtc(sicd, geogrid, opts): geogrid.start_x = np.floor(float(geogrid.start_x) / x_snap) * x_snap geogrid.start_y = np.ceil(float(geogrid.start_y) / y_snap) * y_snap - if isinstance(sicd, SicdSlc): - input_filename = sicd.filepath.parent / (sicd.filepath.stem + '_beta0.tif') - if not input_filename.exists(): - sicd.create_complex_beta0(input_filename) - input_filename = str(input_filename) - elif isinstance(sicd, S1BurstSlc) and opts.apply_thermal_noise or opts.apply_abs_rad: - # Convert to beta0 and apply thermal noise correction - input_filename = str(sicd.filepath.parent / (sicd.filepath.stem + '_beta0.tif')) - apply_slc_corrections( - sicd.source, - sicd.filepath, - input_filename, - flag_output_complex=True, - flag_thermal_correction=opts.apply_thermal_noise, - flag_apply_abs_rad_correction=opts.apply_abs_rad, - ) - else: - input_filename = str(sicd.filepath) - # geocoding optional arguments geocode_kwargs = {} layover_shadow_mask_geocode_kwargs = {} - # get sub_swaths metadata - if isinstance(sicd, S1BurstSlc) and opts.apply_valid_samples_sub_swath_masking: - # Extract burst boundaries and create sub_swaths object to mask - # invalid radar samples - n_subswaths = 1 - sub_swaths = isce3.product.SubSwaths(radar_grid.length, radar_grid.width, n_subswaths) - last_range_sample = min([sicd.last_valid_sample, radar_grid.width]) - valid_samples_sub_swath = np.repeat( - [[sicd.first_valid_sample, last_range_sample + 1]], radar_grid.length, axis=0 - ) - for i in range(sicd.first_valid_line): - valid_samples_sub_swath[i, :] = 0 - for i in range(sicd.last_valid_line, radar_grid.length): - valid_samples_sub_swath[i, :] = 0 - - sub_swaths.set_valid_samples_array(1, valid_samples_sub_swath) + if isinstance(slc, SicdSlc): + input_filename = slc.filepath.parent / (slc.filepath.stem + '_beta0.tif') + if not input_filename.exists(): + slc.create_complex_beta0(input_filename) + input_filename = str(input_filename) + elif isinstance(slc, S1BurstSlc): + input_filename = slc.filepath.parent / (slc.filepath.stem + '_beta0.tif') + if not input_filename.exists(): + slc.create_complex_beta0(input_filename, flag_thermal_correction=opts.apply_thermal_noise) + input_filename = str(input_filename) + sub_swaths = slc.apply_valid_data_masking() geocode_kwargs['sub_swaths'] = sub_swaths layover_shadow_mask_geocode_kwargs['sub_swaths'] = sub_swaths + else: + input_filename = str(slc.filepath) layover_shadow_mask_file = f'{output_dir}/{product_id}_{LAYER_NAME_LAYOVER_SHADOW_MASK}.{raster_extension}' logger.info(f' computing layover shadow mask for {product_id}') @@ -658,7 +406,7 @@ def capella_rtc(sicd, geogrid, opts): geocode_kwargs['out_geo_rtc_gamma0_to_sigma0'] = out_geo_rtc_gamma0_to_sigma0_obj if opts.apply_bistatic_delay or opts.apply_static_tropo: rg_lut, az_lut = compute_correction_lut( - sicd.source, + slc.source, dem_raster, output_dir, opts.correction_lut_range_spacing_in_meters, diff --git a/src/multirtc/multirtc.py b/src/multirtc/multirtc.py index 3fe8f29..4ae8f3d 100644 --- a/src/multirtc/multirtc.py +++ b/src/multirtc/multirtc.py @@ -4,7 +4,7 @@ import isce3 -from multirtc.create_rtc import capella_rtc, run_single_job, umbra_rtc +from multirtc.create_rtc import rtc, umbra_rtc from multirtc.define_geogrid import generate_geogrids from multirtc.prep_burst import prep_burst from multirtc.prep_capella import prep_capella @@ -45,15 +45,14 @@ def opera_rtc_s1_burst(granule: str, resolution: int = 30, work_dir: Optional[Pa apply_static_tropo=True, ) geogrid = generate_geogrids(burst, opts.resolution) - capella_rtc(burst, geogrid, opts) - # run_single_job(granule, burst, geogrid, opts) + rtc(burst, geogrid, opts) -def opera_rtc_umbra_sicd(granule: str, resolution: int = 30, work_dir: Optional[Path] = None) -> None: - """Create an OPERA RTC for an UMBRA SICD file +def opera_rtc_capella_sicd(granule: str, resolution: int = 30, work_dir: Optional[Path] = None) -> None: + """Create an OPERA RTC for an CAPELLA SICD file Args: - granule: Umbra SICD file name to create an RTC for + granule: Capella SICD file name to create an RTC for resolution: Resolution of the output RTC (m) work_dir: Working directory for processing """ @@ -63,18 +62,25 @@ def opera_rtc_umbra_sicd(granule: str, resolution: int = 30, work_dir: Optional[ output_dir = work_dir / 'output' granule_path = input_dir / granule if not granule_path.exists(): - raise FileNotFoundError(f'Umbra SICD must be present in input dir {input_dir} for processing.') + raise FileNotFoundError(f'Capella SICD must be present in input dir {input_dir} for processing.') [d.mkdir(parents=True, exist_ok=True) for d in [input_dir, output_dir]] - umbra_sicd, dem_path = prep_umbra(granule_path, work_dir=input_dir) - geogrid = generate_geogrids(umbra_sicd, resolution, rda=False) - umbra_rtc(umbra_sicd, geogrid, dem_path, output_dir=output_dir) + capella_sicd, dem_path = prep_capella(granule_path, work_dir=input_dir) + opts = RtcOptions( + dem_path=str(dem_path), + output_dir=str(output_dir), + resolution=resolution, + apply_bistatic_delay=False, + apply_static_tropo=False, + ) + geogrid = generate_geogrids(capella_sicd, opts.resolution) + rtc(capella_sicd, geogrid, opts) -def opera_rtc_capella_sicd(granule: str, resolution: int = 30, work_dir: Optional[Path] = None) -> None: - """Create an OPERA RTC for an CAPELLA SICD file +def opera_rtc_umbra_sicd(granule: str, resolution: int = 30, work_dir: Optional[Path] = None) -> None: + """Create an OPERA RTC for an UMBRA SICD file Args: - granule: Capella SICD file name to create an RTC for + granule: Umbra SICD file name to create an RTC for resolution: Resolution of the output RTC (m) work_dir: Working directory for processing """ @@ -84,18 +90,11 @@ def opera_rtc_capella_sicd(granule: str, resolution: int = 30, work_dir: Optiona output_dir = work_dir / 'output' granule_path = input_dir / granule if not granule_path.exists(): - raise FileNotFoundError(f'Capella SICD must be present in input dir {input_dir} for processing.') + raise FileNotFoundError(f'Umbra SICD must be present in input dir {input_dir} for processing.') [d.mkdir(parents=True, exist_ok=True) for d in [input_dir, output_dir]] - capella_sicd, dem_path = prep_capella(granule_path, work_dir=input_dir) - opts = RtcOptions( - dem_path=str(dem_path), - output_dir=str(output_dir), - resolution=resolution, - apply_bistatic_delay=False, - apply_static_tropo=False, - ) - geogrid = generate_geogrids(capella_sicd, opts.resolution) - capella_rtc(capella_sicd, geogrid, opts) + umbra_sicd, dem_path = prep_umbra(granule_path, work_dir=input_dir) + geogrid = generate_geogrids(umbra_sicd, resolution, rda=False) + umbra_rtc(umbra_sicd, geogrid, dem_path, output_dir=output_dir) def main(): diff --git a/src/multirtc/prep_burst.py b/src/multirtc/prep_burst.py index b236753..638446e 100644 --- a/src/multirtc/prep_burst.py +++ b/src/multirtc/prep_burst.py @@ -1,18 +1,20 @@ from pathlib import Path from shutil import make_archive from typing import Optional -from zipfile import ZipFile import isce3 -import lxml.etree as ET +import numpy as np import s1reader from burst2safe.burst2safe import burst2safe -from shapely.geometry import Polygon, box +from osgeo import gdal from multirtc import dem, orbit from multirtc.base import SlcTemplate, from_isce_datetime, to_isce_datetime +gdal.UseExceptions() + + class S1BurstSlc(SlcTemplate): def __init__(self, safe_path, orbit_path, burst_name): _, burst_id, swath, _, polarization, _ = burst_name.split('_') @@ -55,22 +57,52 @@ def __init__(self, safe_path, orbit_path, burst_name): self.last_valid_sample = burst.last_valid_sample self.source = burst + def apply_valid_data_masking(self): + # Extract burst boundaries and create sub_swaths object to mask invalid radar samples + n_subswaths = 1 + sub_swaths = isce3.product.SubSwaths(self.radar_grid.length, self.radar_grid.width, n_subswaths) + last_range_sample = min([self.last_valid_sample, self.radar_grid.width]) + valid_samples_sub_swath = np.repeat( + [[self.first_valid_sample, last_range_sample + 1]], self.radar_grid.length, axis=0 + ) + for i in range(self.first_valid_line): + valid_samples_sub_swath[i, :] = 0 + for i in range(self.last_valid_line, self.radar_grid.length): + valid_samples_sub_swath[i, :] = 0 -def get_s1_granule_bbox(granule_path: Path) -> box: - if granule_path.suffix == '.zip': - with ZipFile(granule_path, 'r') as z: - manifest_path = [x for x in z.namelist() if x.endswith('manifest.safe')][0] - with z.open(manifest_path) as m: - manifest = ET.parse(m).getroot() - else: - manifest_path = granule_path / 'manifest.safe' - manifest = ET.parse(manifest_path).getroot() - - frame_element = [x for x in manifest.findall('.//metadataObject') if x.get('ID') == 'measurementFrameSet'][0] - frame_string = frame_element.find('.//{http://www.opengis.net/gml}coordinates').text - coord_strings = [pair.split(',') for pair in frame_string.split(' ')] - coords = [(float(lon), float(lat)) for lat, lon in coord_strings] - return Polygon(coords) + sub_swaths.set_valid_samples_array(1, valid_samples_sub_swath) + return sub_swaths + + def create_complex_beta0(self, outpath: Path, flag_thermal_correction: bool = True): + """Apply conversion to beta0 and optionally applies a thermal correction.""" + # Load the SLC of the burst + slc_gdal_ds = gdal.Open(str(self.filepath)) + arr_slc_from = slc_gdal_ds.ReadAsArray() + + # Apply thermal noise correction + if flag_thermal_correction: + corrected_image = np.abs(arr_slc_from) ** 2 - self.source.thermal_noise_lut + min_backscatter = 0 + max_backscatter = None + corrected_image = np.clip(corrected_image, min_backscatter, max_backscatter) + else: + corrected_image = np.abs(arr_slc_from) ** 2 + + # Apply absolute radiometric correction + corrected_image = corrected_image / self.source.burst_calibration.beta_naught**2 + + factor_mag = np.sqrt(corrected_image) / np.abs(arr_slc_from) + factor_mag[np.isnan(factor_mag)] = 0.0 + corrected_image = arr_slc_from * factor_mag + dtype = gdal.GDT_CFloat32 + + # Save the corrected image + drvout = gdal.GetDriverByName('GTiff') + raster_out = drvout.Create(outpath, self.shape[1], self.shape[0], 1, dtype) + band_out = raster_out.GetRasterBand(1) + band_out.WriteArray(corrected_image) + band_out.FlushCache() + del band_out def prep_burst(burst_granule: str, work_dir: Optional[Path] = None) -> Path: diff --git a/src/multirtc/s1burst_corrections.py b/src/multirtc/s1burst_corrections.py index 966db9f..d4a9b9b 100644 --- a/src/multirtc/s1burst_corrections.py +++ b/src/multirtc/s1burst_corrections.py @@ -3,7 +3,6 @@ import isce3 import numpy as np from osgeo import gdal -from s1reader.s1_burst_slc import Sentinel1BurstSlc logger = logging.getLogger('rtc_s1') @@ -126,68 +125,3 @@ def compute_correction_lut( ) return rg_lut, az_lut - - -def apply_slc_corrections( - burst: Sentinel1BurstSlc, - path_input_slc_vrt: str, - path_slc_out: str, - flag_output_complex: bool = False, - flag_thermal_correction: bool = True, - flag_apply_abs_rad_correction: bool = True, -): - """Apply thermal correction stored in burst_in. Save the corrected signal - back to ENVI format. Preserves the phase when the output is complex - - Parameters - ---------- - burst_in: Sentinel1BurstSlc - Input burst to apply the correction - path_input_slc_vrt: str - Path to the input burst to apply correction - path_slc_out: str - Path to the output SLC which the corrections are applied - flag_output_complex: bool - `path_slc_out` will be in complex number when this is `True` - Otherwise, the output will be amplitude only. - flag_thermal_correction: bool - flag whether or not to apple the thermal correction. - flag_apply_abs_rad_correction: bool - Flag to apply radiometric calibration - """ - - # Load the SLC of the burst - slc_gdal_ds = gdal.Open(path_input_slc_vrt) - arr_slc_from = slc_gdal_ds.ReadAsArray() - - # Apply thermal noise correction - if flag_thermal_correction: - logger.info(' applying thermal noise correction to burst SLC') - corrected_image = np.abs(arr_slc_from) ** 2 - burst.thermal_noise_lut - min_backscatter = 0 - max_backscatter = None - corrected_image = np.clip(corrected_image, min_backscatter, max_backscatter) - else: - corrected_image = np.abs(arr_slc_from) ** 2 - - # Apply absolute radiometric correction - if flag_apply_abs_rad_correction: - logger.info(' applying absolute radiometric correction to burst SLC') - corrected_image = corrected_image / burst.burst_calibration.beta_naught**2 - - # Output as complex - if flag_output_complex: - factor_mag = np.sqrt(corrected_image) / np.abs(arr_slc_from) - factor_mag[np.isnan(factor_mag)] = 0.0 - corrected_image = arr_slc_from * factor_mag - dtype = gdal.GDT_CFloat32 - else: - dtype = gdal.GDT_Float32 - - # Save the corrected image - drvout = gdal.GetDriverByName('GTiff') - raster_out = drvout.Create(path_slc_out, burst.shape[1], burst.shape[0], 1, dtype) - band_out = raster_out.GetRasterBand(1) - band_out.WriteArray(corrected_image) - band_out.FlushCache() - del band_out From da181a264c1c8b724763d356d8e0379504ee9323 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Wed, 28 May 2025 09:27:51 -0500 Subject: [PATCH 12/12] update docs --- CHANGELOG.md | 5 +++++ README.md | 16 +++++++++++----- 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d87058b..28c9aff 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,11 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [PEP 440](https://www.python.org/dev/peps/pep-0440/) and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [0.2.0] + +### Added +* Support for Capella SICD SLC products + ## [0.1.1] ### Added diff --git a/README.md b/README.md index 45bd5b9..091c8ec 100644 --- a/README.md +++ b/README.md @@ -1,13 +1,19 @@ # MultiRTC -A python library for creating ISCE3-based RTCs from Sentinel-1 Burst and Umbra SICD SLC data. +A python library for creating ISCE3-based RTCs for multiple SAR data sources -**ALL CREDIT FOR THIS LIBRARY'S RTC ALGORITHM GOES TO GUSTAVO SHIROMA AND THE [JPL OPERA TEAM](https://www.jpl.nasa.gov/go/opera). THIS PLUGIN MERELY ALLOWS OTHERS TO USE THEIR ALGORITHM WITH MULTIPLE SENSORS.** +**ALL CREDIT FOR THIS LIBRARY'S RTC ALGORITHM GOES TO GUSTAVO SHIROMA AND THE JPL [OPERA](https://www.jpl.nasa.gov/go/opera/about-opera/) AND [ISCE3](https://github.com/isce-framework/isce3) TEAMS. THIS PLUGIN MERELY ALLOWS OTHERS TO USE THEIR ALGORITHM WITH MULTIPLE SENSORS.** + +The RTC algorithm utilized by this library is described in [Shiroma et al., 2023](https://doi.org/10.1109/TGRS.2022.3147472). ## Usage MultiRTC allows users to create RTC products from SLC data for multiple SAR sensor platforms. Currently this list includes: -- [Sentinel-1 burst SLCs](https://www.earthdata.nasa.gov/data/catalog/alaska-satellite-facility-distributed-active-archive-center-sentinel-1-bursts-version) +Full RTC: +- [Sentinel-1 Burst SLCs](https://www.earthdata.nasa.gov/data/catalog/alaska-satellite-facility-distributed-active-archive-center-sentinel-1-bursts-version) +- [Capella SICD SLCs](https://www.capellaspace.com/earth-observation/data) + +Geocode Only: - [UMBRA SICD SLCs](https://help.umbra.space/product-guide/umbra-products/umbra-product-specifications) To create an RTC, use the `multirtc` CLI entrypoint using the following pattern: @@ -15,7 +21,7 @@ To create an RTC, use the `multirtc` CLI entrypoint using the following pattern: ```bash multirtc PLATFORM SLC-GRANULE --resolution RESOLUTION --work-dir WORK-DIR ``` -Where `PLATFORM` is the name of the satellite platform (currently `S1` or `UMBRA`), `SLC-GRANULE` is the name of the SLC granule, `RESOLUTION` is the desired output resolution of the RTC image in meters, and `WORK-DIR` is the name of the working directory to perform processing in. Inputs such as the SLC data, DEM, and external orbit information are stored in `WORK-DIR/input`, while the RTC image and associated outputs are stored in `WORK-DIR/output` once processing is complete. SLC data that is available in the [Alaska Satellite Facility's data archive](https://search.asf.alaska.edu/#/?maxResults=250) (such as Sentinel-1 burst SLCs) will be automatically downloaded to the input directory, but data not available in this archive (Umbra SICD SLCs) are required to be staged in the input directory prior to processing. +Where `PLATFORM` is the name of the satellite platform (currently `S1`, `CAPELLA` or `UMBRA`), `SLC-GRANULE` is the name of the SLC granule, `RESOLUTION` is the desired output resolution of the RTC image in meters, and `WORK-DIR` is the name of the working directory to perform processing in. Inputs such as the SLC data, DEM, and external orbit information are stored in `WORK-DIR/input`, while the RTC image and associated outputs are stored in `WORK-DIR/output` once processing is complete. SLC data that is available in the [Alaska Satellite Facility's data archive](https://search.asf.alaska.edu/#/?maxResults=250) (such as Sentinel-1 Burst SLCs) will be automatically downloaded to the input directory, but data not available in this archive (commercial datasets) are required to be staged in the input directory prior to processing. Output RTC products are in Gamma0 radiometry. @@ -23,7 +29,7 @@ Output RTC products are in Gamma0 radiometry. Currently, the Umbra processor only supports basic geocoding and not full RTC processing. ISCE3's RTC algorithm is only designed to work with Range Migration Algorithm (RMA) focused SLC products, but Umbra creates their data using the Polar Format Algorithm (PFA). Using an [approach detailed by Piyush Agram](https://arxiv.org/abs/2503.07889v1) to adapt RMA approaches to the PFA image geometry, we have developed a workflow to geocode an Umbra SLC but there is more work to be done to implement full RTC processing. Since full RTC is not yet implemented, Umbra geocoded products are in Sigma0 radiometry. ### DEM options -Currently, only the NISAR DEM is supported. This is a roughly global Height Above Ellipsoid DEM sourced from the [COP-30 DEM](https://portal.opentopography.org/raster?opentopoID=OTSDEM.032021.4326.3). In the future, we hope to support a wider variety of automatically retrieved and user provided DEMs. +Currently, only the OPERA DEM is supported. This is a global Height Above Ellipsoid DEM sourced from the [COP-30 DEM](https://portal.opentopography.org/raster?opentopoID=OTSDEM.032021.4326.3). In the future, we hope to support a wider variety of automatically retrieved and user provided DEMs. ## Developer Setup 1. Ensure that conda is installed on your system (we recommend using [mambaforge](https://github.com/conda-forge/miniforge#mambaforge) to reduce setup times).