Release v0.5.1

.github/workflows/build-and-deploy-test.yml

@@ -8,7 +8,7 @@ jobs:
     name: Build and publish Python distributions TestPyPI
     runs-on: ubuntu-latest
     steps:
-    - uses: actions/checkout@v4
+    - uses: actions/checkout@v5
       with:
         fetch-depth: 0
     - name: Checkout lastest tagged version

.github/workflows/build-and-deploy.yml

@@ -8,7 +8,7 @@ jobs:
     name: Build and publish Python distributions to PyPI
     runs-on: ubuntu-latest
     steps:
-    - uses: actions/checkout@v4
+    - uses: actions/checkout@v5
       with:
         fetch-depth: 0
     - name: Checkout lastest tagged version

CHANGELOG.md

@@ -6,6 +6,12 @@ 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.1]
+
+### Added
+* Support for custom DEMs.
+* Support for prepping external DEMs to ISCE3-compatible format.
+
 ## [0.5.0]
 
 ### Added

README.md

@@ -71,7 +71,7 @@ PROJECT/
         |--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`).
+If you're encountering `permission denied` errors when running the container, make sure that the input and output folders are owned by the same group and user IDs that the container uses (`chown -R 1000:1000 ~/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:
@@ -89,7 +89,24 @@ More information on the metadata images can be found in the OPERA RTC Static Pro
 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. If the low resolution of the default DEM is causing radiometry issues, try using the `geocode` instead of `rtc` workflow.
