diff --git a/CHANGELOG.md b/CHANGELOG.md index c88db45..8ae14a5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,12 +8,16 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [0.5.0] +### Added +* Geocode-only entrypoint. + ### Changed * Update readme with docker instructions and layer info. ### Removed * Output of area projection metadata. * Prototype polar grid support in favor of full support via the new docker image. +* Old per-file entrypoints. ### Fixed * Plotting issues with the cal/val submodule. diff --git a/README.md b/README.md index 97a9bb9..5ed8c16 100644 --- a/README.md +++ b/README.md @@ -45,6 +45,14 @@ Where `PLATFORM` is the name of the satellite platform (currently `S1`, `CAPELLA Output RTC pixel values represent gamma0 power. +To create an image that is geocoded but not radiometricly corrected, use the `geocoded` flag instead: + +```bash +multirtc geocode PLATFORM SLC-GRANULE --resolution RESOLUTION --work-dir WORK-DIR +``` + +Output geocoded pixel values represent sigma0 power. + ### Running via Docker In addition to the main python interface, I've also provided an experimental docker container that contains full support for polar grid format SICD data. Encapsulating this functionality in a docker container is ncessary for now because it requires re-compiling a development version of ISCE3. The docker container can be run using a similar interface, with exception of needing to pass your EarthData credentials and the need to pass a mounted volume with an `input` and `output` directory inside: @@ -78,8 +86,10 @@ MultiRTC outputs one main RTC image and seven metadata images as GeoTIFFs. All l More information on the metadata images can be found in the OPERA RTC Static Product guide on the [OPERA RTC Product website](https://www.jpl.nasa.gov/go/opera/products/rtc-product/). +All metadata images other than `FILEID_mask.tif`, and `FILEID_number_of_looks.tif` are omitted for geocode-only products. + ### DEM options -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. +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. If the low resolution of the default DEM is causing radiometry issues, try using the `geocode` instead of `rtc` workflow. ## Calibration & Validation Subcommands diff --git a/src/multirtc/__main__.py b/src/multirtc/__main__.py index e6adad7..efbbdd4 100644 --- a/src/multirtc/__main__.py +++ b/src/multirtc/__main__.py @@ -1,6 +1,6 @@ import argparse -from multirtc import multirtc +from multirtc import geocode, multirtc from multirtc.multimetric import ale, point_target, rle @@ -12,8 +12,11 @@ def main(): ) subparsers = global_parser.add_subparsers(title='command', help='MultiRTC sub-commands') - multirtc_parser = multirtc.create_parser(subparsers.add_parser('rtc', help=multirtc.__doc__)) - multirtc_parser.set_defaults(func=multirtc.run) + rtc_parser = multirtc.create_parser(subparsers.add_parser('rtc', help=multirtc.__doc__)) + rtc_parser.set_defaults(func=multirtc.run) + + geocode_parser = multirtc.create_parser(subparsers.add_parser('geocode', help=geocode.__doc__)) + geocode_parser.set_defaults(func=geocode.run) ale_parser = ale.create_parser(subparsers.add_parser('ale', help=ale.__doc__)) ale_parser.set_defaults(func=ale.run) diff --git a/src/multirtc/create_rtc.py b/src/multirtc/create_rtc.py index 2e768cf..dfc6982 100755 --- a/src/multirtc/create_rtc.py +++ b/src/multirtc/create_rtc.py @@ -518,11 +518,16 @@ def rtc(slc, geogrid, opts): 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 + out_geo_rtc_obj = None + if opts.apply_rtc: + 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 + if opts.apply_bistatic_delay or opts.apply_static_tropo: rg_lut, az_lut = compute_correction_lut( slc.source, @@ -611,26 +616,27 @@ def rtc(slc, geogrid, opts): out_geo_nlooks_obj.close_dataset() del out_geo_nlooks_obj - out_geo_rtc_obj.close_dataset() - del out_geo_rtc_obj + if opts.apply_rtc: + out_geo_rtc_gamma0_to_sigma0_obj.close_dataset() + del out_geo_rtc_gamma0_to_sigma0_obj - out_geo_rtc_gamma0_to_sigma0_obj.close_dataset() - del out_geo_rtc_gamma0_to_sigma0_obj + out_geo_rtc_obj.close_dataset() + del out_geo_rtc_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, - doppler=doppler, - ) + 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, + doppler=doppler, + ) t_end = time.time() logger.info(f'elapsed time: {t_end - t_start}') diff --git a/src/multirtc/geocode.py b/src/multirtc/geocode.py new file mode 100644 index 0000000..f12daab --- /dev/null +++ b/src/multirtc/geocode.py @@ -0,0 +1,19 @@ +"""Create a geocoded dataset for a multiple satellite platforms""" + +from pathlib import Path + +from multirtc.multirtc import SUPPORTED, run_multirtc + + +def create_parser(parser): + parser.add_argument('platform', choices=SUPPORTED, help='Platform to create geocoded dataset for') + parser.add_argument('granule', help='Data granule to create geocoded dataset for.') + parser.add_argument('--resolution', type=float, help='Resolution of the output dataset (m)') + parser.add_argument('--work-dir', type=Path, default=None, help='Working directory for processing') + return parser + + +def run(args): + if args.work_dir is None: + args.work_dir = Path.cwd() + run_multirtc(args.platform, args.granule, args.resolution, args.work_dir, apply_rtc=False) diff --git a/src/multirtc/multimetric/ale.py b/src/multirtc/multimetric/ale.py index 65b5552..78d9ad5 100644 --- a/src/multirtc/multimetric/ale.py +++ b/src/multirtc/multimetric/ale.py @@ -1,6 +1,5 @@ """Absolute Location Error (ALE) analysis""" -from argparse import ArgumentParser from datetime import datetime from pathlib import Path @@ -243,14 +242,3 @@ def run(args): assert args.filepath.exists(), f'File {args.filepath} does not exist.' ale(args.filepath, args.date, args.azmangle, args.project, basedir=args.basedir) - - -def main(): - parser = ArgumentParser(description=__doc__) - parser = create_parser(parser) - args = parser.parse_args() - run(args) - - -if __name__ == '__main__': - main() diff --git a/src/multirtc/multimetric/point_target.py b/src/multirtc/multimetric/point_target.py index a2235e1..2f1c579 100644 --- a/src/multirtc/multimetric/point_target.py +++ b/src/multirtc/multimetric/point_target.py @@ -1,6 +1,5 @@ """Point target (resolution, PSLR, and ISLR) analysis""" -from argparse import ArgumentParser from pathlib import Path import matplotlib.pyplot as plt @@ -111,14 +110,3 @@ def run(args): assert args.filepath.exists() args.basedir = Path(args.basedir).expanduser() analyze_point_targets(args.platform, args.filepath, args.project, basedir=args.basedir) - - -def main(): - parser = ArgumentParser(description='Point target (resolution, PSLR, and ISLR) analysis') - parser = create_parser(parser) - args = parser.parse_args() - run(args) - - -if __name__ == '__main__': - main() diff --git a/src/multirtc/multimetric/rle.py b/src/multirtc/multimetric/rle.py index 01633fc..998e8a6 100644 --- a/src/multirtc/multimetric/rle.py +++ b/src/multirtc/multimetric/rle.py @@ -1,6 +1,5 @@ """Relative Location Error (RLE) analysis""" -from argparse import ArgumentParser from dataclasses import dataclass from pathlib import Path from tempfile import NamedTemporaryFile @@ -196,14 +195,3 @@ def run(args): args.basedir = Path(args.basedir) assert args.basedir.exists(), f'Base directory {args.basedir} does not exist' rle(args.reference, args.secondary, project=args.project, basedir=args.basedir) - - -def main(): - parser = ArgumentParser(description=__doc__) - parser = create_parser(parser) - args = parser.parse_args() - run(args) - - -if __name__ == '__main__': - main() diff --git a/src/multirtc/multirtc.py b/src/multirtc/multirtc.py index b341745..a8bc71c 100644 --- a/src/multirtc/multirtc.py +++ b/src/multirtc/multirtc.py @@ -1,6 +1,5 @@ """Create an RTC dataset for a multiple satellite platforms""" -import argparse from pathlib import Path from burst2safe.burst2safe import burst2safe @@ -61,7 +60,7 @@ def get_slc(platform: str, granule: str, input_dir: Path) -> Slc: return slc -def run_multirtc(platform: str, granule: str, resolution: int, work_dir: Path) -> None: +def run_multirtc(platform: str, granule: str, resolution: int, work_dir: Path, apply_rtc=True) -> None: """Create an RTC or Geocoded dataset using the OPERA algorithm. Args: @@ -69,6 +68,7 @@ def run_multirtc(platform: str, granule: str, resolution: int, work_dir: Path) - granule: Granule name if data is available in ASF archive, or filename if granule is already downloaded. resolution: Resolution of the output RTC (in meters). work_dir: Working directory for processing. + apply_rtc: If True perform radiometric correction; if False, only geocode. """ input_dir, output_dir = prep_dirs(work_dir) slc = get_slc(platform, granule, input_dir) @@ -79,6 +79,7 @@ def run_multirtc(platform: str, granule: str, resolution: int, work_dir: Path) - opts = RtcOptions( dem_path=str(dem_path), output_dir=str(output_dir), + apply_rtc=apply_rtc, resolution=resolution, apply_bistatic_delay=slc.supports_bistatic_delay, apply_static_tropo=slc.supports_static_tropo, @@ -102,20 +103,4 @@ def create_parser(parser): def run(args): if args.work_dir is None: args.work_dir = Path.cwd() - run_multirtc(args.platform, args.granule, args.resolution, args.work_dir) - - -def main(): - """Create a RTC or geocoded dataset for a multiple satellite platforms - - Example command: - multirtc UMBRA umbra_image.ntif --resolution 40 - """ - parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.ArgumentDefaultsHelpFormatter) - parser = create_parser(parser) - args = parser.parse_args() - run(args) - - -if __name__ == '__main__': - main() + run_multirtc(args.platform, args.granule, args.resolution, args.work_dir, apply_rtc=True) diff --git a/src/multirtc/rtc_options.py b/src/multirtc/rtc_options.py index 28b7f86..d850603 100644 --- a/src/multirtc/rtc_options.py +++ b/src/multirtc/rtc_options.py @@ -10,11 +10,13 @@ class RtcOptions: output_dir: str dem_path: str - apply_rtc: bool = True + apply_rtc: bool + resolution: float apply_thermal_noise: bool = True apply_abs_rad: bool = True apply_bistatic_delay: bool = True apply_static_tropo: bool = True + terrain_radiometry: str = 'gamma0' # 'gamma0' or 'sigma0' apply_valid_samples_sub_swath_masking: bool = True apply_shadow_masking: bool = True dem_interpolation_method: str = 'biquintic' @@ -28,33 +30,22 @@ class RtcOptions: clip_min: float = np.nan clip_max: float = np.nan upsample_radar_grid: bool = False - terrain_radiometry: str = 'gamma0' # 'gamma0' or 'sigma0' rtc_algorithm_type: str = 'area_projection' # 'area_projection' or 'bilinear_distribution' input_terrain_radiometry: str = 'beta0' - rtc_min_value_db: int = -30.0 + rtc_min_value_db: float = -30.0 rtc_upsampling: int = 2 rtc_area_beta_mode: str = 'auto' geo2rdr_threshold: float = 1.0e-7 geo2rdr_numiter: int = 50 rdr2geo_threshold: float = 1.0e-7 rdr2geo_numiter: int = 25 - output_epsg: int = None - resolution: int = 30 + output_epsg: int | None = None def __post_init__(self): if not self.apply_rtc: - if self.save_rtc_anf: - raise ValueError('RTC ANF flags are only available with RTC enabled') - if self.save_rtc_anf_gamma0_to_sigma0: - raise ValueError('RTC ANF gamma0 to sigma0 flags are only available with RTC enabled') - - if self.terrain_radiometry == 'sigma0' and self.save_rtc_anf_gamma0_to_sigma0: - raise ValueError('RTC ANF gamma0 to sigma0 flags are only available with output type set to gamma0') + self.terrain_radiometry: str = 'sigma0' - if self.apply_rtc: - self.layer_name_rtc_anf = f'rtc_anf_{self.terrain_radiometry}_to_{self.input_terrain_radiometry}' - else: - self.layer_name_rtc_anf = '' + self.layer_name_rtc_anf = f'rtc_anf_{self.terrain_radiometry}_to_{self.input_terrain_radiometry}' if self.dem_interpolation_method == 'biquintic': self.dem_interpolation_method_isce3 = isce3.core.DataInterpMethod.BIQUINTIC @@ -94,11 +85,6 @@ def __post_init__(self): else: raise ValueError(f'Invalid terrain radiometry: {self.terrain_radiometry}') - if self.apply_rtc: - self.layer_name_rtc_anf = f'rtc_anf_{self.terrain_radiometry}_to_{self.input_terrain_radiometry}' - else: - self.layer_name_rtc_anf = '' - if self.rtc_area_beta_mode == 'pixel_area': self.rtc_area_beta_mode_isce3 = isce3.geometry.RtcAreaBetaMode.PIXEL_AREA elif self.rtc_area_beta_mode == 'projection_angle':