From 29f0facf5dea2ad04fa95d67ba2b1df99776121b Mon Sep 17 00:00:00 2001
From: Forrest Williams <ffwilliams2@alaska.edu>
Date: Tue, 1 Jul 2025 13:46:29 -0500
Subject: [PATCH 01/11] starting work on custom dems

---
 scripts/nad83_to_wgs84.py |  72 +++++++++++++++++++++++
 src/multirtc/dem.py       | 121 +++++++++++++++++++++++++++++++++++++-
 src/multirtc/multirtc.py  |  14 +++--
 3 files changed, 202 insertions(+), 5 deletions(-)
 create mode 100644 scripts/nad83_to_wgs84.py

diff --git a/scripts/nad83_to_wgs84.py b/scripts/nad83_to_wgs84.py
new file mode 100644
index 0000000..776ab22
--- /dev/null
+++ b/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'
+)
diff --git a/src/multirtc/dem.py b/src/multirtc/dem.py
index b444ffa..a1daa21 100644
--- a/src/multirtc/dem.py
+++ b/src/multirtc/dem.py
@@ -1,16 +1,68 @@
 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.tif',
+        box(-130.0083313, 23.9916667, -59.9920366, 58.0083351),
+    ],
+    'ak_east': [
+        '/vsicurl/https://asf-dem-west.s3.amazonaws.com/GEOID/us_noaa_g2012ba0_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_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):
+        raise ValueError(f'DEM does not fully cover the footprint: {footprint}')
 
 
 def check_antimeridean(poly: Polygon) -> list[Polygon]:
@@ -126,3 +178,70 @@ 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.lower() == 'egm2008':
+        return EGM2008_GEOID['world'][0]
+
+    if input_datum.lower() == '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.lower() == '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()
+        dem_data += geoid_data
+        dem_ds.GetRasterBand(1).WriteArray(dem_data)
+        dem_ds.FlushCache()
+        del dem_ds
+
+
+def prep_dem(input_path: Path, input_datum: str, output_path: Path) -> 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. If None, the NISAR DEM will be downloaded.
+        input_datum: Datum of the input DEM, either 'WGS84', 'EGM2008', 'NAD88'.
+        output_path: Path where the prepared DEM will be saved.
+    """
+    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.lower())
diff --git a/src/multirtc/multirtc.py b/src/multirtc/multirtc.py
index 9677e7a..bbe338b 100644
--- a/src/multirtc/multirtc.py
+++ b/src/multirtc/multirtc.py
@@ -61,19 +61,22 @@ 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, dem_path: Path | None, work_dir: Path) -> None:
     """Create an RTC or Geocoded dataset using the OPERA algorithm.
 
     Args:
         platform: Platform type (e.g., 'UMBRA').
         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).
