Merge pull request #23 from forrestfwilliams/develop

Release v0.3.1

Author: forrestfwilliams

CHANGELOG.md

@@ -6,6 +6,14 @@ 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.3.1]
+
+### Changed
+* Loading of SICD data during beta0/sigma0 creation to a chunked strategy to reduce memory requirements
+
+### Fixed
+* Geolocation issue for prototype Umbra workflow related to switching to local UTM zone during processing
+
 ## [0.3.0]
 
 ### Changed

README.md

@@ -2,9 +2,11 @@
 
 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](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.**
+> [!WARNING]
+> This package is still in early development. Users are encouraged to not use this package in production or other critical contexts until the v1.0.0 release.
 
-The RTC algorithm utilized by this library is described in [Shiroma et al., 2023](https://doi.org/10.1109/TGRS.2022.3147472).
+> [!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
 MultiRTC allows users to create RTC products from SLC data for multiple SAR sensor platforms. Currently this list includes:
@@ -23,14 +25,21 @@ multirtc PLATFORM SLC-GRANULE --resolution RESOLUTION --work-dir WORK-DIR
 ```
 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.
+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 products are in Sigma0 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 pixel values represent sigma0 power.
 
 ### 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.
 
+## When will support for [insert SAR provider here] products be added?
+We're currently working on this package on a "best effort" basis with no specific timeline for any particular dataset. We would love to add support for every SAR dataset ASAP, but we only have so much time to devote to this package. If you want a particular dataset to be prioritized there are several things you can do:
+
+- [Open an issue](https://github.com/forrestfwilliams/multirtc/issues/new) requesting support for your dataset and encourage others to like or comment on it.
+- Provides links to example datasets over the Rosamond, California corner reflector site (Lat/Lon 34.799,-118.095) for performing cal/val.
+- Reach out to us about funding the development required to add your dataset.
+
 ## 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).
 2. Download a local version of the `multirtc` repository (`git clone https://github.com/forrestfwilliams/multirtc.git`)

src/multirtc/create_rtc.py

@@ -10,6 +10,7 @@
 from scipy import ndimage
 from tqdm import tqdm
 
+from multirtc.define_geogrid import get_point_epsg
 from multirtc.prep_burst import S1BurstSlc
 from multirtc.s1burst_corrections import compute_correction_lut
 from multirtc.sicd import SicdSlc
@@ -360,13 +361,11 @@ def rtc(slc, geogrid, opts):
 
     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)
+        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)
+        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
@@ -522,8 +521,8 @@ def pfa_prototype_geocode(sicd, geogrid, dem_path, output_dir):
     slc_lut = isce3.core.LUT2d(
         np.arange(sigma0_data.shape[1]), np.arange(sigma0_data.shape[0]), sigma0_data, interp_method
     )
-    local2ecef = pyproj.Transformer.from_crs(f'EPSG:{geogrid.epsg}', 'EPSG:4978', always_xy=True)
-    local2ll = pyproj.Transformer.from_crs(f'EPSG:{geogrid.epsg}', 'EPSG:4326', always_xy=True)
+    assert geogrid.epsg == 4326, 'Only EPSG:4326 is supported for PFA prototype geocoding'
+    ll2ecef = pyproj.Transformer.from_crs('EPSG:4326', 'EPSG:4978', always_xy=True)
     dem_raster = isce3.io.Raster(str(dem_path))
     dem = isce3.geometry.DEMInterpolator()
     dem.load_dem(dem_raster)
@@ -535,9 +534,8 @@ def pfa_prototype_geocode(sicd, geogrid, dem_path, output_dir):
     for i, j in tqdm(itertools.product(range(geogrid.width), range(geogrid.length)), total=n_iters):
         x = geogrid.start_x + (i * geogrid.spacing_x)
         y = geogrid.start_y + (j * geogrid.spacing_y)
-        lon, lat = local2ll.transform(x, y)
         hae = dem.interpolate_lonlat(np.deg2rad(x), np.deg2rad(y))  # ISCE3 expects lat/lon to be in radians!
-        ecef_x, ecef_y, ecef_z = local2ecef.transform(x, y, hae)
+        ecef_x, ecef_y, ecef_z = ll2ecef.transform(x, y, hae)
         row, col = sicd.geo2rowcol(np.array([(ecef_x, ecef_y, ecef_z)]))[0]
         if slc_lut.contains(row, col):
             output[j, i] = slc_lut.eval(row, col)
@@ -556,3 +554,6 @@ def pfa_prototype_geocode(sicd, geogrid, dem_path, output_dir):
     out_ds.GetRasterBand(1).SetNoDataValue(np.nan)
     out_ds.SetMetadata({'AREA_OR_POINT': 'Area'})
     out_ds = None
+
+    local_epsg = get_point_epsg(geogrid.start_y, geogrid.start_x)
+    gdal.Warp(str(output_path), str(output_path), dstSRS=f'EPSG:{local_epsg}', format='GTiff')

src/multirtc/sicd.py

@@ -9,8 +9,8 @@
 from sarpy.io.complex.sicd import SICDReader
 from shapely.geometry import Point, Polygon
 
-from multirtc import define_geogrid
 from multirtc.base import SlcTemplate, to_isce_datetime
+from multirtc.define_geogrid import get_point_epsg
 
 
 def check_poly_order(poly):
@@ -19,8 +19,8 @@ def check_poly_order(poly):
 
 class SicdSlc:
     def __init__(self, sicd_path):
-        reader = SICDReader(str(sicd_path.expanduser().resolve()))
-        sicd = reader.get_sicds_as_tuple()[0]
+        self.reader = SICDReader(str(sicd_path.expanduser().resolve()))
+        sicd = self.reader.get_sicds_as_tuple()[0]
         self.source = sicd
         self.id = Path(sicd_path).with_suffix('').name
         self.filepath = Path(sicd_path)
@@ -50,7 +50,7 @@ def __init__(self, sicd_path):
         self.raw_time_coa_poly = sicd.Grid.TimeCOAPoly
         last_line_time = self.raw_time_coa_poly(0, self.shape[1] - self.shift[1])
         first_line_time = self.raw_time_coa_poly(0, -self.shift[1])
-        self.az_reversed = last_line_time > first_line_time
+        self.az_reversed = last_line_time < first_line_time
         self.arp_pos = sicd.SCPCOA.ARPPos.get_array()
         self.scp_pos = sicd.GeoData.SCP.ECF.get_array()
         azimuth_angle, elevation_angle = self.calculate_look_angles()
@@ -88,71 +88,69 @@ def calculate_look_angles(self):
         elevation = np.arcsin(up / np.linalg.norm(topocentric))
         return np.rad2deg(azimuth), np.rad2deg(elevation)
 
-    def get_xrow_ycol(self) -> np.ndarray:
+    def get_xrow_ycol(self, rowrange=None, colrange=None) -> 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]
+        rowlen = self.shape[0] if rowrange is None else rowrange[1] - rowrange[0]
+        collen = self.shape[1] if colrange is None else colrange[1] - colrange[0]
+        rowoffset = self.scp_index[0] if rowrange is None else self.scp_index[0] + rowrange[0]
+        coloffset = self.scp_index[1] if colrange is None else self.scp_index[1] + colrange[0]
+
+        irow = np.tile(np.arange(rowlen), (collen, 1)).T
+        irow -= rowoffset
         xrow = irow * self.spacing[0]
 
