Release v0.5.0

CHANGELOG.md

@@ -6,6 +6,22 @@ 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.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.
+
 ## [0.4.0]
 
 ### Changed

README.md

@@ -8,18 +8,33 @@ A python library for creating ISCE3-based RTCs for multiple SAR data sources
 > [!IMPORTANT]
 > 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 package merely allows others to use their algorithm with a wider set of SAR data sources. The RTC algorithm utilized by this package is described in [Shiroma et al., 2023](https://doi.org/10.1109/TGRS.2022.3147472).
 
-## Usage
+## Dataset Support
 MultiRTC allows users to create RTC products from SLC data for multiple SAR sensor platforms, and provides utilities for assessing the resulting products. All utilities can be accessed via CLI pattern `multirtc SUBCOMMAND ARGS`, with the primary subcommand `multirtc rtc`.
 
-Currently the list of supported datasets includes:
-
-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)
-- [ICEYE SICD SLCs](https://sar.iceye.com/5.0/productFormats/slc/)
+Below is a list of relevant SAR data sources and their support status:
+
+| Mission    | File Format | Image Mode        | Image Grid Type | Status      |
+|------------|-------------|-------------------|-----------------|-------------|
+| Sentinel-1 | SAFE        | Burst IW          | Range Doppler   | Supported   |
+| Sentinel-1 | SAFE        | Full-frame IW     | Range Doppler   | Unsupported |
+| Sentinel-1 | SAFE        | Burst EW          | Range Doppler   | Unsupported |
+| Sentinel-1 | SAFE        | Full-frame EW     | Range Doppler   | Unsupported |
+| Capella    | SICD        | Spotlight         | Polar           | Supported\* |
+| Capella    | SICD        | Sliding Spotlight | Range Doppler   | Supported   |
+| Capella    | SICD        | Stripmap          | Range Doppler   | Supported   | 
+| Iceye      | SICD        | Dwell             | Range Doppler   | Supported   | 
+| Iceye      | SICD        | Spotlight         | Range Doppler   | Supported   | 
+| Iceye      | SICD        | Sliding Spotlight | Range Doppler   | Supported   |
+| Iceye      | SICD        | Stripmap          | Range Doppler   | Supported   |
+| Iceye      | SICD        | Scan              | Range Doppler   | Supported   |
+| Umbra      | SICD        | Dwell             | Polar           | Supported\* |
+| Umbra      | SICD        | Spotlight         | Polar           | Supported\* |
+
+I have done my best to accurately reflect the support status of each SAR image type, but please let me know if I have made any mistakes. Note that some commercial datasets used to use polar instead of range doppler image grids for specific images modes. This table is based on the image grid types currently being used.
+
+\*Polar image grid support is implemented via the [approach detailed by Piyush Agram](https://arxiv.org/abs/2503.07889v1) in his recent technical note. I have implemented his method in a fork of the main ISCE3 repo, which you can view [here](https://github.com/forrestfwilliams/isce3/tree/pfa). The long-term plan is to merge this into the main ISCE3 repo but until that is complete, polar grid support is only available via this project's `pfa`-suffixed docker containers. See the running via docker section for more details.
 
-Geocode Only:
-- [UMBRA SICD SLCs](https://help.umbra.space/product-guide/umbra-products/umbra-product-specifications)
+## Usage
 
 To create an RTC, use the `multirtc` CLI entrypoint using the following pattern:
 
@@ -30,13 +45,57 @@ Where `PLATFORM` is the name of the satellite platform (currently `S1`, `CAPELLA
 
 Output RTC pixel values represent gamma0 power.
 
-### Current Umbra Implementation
-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 pixel values represent sigma0 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:
+
+```bash
+docker run -it --rm \
+    -e EARTHDATA_USERNAME=YOUR_USERNAME_HERE \
+    -e EARTHDATA_PASSWORD=YOUR_PASSWORD_HERE \
+    -v ~/LOCAL_PATH/PROJECT:/home/conda/PROJECT \
+    ghcr.io/forrestfwilliams/multirtc:VERSION.pfa \
+    rtc PLATFORM SLC-GRANULE --resolution RESOLUTION --work-dir PROJECT
+```
+The local `project1` directory can be a name of your choosing and should have the structure:
+```
+PROJECT/
+    |--input/
+        |--input.slc (if needed)
+    |--output/
+```
+If you're encountering `permission denied` errors when running the container, make sure other users are allowed to read/write to your project directory (`chmod -R a+rwX ~/LOCAL_PATH/PROJECT`).
+
+### Output Layers
+MultiRTC outputs one main RTC image and seven metadata images as GeoTIFFs. All layers follow the naming schema `{FILEID}_{DATASET}.tif`, with the main RTC image omiting the `_{DATASET}` component. The layers are as follows:
+- `FILEID.tif`: The radiometric and terrain corrected backscatter data in gamma0 radiometry.
+- `FILEID_incidence_angle.tif`: The angle between the LOS vector and the ellipsoid normal at the target.
+- `FILEID_interpolated_dem.tif`: The DEM used of calculating layover/shadow.
+- `FILEID_local_incidence_angle.tif`: The angle between the LOS vector and terrain normal vector at the target.
+- `FILEID_mask.tif`: The layover/shadow mask. `0` is no shadow or shadow, `1` is shadow, `2` is layover and `3` is layover and shadow.
+- `FILEID_number_of_looks.tif`: The number of radar samples used to compute each output image pixel.
+- `FILEID_rtc_anf_gamma0_to_beta0.tif`: The conversion values needed to normalize the gamma0 backscatter to beta0.
+- `FILEID_rtc_anf_gamma0_to_sigma0.tif`: The conversion values needed to normalize the gamma0 backscatter to sigma0.
+
+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
+
+> [!WARNING]
+> This submodule currently only support Umbra SICD data! Reach out if you would like to see this submodule expanded to other datasets.
+
 MultiRTC includes three calibration and validation (cal/val) subcommands for assessing the geometric and radiometric quality of SAR products. These tools are useful for analyzing geolocation, co-registration, and impulse response performance.
 
 ### `ale` Absolute Location Error

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)

src/multirtc/create_rtc.py

@@ -2,6 +2,7 @@
 import logging
 import os
 import time
+from pathlib import Path
 
 import isce3
 import numpy as np
@@ -143,7 +144,8 @@ def compute_correction_lut(
     rg_lut = isce3.core.LUT2d(
         bistatic_delay.x_start, bistatic_delay.y_start, bistatic_delay.x_spacing, bistatic_delay.y_spacing, tropo
     )
-
+    [x.unlink() for x in Path(scratch_path).glob('*.hdr') if x.is_file()]
+    [x.unlink() for x in Path(scratch_path).glob('*.rdr') if x.is_file()]
     return rg_lut, az_lut
 
 
@@ -388,10 +390,10 @@ def save_intermediate_geocode_files(
     names = [
         LAYER_NAME_LOCAL_INCIDENCE_ANGLE,
         LAYER_NAME_INCIDENCE_ANGLE,
-        LAYER_NAME_PROJECTION_ANGLE,
-        LAYER_NAME_RTC_ANF_PROJECTION_ANGLE,
-        # LAYER_NAME_RANGE_SLOPE, # FIXME
         LAYER_NAME_DEM,
+        # LAYER_NAME_PROJECTION_ANGLE,
+        # LAYER_NAME_RTC_ANF_PROJECTION_ANGLE,
+        # LAYER_NAME_RANGE_SLOPE, # FIXME
     ]
     raster_objs = []
     for name in names:
@@ -408,10 +410,10 @@ def save_intermediate_geocode_files(
     (
         local_incidence_angle_raster,
         incidence_angle_raster,
-        projection_angle_raster,
-        rtc_anf_projection_angle_raster,
-        # range_slope_raster, # FIXME
         interpolated_dem_raster,
+        # projection_angle_raster,
+        # rtc_anf_projection_angle_raster,
+        # range_slope_raster, # FIXME
     ) = raster_objs
 
     # TODO review this (Doppler)!!!
@@ -432,9 +434,9 @@ def save_intermediate_geocode_files(
         dem_interp_method_enum,
         incidence_angle_raster=incidence_angle_raster,
         local_incidence_angle_raster=local_incidence_angle_raster,
-        projection_angle_raster=projection_angle_raster,
-        simulated_radar_brightness_raster=rtc_anf_projection_angle_raster,
         interpolated_dem_raster=interpolated_dem_raster,
+        # projection_angle_raster=projection_angle_raster,
+        # simulated_radar_brightness_raster=rtc_anf_projection_angle_raster,
         # range_slope_angle_raster=range_slope_raster, # FIXME
     )
     for obj in output_obj_list:
@@ -516,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,
@@ -609,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}')
 

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)

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
 
@@ -17,6 +16,10 @@
 gdal.UseExceptions()
 
 
+def db(data):
+    return 10 * np.log10(data)
+
+
 def gaussfit(x, y, A, x0, y0, sigma_x, sigma_y, theta):
     theta = np.radians(theta)
     sigx2 = sigma_x**2
@@ -52,7 +55,10 @@ def plot_crs_on_image(cr_df, data, project, outdir):
     max_y = cr_df['yloc'].max() + buffer
 
     fig, ax = plt.subplots(figsize=(15, 7))
-    ax.imshow(data, cmap='gray', interpolation='bilinear', vmin=0.3, vmax=1.7, origin='upper')
+    data_db = db(data)
+    vmin = np.nanpercentile(data_db, 2)
+    vmax = np.nanpercentile(data_db, 98)
+    ax.imshow(data_db, cmap='gray', interpolation='bilinear', vmin=vmin, vmax=vmax, origin='upper')
     ax.set_xlim(min_x, max_x)
     ax.set_ylim(min_y, max_y)
     ax.axis('off')
@@ -130,7 +136,7 @@ def calculate_ale_for_cr(point, data, project, outdir, search_window=100, oversa
 
     plt.rcParams.update({'font.size': 14})
     fig, ax = plt.subplots(1, 3, figsize=(15, 7))
-    ax[0].imshow(centered_data, cmap='gray', interpolation=None, origin='upper')
+    ax[0].imshow(db(centered_data), cmap='gray', interpolation=None, origin='upper')
     ax[0].plot(xpeak_centered, ypeak_centered, 'r+', label='Return Peak')
     ax[0].plot(xreal_centered, yreal_centered, 'b+', label='CR Location')
     ax[0].legend()
@@ -144,7 +150,7 @@ def calculate_ale_for_cr(point, data, project, outdir, search_window=100, oversa
     ax[2].set_title(f'Gaussian Fit Corner Reflector (ID {int(point["ID"])})')
     [axi.axis('off') for axi in ax]
     fig.tight_layout()
-    fig.savefig(outdir / f'{project}_CR_{point["ID"]}.png', dpi=300, bbox_inches='tight')
+    fig.savefig(outdir / f'{project}_CR_{int(point["ID"])}.png', dpi=300, bbox_inches='tight')
 
     return point
 
@@ -236,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()

src/multirtc/multimetric/corner_reflector.py

@@ -83,7 +83,7 @@ def get_cr_df(bounds, date, azmangle, outdir):
 def add_geo_image_location(cr_df, geotransform, shape, epsg):
     bounds = get_bounds(geotransform, shape)
     x_start = bounds.bounds[0]
-    y_start = bounds.bounds[1]
+    y_start = bounds.bounds[3]
     x_spacing = geotransform[1]
     y_spacing = geotransform[5]
     blank = [np.nan] * cr_df.shape[0]
@@ -124,6 +124,10 @@ def add_rdr_image_location(slc, cr_df, search_radius):
         row_guess, col_guess = int(round(row_guess)), int(round(col_guess))
         row_range = (row_guess - search_radius, row_guess + search_radius)
         col_range = (col_guess - search_radius, col_guess + search_radius)
+        if row_range[0] < 0 or row_range[1] >= slc.shape[0] or col_range[0] < 0 or col_range[1] >= slc.shape[1]:
+            no_peak.append(idx)
+            print(f'CR {int(row["ID"])} outside of SLC bounds, skipping.')
+            continue
         data = slc.load_data(row_range, col_range)
         in_db = 10 * np.log10(np.abs(data))
         row_peak, col_peak = np.unravel_index(np.argmax(in_db, axis=None), in_db.shape)