+        dem_path: Path to the DEM to use for processing. If None, the NISAR DEM will be downloaded.
         work_dir: Working directory for processing.
     """
     input_dir, output_dir = prep_dirs(work_dir)
     slc = get_slc(platform, granule, input_dir)
-    dem_path = input_dir / 'dem.tif'
-    dem.download_opera_dem_for_footprint(dem_path, slc.footprint)
+    if dem_path is None:
+        dem_path = input_dir / 'dem.tif'
+        dem.download_opera_dem_for_footprint(dem_path, slc.footprint)
+    dem.validate_dem(dem_path, slc.footprint)
     geogrid = slc.create_geogrid(spacing_meters=resolution)
     if slc.supports_rtc:
         opts = RtcOptions(
@@ -92,14 +95,17 @@ def create_parser(parser):
     parser.add_argument('platform', choices=SUPPORTED, help='Platform to create RTC for')
     parser.add_argument('granule', help='Data granule to create an RTC for.')
     parser.add_argument('--resolution', type=float, help='Resolution of the output RTC (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)
+    run_multirtc(args.platform, args.granule, args.resolution, args.dem, args.work_dir)
 
 
 def main():

From 98fc60c368212f0e50ff14f2531c05c2f389769e Mon Sep 17 00:00:00 2001
From: Forrest Williams <ffwilliams2@alaska.edu>
Date: Mon, 11 Aug 2025 13:30:21 -0500
Subject: [PATCH 02/11] fix docker text

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 5ed8c16..a3b415a 100644
--- a/README.md
+++ b/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:

From 582d771b5114eeaee7210e4b34944f721c500799 Mon Sep 17 00:00:00 2001
From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com>
Date: Tue, 12 Aug 2025 03:48:31 +0000
Subject: [PATCH 03/11] Bump actions/checkout from 4 to 5

Bumps [actions/checkout](https://github.com/actions/checkout) from 4 to 5.
- [Release notes](https://github.com/actions/checkout/releases)
- [Changelog](https://github.com/actions/checkout/blob/main/CHANGELOG.md)
- [Commits](https://github.com/actions/checkout/compare/v4...v5)

---
updated-dependencies:
- dependency-name: actions/checkout
  dependency-version: '5'
  dependency-type: direct:production
  update-type: version-update:semver-major
...

Signed-off-by: dependabot[bot] <support@github.com>
---
 .github/workflows/build-and-deploy-test.yml | 2 +-
 .github/workflows/build-and-deploy.yml      | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/build-and-deploy-test.yml b/.github/workflows/build-and-deploy-test.yml
index a4277f1..4954161 100644
--- a/.github/workflows/build-and-deploy-test.yml
+++ b/.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
diff --git a/.github/workflows/build-and-deploy.yml b/.github/workflows/build-and-deploy.yml
index d4818be..6b2ceef 100644
--- a/.github/workflows/build-and-deploy.yml
+++ b/.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

From 987f280abef7a7691b7870b656d0d635a0b846b1 Mon Sep 17 00:00:00 2001
From: Forrest Williams <ffwilliams2@alaska.edu>
Date: Mon, 18 Aug 2025 13:19:30 -0500
Subject: [PATCH 04/11] fix ruff

---
 src/multirtc/multirtc.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/multirtc/multirtc.py b/src/multirtc/multirtc.py
index 0479031..c8bb43d 100644
--- a/src/multirtc/multirtc.py
+++ b/src/multirtc/multirtc.py
@@ -60,7 +60,9 @@ 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, dem_path: Path | None = None, apply_rtc=True) -> None:
+def run_multirtc(
+    platform: str, granule: str, resolution: int, work_dir: Path, dem_path: Path | None = None, apply_rtc=True
+) -> None:
     """Create an RTC or Geocoded dataset using the OPERA algorithm.
 
     Args:

From 8b9e5b74621cde749c41e32cde9833a43a4a260c Mon Sep 17 00:00:00 2001
From: Forrest Williams <ffwilliams2@alaska.edu>
Date: Mon, 18 Aug 2025 13:20:49 -0500
Subject: [PATCH 05/11] update changelog

---
 CHANGELOG.md | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8ae14a5..a81c3bf 100644
--- a/CHANGELOG.md
+++ b/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

From 2f852fc7f42fa3240af494e6271648c591c83120 Mon Sep 17 00:00:00 2001
From: Forrest Williams <ffwilliams2@alaska.edu>
Date: Mon, 18 Aug 2025 14:19:06 -0500
Subject: [PATCH 06/11] add dem entrypoint

---
 src/multirtc/__main__.py |  5 ++++-
 src/multirtc/dem.py      | 22 ++++++++++++++++++++--
 2 files changed, 24 insertions(+), 3 deletions(-)

diff --git a/src/multirtc/__main__.py b/src/multirtc/__main__.py
index efbbdd4..d0273fb 100644
--- a/src/multirtc/__main__.py
+++ b/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)
 
diff --git a/src/multirtc/dem.py b/src/multirtc/dem.py
index a1daa21..89659e4 100644
--- a/src/multirtc/dem.py
+++ b/src/multirtc/dem.py
@@ -1,3 +1,4 @@
+"""Prepare an external DEM for use in MultiRTC"""
 from concurrent.futures import ThreadPoolExecutor
 from itertools import product
 from pathlib import Path