-        icol = np.tile(np.arange(self.shape[1]), (self.shape[0], 1))
-        icol -= self.scp_index[1]
+        icol = np.tile(np.arange(collen), (rowlen, 1))
+        icol -= coloffset
         ycol = icol * self.spacing[1]
         return xrow, ycol
 
-    def load_data(self):
-        reader = SICDReader(str(self.filepath))
-        data = reader[:, :]
-        return data
-
-    def load_scaled_data(self, scale, power=False):
+    def load_scaled_data(self, scale, power=False, rowrange=None, colrange=None):
         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()
+        xrow, ycol = self.get_xrow_ycol(rowrange=rowrange, colrange=colrange)
+        if colrange is not None and rowrange is not None:
+            data = self.reader[rowrange[0] : rowrange[1], colrange[0] : colrange[1]]
+        elif colrange is None and rowrange is None:
+            data = self.reader[:, :]
+        else:
+            raise ValueError('Both xrange and yrange must be provided or neither.')
+
         scale_factor = polyval2d(xrow, ycol, coeff)
-        data = self.load_data()
+        del xrow, ycol  # deleting for memory management
+
         if power:
-            scaled_data = (data.real**2 + data.imag**2) * scale_factor
+            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
+            data = data * np.sqrt(scale_factor)
+        return data
 