+Currently, only the NISAR DEM can be automatically downloaded and it is the default DEM option. This is a global height above the WGS84 ellipsoid DEM sourced from the [COP-30 DEM](https://portal.opentopography.org/raster?opentopoID=OTSDEM.032021.4326.3).
+
+MultiRTC also supports custom DEMs via the use of the `--dem` option. For example:
+```bash
+multirtc rtc PLATFORM SLC-GRANULE --dem DEM-PATH --resolution RESOLUTION --work-dir WORK-DIR
+```
+We recommend using [OpenTopography](https://opentopography.org) to locate high-resolution DEMs for your area of interest. Note that custom DEM must meet three criteria to be used for processing:
+1. The DEM must fully cover the spatial extent of your input SLC granule.
+1. The DEM must be in the EPSG:4326 (lat/lon) spatial projection.
+1. The height values of the DEM must represent the height above the WGS84 ellipsoid.
+
+The `prepdem` subroutine can be used to reformat an existing DEM to meet the final two requirements:
+```bash
+multirtc prepdem INPUT-DEM-PATH OUTPUT-DEM-PATH VERTICAL-DATUM
+```
+Where `VERTICAL-DATUM` is the vertical datum of your input DEM. Currently the `WGS84`, `EGM2008`, and `NAD88` (over the United States) datums are supported.
+
+If the low resolution of the available DEM options is causing radiometry issues, try using the `geocode` workflow instead of `rtc`.
 
 ## Calibration & Validation Subcommands
 

scripts/nad83_to_wgs84.py

@@ -0,0 +1,72 @@
+from pathlib import Path
+
+import numpy as np
+from osgeo import gdal
+from pyproj import Transformer
+
+
+gdal.UseExceptions()
+
+
+def write_epsg4326_geotiff(output_path: Path, data: np.ndarray, transform: tuple) -> None:
+    """Write a numpy array to a GeoTIFF file.
+
+    Args:
+        output_path: Path to the output GeoTIFF file.
+        data: 2D numpy array containing the data to write.
+        transform: GeoTransform tuple for the raster.
+    """
+    driver = gdal.GetDriverByName('GTiff')
+    ds = driver.Create(str(output_path), data.shape[1], data.shape[0], 1, gdal.GDT_Float32)
+    ds.SetGeoTransform(transform)
+    ds.SetProjection('EPSG:4326')  # WGS84
+    ds.GetRasterBand(1).WriteArray(data)
+    ds.FlushCache()
+    ds = None
+
+
+def convert_from_nad83_to_wgs84(navd83: str) -> None:
+    """Convert NAD83(2011) ellipsoidal heights to WGS84 ellipsoidal heights.
+
+    Args:
+        navd83: Path to the NAD83(2011) geoid model.
+    """
+    ds = gdal.Open(navd83, gdal.GA_ReadOnly)
+    nad83_heights = ds.GetRasterBand(1).ReadAsArray()
+    transform = ds.GetGeoTransform()
+    lons = np.zeros_like(nad83_heights)
+    lats = np.zeros_like(nad83_heights)
+    wgs84_heights = np.zeros_like(nad83_heights)
+    ds = None
+
+    # NAD83(2011) ellipsoidal heights → WGS84 ellipsoidal heights
+    transformer = Transformer.from_crs('EPSG:6319', 'EPSG:4979', always_xy=True)
+    for i in range(nad83_heights.shape[0]):
+        for j in range(nad83_heights.shape[1]):
+            lon, lat = transform[0] + j * transform[1], transform[3] + i * transform[5]
+            lon_wgs84, lat_wgs84, h_wgs84 = transformer.transform(lon, lat, nad83_heights[i, j])
+            lons[i, j] = lon_wgs84
+            lats[i, j] = lat_wgs84
+            wgs84_heights[i, j] = h_wgs84
+
+    wgs84_lon_start = lons[0, 0]
+    wgs84_lon_step = np.diff(lons, axis=1).mean()
+    wgs84_lat_start = lats[0, 0]
+    wgs84_lat_step = np.diff(lats, axis=0).mean()
+    trans = [wgs84_lon_start, wgs84_lon_step, 0, wgs84_lat_start, 0, wgs84_lat_step]
+    if trans[0] > 180:
+        trans[0] -= 360
+    if (trans[0] + (wgs84_lon_step * wgs84_heights.shape[1])) > 180:
+        write_epsg4326_geotiff(f'{Path(navd83).stem}_wgs84_east.tif', wgs84_heights, trans)
+        west_trans = [trans[0] - 360, trans[1], trans[2], trans[3], trans[4], trans[5]]
+        write_epsg4326_geotiff(f'{Path(navd83).stem}_wgs84_west.tif', wgs84_heights, west_trans)
+    else:
+        write_epsg4326_geotiff(f'{Path(navd83).stem}_wgs84.tif', wgs84_heights, trans)
+
+
+convert_from_nad83_to_wgs84(
+    '/vsicurl/https://github.com/OSGeo/PROJ-data/raw/refs/heads/master/us_noaa/us_noaa_g2012ba0.tif'
+)
+convert_from_nad83_to_wgs84(
+    '/vsicurl/https://github.com/OSGeo/PROJ-data/raw/refs/heads/master/us_noaa/us_noaa_g2012bu0.tif'
+)

src/multirtc/__main__.py

@@ -1,6 +1,6 @@
 import argparse
 
-from multirtc import geocode, multirtc
+from multirtc import dem, geocode, multirtc
 from multirtc.multimetric import ale, point_target, rle
 
 
@@ -18,6 +18,9 @@ def main():
     geocode_parser = multirtc.create_parser(subparsers.add_parser('geocode', help=geocode.__doc__))
     geocode_parser.set_defaults(func=geocode.run)
 
+    dem_parser = dem.create_parser(subparsers.add_parser('prepdem', help=dem.__doc__))
+    dem_parser.set_defaults(func=dem.run)
+
     ale_parser = ale.create_parser(subparsers.add_parser('ale', help=ale.__doc__))
     ale_parser.set_defaults(func=ale.run)
 

src/multirtc/dem.py

@@ -1,16 +1,74 @@
+"""Prepare an external DEM for use in MultiRTC"""
+
 from concurrent.futures import ThreadPoolExecutor
 from itertools import product
 from pathlib import Path
+from shutil import copyfile
+from tempfile import NamedTemporaryFile
 
 import numpy as np
 import shapely
 from hyp3lib.fetch import download_file
-from osgeo import gdal
+from osgeo import gdal, osr
 from shapely.geometry import LinearRing, Polygon, box
 
 
 gdal.UseExceptions()
 URL = 'https://nisar.asf.earthdatacloud.nasa.gov/STATIC/DEM/v1.1/EPSG4326'
+EGM2008_GEOID = {
+    'WORLD': [
+        '/vsicurl/https://asf-dem-west.s3.amazonaws.com/GEOID/us_nga_egm2008_1.tif',
+        box(-180.0083333, -90.0083333, 180.0083333, 90.0083333),
+    ]
+}
+NAD88 = {
+    'CONUS': [
+        '/vsicurl/https://asf-dem-west.s3.amazonaws.com/GEOID/us_noaa_g2012bu0_wgs84.tif',
+        box(-130.0083313, 23.9916667, -59.9920366, 58.0083351),
+    ],
+    'AK_EAST': [
+        '/vsicurl/https://asf-dem-west.s3.amazonaws.com/GEOID/us_noaa_g2012ba0_wgs84_east.tif',
+        box(171.9916687, 48.9916583, 234.008, 72.008),
+    ],
+    'AK_WEST': [
+        '/vsicurl/https://asf-dem-west.s3.amazonaws.com/GEOID/us_noaa_g2012ba0_wgs84_west.tif',
+        box(-188.008, 48.992, -125.9915921, 72.0083313),
+    ],
+}
+VALID_DATUMS = ['WGS84', 'EGM2008', 'NAD88']
+
+
+def get_bbox_from_info(info: dict) -> Polygon:
+    minx = info['cornerCoordinates']['lowerLeft'][0]
+    miny = info['cornerCoordinates']['lowerLeft'][1]
+    maxx = info['cornerCoordinates']['upperRight'][0]
+    maxy = info['cornerCoordinates']['upperRight'][1]
+    return box(minx, miny, maxx, maxy)
+
+
+def validate_dem(dem_path: Path, footprint: Polygon) -> None:
+    """Validate that the DEM file is in EPSG:4326 and contains the given footprint.
+
+    Args:
+        dem_path: Path to the DEM file.
+        footprint: Polygon representing the area of interest.
+
+    Raises:
+        ValueError: If the DEM does not cover the footprint.
+    """
+    info = gdal.Info(str(dem_path), format='json')
+
+    srs = osr.SpatialReference()
+    srs.ImportFromWkt(info['coordinateSystem']['wkt'])
+    assert int(srs.GetAttrValue('AUTHORITY', 1)) == 4326, f'DEM file {dem_path} is not in EPSG:4326 projection.'
+
+    dem_extent = get_bbox_from_info(info)
+    if not dem_extent.contains(footprint):
+        dem_bound_str = ', '.join([str(round(x, 3)) for x in dem_extent.bounds])
+        footprint_bound_str = ', '.join([str(round(x, 3)) for x in footprint.bounds])
+        raise ValueError(
+            f'DEM does not fully cover the footprint: ({dem_bound_str}) for DEM, vs ({footprint_bound_str})'
+        )
 
 
 def check_antimeridean(poly: Polygon) -> list[Polygon]:
@@ -126,3 +184,112 @@ def download_opera_dem_for_footprint(output_path: Path, footprint: Polygon, buff
 
     ds = None
     [Path(f).unlink() for f in input_files + [vrt_filepath]]
+
+
+def get_correction_geoid(bbox: Polygon, input_datum: str) -> str:
+    """Get the path to the geoid correction file based on the bounding box and input datum."""
+    if input_datum.upper() == 'EGM2008':
+        return EGM2008_GEOID['WORLD'][0]
+
+    if input_datum.upper() == 'NAD88':
+        for region, (correction_path, region_bbox) in NAD88.items():
+            if bbox.intersects(region_bbox):
+                return correction_path
+
+    raise ValueError(f'No suitable geoid correction found for datum {input_datum} and bounding box {bbox}.')
+
+
+def convert_to_height_above_ellipsoid(dem_file: Path, input_datum) -> None:
+    assert input_datum in VALID_DATUMS, f'Input datum must be one of {VALID_DATUMS}, got {input_datum}.'
+    if input_datum.upper() == 'WGS84':
+        return
+    dem_info = gdal.Info(str(dem_file), format='json')
+    dem_bbox = get_bbox_from_info(dem_info)
+    correction_path = get_correction_geoid(dem_bbox, input_datum)
+    with NamedTemporaryFile() as geoid_file:
+        gdal.Warp(
+            geoid_file.name,
+            correction_path,
+            dstSRS=dem_info['coordinateSystem']['wkt'],
+            outputBounds=dem_bbox.bounds,
+            width=dem_info['size'][0],
+            height=dem_info['size'][1],
+            resampleAlg='cubic',
+            multithread=True,
+            format='GTiff',
+        )
+        geoid_ds = gdal.Open(geoid_file.name)
+        geoid_data = geoid_ds.GetRasterBand(1).ReadAsArray()
+        del geoid_ds
+
+        dem_ds = gdal.Open(str(dem_file), gdal.GA_Update)
+        dem_data = dem_ds.GetRasterBand(1).ReadAsArray()
+        nan_value = dem_ds.GetRasterBand(1).GetNoDataValue()
+        if nan_value is not None:
+            nan_mask = dem_data == nan_value
+            dem_data += geoid_data
+            dem_data[nan_mask] = nan_value
+        else:
+            dem_data += geoid_data
+        dem_ds.GetRasterBand(1).WriteArray(dem_data)
+        dem_ds.FlushCache()
+        del dem_ds
+
+
+def set_nodata_value(dem_file: Path, nodata_value: float) -> None:
+    """Set the NoData value for the DEM file.
+
+    Args:
+        dem_file: Path to the DEM file.
+        nodata_value: Value to set as NoData.
+    """
+    dem_ds = gdal.Open(str(dem_file), gdal.GA_Update)
+    band = dem_ds.GetRasterBand(1)
+    if band.GetNoDataValue() is not None:
+        data = band.ReadAsArray()
+        data[data == band.GetNoDataValue()] = nodata_value
+        band.WriteArray(data)
+    band.SetNoDataValue(nodata_value)
+    band.FlushCache()
+    del dem_ds
+
+
+def prep_dem(input_path: Path, output_path: Path, input_datum: str) -> None:
+    """Prepare the DEM for processing by reprojecting it to EPSG:4326 and (optionally) converting it to height above ellipsoid.
+
+    Args:
+        dem_path: Path to the DEM file.
+        input_datum: Datum of the input DEM, either 'WGS84', 'EGM2008', 'NAD88'.
+        output_path: Path where the prepared DEM will be saved.
+    """
+    assert input_datum in VALID_DATUMS, f'Input datum must be one of {VALID_DATUMS}, got {input_datum}.'
+    copyfile(input_path, output_path)
+    info = gdal.Info(str(input_path), format='json')
+    srs = osr.SpatialReference()
+    srs.ImportFromWkt(info['coordinateSystem']['wkt'])
+    if int(srs.GetAttrValue('AUTHORITY', 1)) != 4326:
+        gdal.Warp(
+            str(output_path),
+            str(output_path),
+            dstSRS='EPSG:4326',
+            resampleAlg='cubic',
+            multithread=True,
+        )
+    convert_to_height_above_ellipsoid(output_path, input_datum.upper())
+    set_nodata_value(output_path, 0)
+
+
+def create_parser(parser):
+    parser.add_argument('input_dem', help='Path to the input DEM file to prepare for processing.')
+    parser.add_argument('output_dem', help='Path where the prepared DEM will be saved.')
+    parser.add_argument(
+        'input_datum', choices=VALID_DATUMS, help='Datum of the input DEM (supported: WGS84, EGM2008, NAD88).'
+    )
+    return parser
+
+
+def run(args):
+    args.input_dem = Path(args.input_dem)
+    assert args.input_dem.exists(), f'DEM file {args.input_dem} does not exist.'
+    args.output_dem = Path(args.output_dem)
+    prep_dem(args.input_dem, args.output_dem, args.input_datum)

src/multirtc/geocode.py

@@ -9,11 +9,14 @@ 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('--dem', type=Path, default=None, help='Path to the DEM to use for processing')
     parser.add_argument('--work-dir', type=Path, default=None, help='Working directory for processing')
     return parser
 
 
 def run(args):
+    if args.dem is not None:
+        assert args.dem.exists(), f'DEM file {args.dem} does not exist.'
     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)
+    run_multirtc(args.platform, args.granule, args.resolution, args.work_dir, args.dem, apply_rtc=False)