@@ -224,14 +225,15 @@ def convert_to_height_above_ellipsoid(dem_file: Path, input_datum) -> None:
         del dem_ds
 
 
-def prep_dem(input_path: Path, input_datum: str, output_path: Path) -> None:
+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. If None, the NISAR DEM will be downloaded.
+        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()
@@ -245,3 +247,19 @@ def prep_dem(input_path: Path, input_datum: str, output_path: Path) -> None:
             multithread=True,
         )
     convert_to_height_above_ellipsoid(output_path, input_datum.lower())
+
+
+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)

From 112ae5073f46357fbc99cf3b1931be9d3e1efb66 Mon Sep 17 00:00:00 2001
From: Forrest Williams <ffwilliams2@alaska.edu>
Date: Tue, 19 Aug 2025 10:46:41 -0500
Subject: [PATCH 07/11] fix paths and names

---
 src/multirtc/dem.py | 30 +++++++++++++++---------------
 1 file changed, 15 insertions(+), 15 deletions(-)

diff --git a/src/multirtc/dem.py b/src/multirtc/dem.py
index 89659e4..aad7d14 100644
--- a/src/multirtc/dem.py
+++ b/src/multirtc/dem.py
@@ -15,22 +15,22 @@
 gdal.UseExceptions()
 URL = 'https://nisar.asf.earthdatacloud.nasa.gov/STATIC/DEM/v1.1/EPSG4326'
 EGM2008_GEOID = {
-    'world': [
+    '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.tif',
+    '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_east.tif',
+    '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_west.tif',
+    '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),
     ],
 }
@@ -183,10 +183,10 @@ def download_opera_dem_for_footprint(output_path: Path, footprint: Polygon, buff
 
 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.lower() == 'egm2008':
-        return EGM2008_GEOID['world'][0]
+    if input_datum.upper() == 'EGM2008':
+        return EGM2008_GEOID['WORLD'][0]
 
-    if input_datum.lower() == 'nad88':
+    if input_datum.upper() == 'NAD88':
         for region, (correction_path, region_bbox) in NAD88.items():
             if bbox.intersects(region_bbox):
                 return correction_path
@@ -196,7 +196,7 @@ def get_correction_geoid(bbox: Polygon, input_datum: str) -> str:
 
 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.lower() == 'wgs84':
+    if input_datum.upper() == 'WGS84':
         return
     dem_info = gdal.Info(str(dem_file), format='json')
     dem_bbox = get_bbox_from_info(dem_info)
@@ -246,14 +246,14 @@ def prep_dem(input_path: Path, output_path: Path, input_datum: str) -> None:
             resampleAlg='cubic',
             multithread=True,
         )
-    convert_to_height_above_ellipsoid(output_path, input_datum.lower())
+    convert_to_height_above_ellipsoid(output_path, input_datum.upper())
 
 
 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_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).'
+        'input_datum', choices=VALID_DATUMS, help='Datum of the input DEM (supported: WGS84, EGM2008, NAD88).'
     )
     return parser
 

From c31ae8867c89db29c43926bc37f3b841da944157 Mon Sep 17 00:00:00 2001
From: Forrest Williams <ffwilliams2@alaska.edu>
Date: Tue, 19 Aug 2025 11:36:02 -0500
Subject: [PATCH 08/11] fix ruff

---
 src/multirtc/dem.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/multirtc/dem.py b/src/multirtc/dem.py
index aad7d14..58bd957 100644
--- a/src/multirtc/dem.py
+++ b/src/multirtc/dem.py
@@ -1,4 +1,5 @@
 """Prepare an external DEM for use in MultiRTC"""
+
 from concurrent.futures import ThreadPoolExecutor
 from itertools import product
 from pathlib import Path