+    def create_complex_beta0(self, outpath, row_iter=256):
         driver = gdal.GetDriverByName('GTiff')
-        ds = driver.Create(str(outpath), scaled_data.shape[1], scaled_data.shape[0], 1, gdal.GDT_CFloat32)
+        # Shape transposed for ISCE3 expectations
+        ds = driver.Create(str(outpath), self.shape[0], self.shape[1], 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:
+        n_chunks = int(np.floor(self.shape[0] // row_iter)) + 1
+        for i in range(n_chunks):
+            start_row = i * row_iter
+            end_row = min((i + 1) * row_iter, self.shape[0])
+            rowrange = [start_row, end_row]
+            colrange = [0, self.shape[1]]
+            scaled_data = self.load_scaled_data('beta0', power=False, rowrange=rowrange, colrange=colrange)
+            # Shape transposed for ISCE3 expectations
             if self.az_reversed:
                 scaled_data = scaled_data[:, ::-1].T
             else:
                 scaled_data = scaled_data.T
+            # Offset transposed to match ISCE3 expectations
+            band.WriteArray(scaled_data, xoff=start_row, yoff=0)
 
-        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.FlushCache()
         ds = None
 
 
@@ -304,17 +302,24 @@ def geo2rowcol(self, xyz: np.ndarray) -> np.ndarray:
         return row_col
 
     def create_geogrid(self, spacing_meters):
-        epsg_local = define_geogrid.get_point_epsg(self.center.y, self.center.x)
         ecef = pyproj.CRS(4978)  # ECEF on WGS84 Ellipsoid
-        local = pyproj.CRS(epsg_local)
-        ecef2local = pyproj.Transformer.from_crs(ecef, local, always_xy=True)
-        x_spacing = spacing_meters
-        y_spacing = -1 * spacing_meters
+        lla = pyproj.CRS(4979)  # WGS84 lat/lon/ellipsoid height
+        local_utm = pyproj.CRS(get_point_epsg(self.center.y, self.center.x))
+        lla2utm = pyproj.Transformer.from_crs(lla, local_utm, always_xy=True)
+        utm2lla = pyproj.Transformer.from_crs(local_utm, lla, always_xy=True)
+        ecef2lla = pyproj.Transformer.from_crs(ecef, lla, always_xy=True)
+
+        lla_point = (self.center.x, self.center.y)
+        utm_point = lla2utm.transform(*lla_point)
+        utm_point_shift = (utm_point[0] + spacing_meters, utm_point[1])
+        lla_point_shift = utm2lla.transform(*utm_point_shift)
+        x_spacing = lla_point_shift[0] - lla_point[0]
+        y_spacing = -1 * x_spacing
 
         points = np.array([(0, 0), (0, self.shape[1]), self.shape, (self.shape[0], 0)])
-        geos = self.rowcol2geo(points, hae=self.scp_hae)
+        geos = self.rowcol2geo(points, self.scp_hae)
 
-        points = np.vstack(ecef2local.transform(geos[:, 0], geos[:, 1], geos[:, 2])).T
+        points = np.vstack(ecef2lla.transform(geos[:, 0], geos[:, 1], geos[:, 2])).T
         minx, maxx = np.min(points[:, 0]), np.max(points[:, 0])
         miny, maxy = np.min(points[:, 1]), np.max(points[:, 1])
 
@@ -327,6 +332,6 @@ def create_geogrid(self, spacing_meters):
             spacing_y=float(y_spacing),
             length=int(length),
             width=int(width),
-            epsg=epsg_local,
+            epsg=4326,
         )
         return geogrid