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: >- 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 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). diff --git a/src/multirtc/base.py b/src/multirtc/base.py new file mode 100644 index 0000000..e1257e4 --- /dev/null +++ b/src/multirtc/base.py @@ -0,0 +1,58 @@ +from abc import ABC +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, + '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 862002e..b165c2e 100755 --- a/src/multirtc/create_rtc.py +++ b/src/multirtc/create_rtc.py @@ -7,12 +7,12 @@ 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.rtc_options import RtcOptions -from multirtc.s1burst_corrections import apply_slc_corrections, compute_correction_lut +from multirtc.prep_burst import S1BurstSlc +from multirtc.s1burst_corrections import compute_correction_lut +from multirtc.sicd import SicdSlc logger = logging.getLogger('rtc_s1') @@ -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) @@ -329,86 +322,58 @@ 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 - """ +def rtc(slc, geogrid, opts): # Common initializations t_start = time.time() output_dir = str(opts.output_dir) + product_id = slc.id 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}' + geo_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 + 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() - zero_doppler = isce3.core.LUT2d() - exponent = 1 if (opts.apply_thermal_noise or opts.apply_ads_rad) else 2 + doppler = slc.doppler_centroid_grid + exponent = 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) + 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) - # 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 @@ -416,7 +381,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, @@ -428,6 +392,7 @@ def run_single_job(product_id: str, burst: Sentinel1BurstSlc, geogrid, opts: Rtc numiter_geo2rdr=opts.geo2rdr_numiter, memory_mode=opts.memory_mode_isce3, geocode_options=layover_shadow_mask_geocode_kwargs, + doppler=doppler, ) logger.info(f'file saved: {layover_shadow_mask_file}') if opts.apply_shadow_masking: @@ -436,19 +401,12 @@ def run_single_job(product_id: str, burst: Sentinel1BurstSlc, geogrid, opts: Rtc 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, + 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, + slc.source, dem_raster, output_dir, opts.correction_lut_range_spacing_in_meters, @@ -460,176 +418,6 @@ def run_single_job(product_id: str, burst: Sentinel1BurstSlc, geogrid, opts: Rtc 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 umbra_rtc_with_radargrid(umbra_sicd, geogrid, opts): - # Common initializations - t_start = time.time() - output_dir = str(opts.output_dir) - product_id = umbra_sicd.id - os.makedirs(output_dir, exist_ok=True) - - raster_format = 'GTiff' - raster_extension = 'tif' - - # Filenames - geo_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 = umbra_sicd.as_isce3_radargrid() - orbit = umbra_sicd.orbit - wavelength = umbra_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() - exponent = 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 - - # input_filename = save_as_beta0(umbra_sicd, Path(opts.output_dir)) - input_filename = 'beta0.tif' - input_filename = str(input_filename) - - # geocoding optional arguments - geocode_kwargs = {} - layover_shadow_mask_geocode_kwargs = {} - - 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, - umbra_sicd, - 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, - doppler=doppler, - ) - 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 - rdr_raster = isce3.io.Raster(input_filename) # Generate output geocoded burst raster geo_raster = isce3.io.Raster( 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/multirtc.py b/src/multirtc/multirtc.py index 584a442..4ae8f3d 100644 --- a/src/multirtc/multirtc.py +++ b/src/multirtc/multirtc.py @@ -2,13 +2,26 @@ from pathlib import Path from typing import Optional -from multirtc.create_rtc import run_single_job, umbra_rtc +import isce3 + +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 from multirtc.prep_umbra import prep_umbra from multirtc.rtc_options import RtcOptions +def print_wkt(sicd): + 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( + 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 @@ -24,9 +37,43 @@ 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) + rtc(burst, 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 + + 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) + 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_umbra_sicd(granule: str, resolution: int = 30, work_dir: Optional[Path] = None) -> None: @@ -57,7 +104,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 +114,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_burst.py b/src/multirtc/prep_burst.py index dd01e99..638446e 100644 --- a/src/multirtc/prep_burst.py +++ b/src/multirtc/prep_burst.py @@ -1,35 +1,111 @@ from pathlib import Path from shutil import make_archive from typing import Optional -from zipfile import ZipFile -import lxml.etree as ET +import isce3 +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 -def get_s1_granule_bbox(granule_path: Path, buffer: float = 0.025) -> 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() +gdal.UseExceptions() + + +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 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 + + 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 - 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] - footprint = Polygon(coords).buffer(buffer) - return box(*footprint.bounds) + # 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(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 +119,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/prep_capella.py b/src/multirtc/prep_capella.py new file mode 100644 index 0000000..ad8bcfd --- /dev/null +++ b/src/multirtc/prep_capella.py @@ -0,0 +1,25 @@ +from pathlib import Path +from typing import Optional + +from osgeo import gdal + +from multirtc import dem +from multirtc.sicd import SicdRzdSlc + + +gdal.UseExceptions() + + +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() + 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/s1burst_corrections.py b/src/multirtc/s1burst_corrections.py index 6e9a51c..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,69 +125,3 @@ def compute_correction_lut( ) return rg_lut, az_lut - - -def apply_slc_corrections( - burst: Sentinel1BurstSlc, - path_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_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 - burst.slc_to_vrt_file(path_slc_vrt) - slc_gdal_ds = gdal.Open(path_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 diff --git a/src/multirtc/sicd.py b/src/multirtc/sicd.py new file mode 100644 index 0000000..685e45e --- /dev/null +++ b/src/multirtc/sicd.py @@ -0,0 +1,156 @@ +from datetime import 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, to_isce_datetime + + +def check_poly_order(poly): + assert len(poly.Coefs) == poly.order1 + 1, 'Polynomial order does not match number of coefficients' + + +class SicdSlc: + 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, 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' + 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