From 27f8f6a5d3f3a21b1fa7fd384a36391a1e6c7722 Mon Sep 17 00:00:00 2001
From: Forrest Williams <ffwilliams2@alaska.edu>
Date: Fri, 22 Aug 2025 09:47:20 -0500
Subject: [PATCH 09/11] update readme

---
 README.md | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/README.md b/README.md
index a3b415a..cd903cb 100644
--- a/README.md
+++ b/README.md
@@ -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
 

From 4750589ccc467585f73300d63d047ca9d19074e2 Mon Sep 17 00:00:00 2001
From: Forrest Williams <ffwilliams2@alaska.edu>
Date: Fri, 22 Aug 2025 09:53:30 -0500
Subject: [PATCH 10/11] actually use custom dem for processing

---
 src/multirtc/geocode.py  | 5 ++++-
 src/multirtc/multirtc.py | 2 +-
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/multirtc/geocode.py b/src/multirtc/geocode.py
index f12daab..2fea523 100644
--- a/src/multirtc/geocode.py
+++ b/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.dem, args.work_dir, apply_rtc=False)
diff --git a/src/multirtc/multirtc.py b/src/multirtc/multirtc.py
index c8bb43d..d7d70c5 100644
--- a/src/multirtc/multirtc.py
+++ b/src/multirtc/multirtc.py
@@ -111,4 +111,4 @@ def run(args):
         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=True)
+    run_multirtc(args.platform, args.granule, args.resolution, args.dem, args.work_dir, apply_rtc=True)

From 1191482db9b89abf8c46ff6a53f047d5006f233d Mon Sep 17 00:00:00 2001
From: Forrest Williams <ffwilliams2@alaska.edu>
Date: Fri, 22 Aug 2025 17:26:09 -0500
Subject: [PATCH 11/11] fix dem issues

---
 src/multirtc/dem.py      | 33 +++++++++++++++++++++++++++++++--
 src/multirtc/geocode.py  |  2 +-
 src/multirtc/multirtc.py |  2 +-
 3 files changed, 33 insertions(+), 4 deletions(-)

diff --git a/src/multirtc/dem.py b/src/multirtc/dem.py
index 58bd957..95ba69d 100644
--- a/src/multirtc/dem.py
+++ b/src/multirtc/dem.py
@@ -64,7 +64,11 @@ def validate_dem(dem_path: Path, footprint: Polygon) -> None:
 
     dem_extent = get_bbox_from_info(info)
     if not dem_extent.contains(footprint):
-        raise ValueError(f'DEM does not fully cover the footprint: {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]:
@@ -220,12 +224,36 @@ def convert_to_height_above_ellipsoid(dem_file: Path, input_datum) -> None:
 
         dem_ds = gdal.Open(str(dem_file), gdal.GA_Update)
         dem_data = dem_ds.GetRasterBand(1).ReadAsArray()
-        dem_data += geoid_data
+        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.
 
@@ -248,6 +276,7 @@ def prep_dem(input_path: Path, output_path: Path, input_datum: str) -> None:
             multithread=True,
         )
     convert_to_height_above_ellipsoid(output_path, input_datum.upper())
+    set_nodata_value(output_path, 0)
 
 
 def create_parser(parser):
diff --git a/src/multirtc/geocode.py b/src/multirtc/geocode.py
index 2fea523..81cad33 100644
--- a/src/multirtc/geocode.py
+++ b/src/multirtc/geocode.py
@@ -19,4 +19,4 @@ def run(args):
         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.dem, args.work_dir, apply_rtc=False)
+    run_multirtc(args.platform, args.granule, args.resolution, args.work_dir, args.dem, apply_rtc=False)
diff --git a/src/multirtc/multirtc.py b/src/multirtc/multirtc.py
index d7d70c5..7c6a433 100644
--- a/src/multirtc/multirtc.py
+++ b/src/multirtc/multirtc.py
@@ -111,4 +111,4 @@ def run(args):
         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.dem, args.work_dir, apply_rtc=True)
+    run_multirtc(args.platform, args.granule, args.resolution, args.work_dir, args.dem, apply_rtc=True)