Release v0.1.4

.github/dependabot.yml

@@ -11,9 +11,17 @@ updates:
       interval: weekly
     labels:
       - bumpless
+    groups:
+      pip-deps:
+        patterns:
+          - "*"
   - package-ecosystem: github-actions
     directory: /
     schedule:
       interval: weekly
     labels:
       - bumpless
+    groups:
+      github-actions-deps:
+        patterns:
+          - "*"

CHANGELOG.md

@@ -6,6 +6,18 @@ 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.1.4]
+
+### Added
+- Added functionality to `prep_rtc.py` and `upload_rtc.py` to accept SLCs or co-pol bursts.
+- Add `--num-workers` as `prep_rtc` cli param
+
+### Changed
+- Updated DEM fetching/tiling strategy to match OPERA's.
+- Updated DEM bounds buffer to 100 km from 0.025 degrees
+- Updated burst database to the OPERA-provided burst_db_0.2.0_230831-bbox-only.sqlite file
+- Updated our [PGE RunConfig template](./src/hyp3_opera_rtc/templates/pge.yml.j2) to more closely align with the [upstream version](https://github.com/nasa/opera-sds-pcm/blob/9bd74458957197b0c6680540c8d09c26ffab81df/conf/RunConfig.yaml.L2_RTC_S1.jinja2.tmpl).
+
 ## [0.1.3]
 
 ### Changed

Dockerfile

@@ -15,7 +15,7 @@ USER root
 RUN chown rtc_user:rtc_user /home/rtc_user/scratch
 USER rtc_user
 
-RUN curl https://asf-dem-west.s3.amazonaws.com/AUX/opera-burst-bbox-only.sqlite3 -o /home/rtc_user/opera-burst-bbox-only.sqlite3
+RUN curl https://asf-dem-west.s3.amazonaws.com/AUX/burst_db_0.2.0_230831-bbox-only.sqlite -o /home/rtc_user/burst_db_0.2.0_230831-bbox-only.sqlite
 
 COPY --chown=rtc_user:rtc_user . /home/rtc_user/hyp3-opera-rtc/
 RUN conda env create -f /home/rtc_user/hyp3-opera-rtc/environment.yml && \

README.md

@@ -38,7 +38,7 @@ docker run -it --rm \
     [CO_POL_GRANULE]
 ```
 
-Where you replace `[CO_POL_GRANULE]` with the name of the Sentinel-1 co-pol burst SLC scene
+Where you replace `[CO_POL_GRANULE]` with the name of the Sentinel-1 co-pol burst or SLC scene
 for which to generate the OPERA RTC product.
 Here are some useful examples:
 
@@ -49,7 +49,7 @@ Here are some useful examples:
 
 ## Architecture
 
-The plugin is composed of three nested docker environments that depend on eachother. They are laid out as below:
+The plugin is composed of three nested docker environments that depend on each other. They are laid out as below:
 
 ```
 +-------------------------+

requirements-static.txt

@@ -1,5 +1,5 @@
-ruff==0.11.11
-mypy==1.15.0
+ruff==0.12.2
+mypy==1.16.1
 lxml-stubs
 types-shapely
 types-requests

src/hyp3_opera_rtc/dem.py

@@ -1,25 +1,59 @@
-from concurrent.futures import ThreadPoolExecutor
-from itertools import product
+from collections.abc import Callable
 from pathlib import Path
-from tempfile import TemporaryDirectory
+from tempfile import NamedTemporaryFile
 
 import numpy as np
 import shapely
 import shapely.ops
 import shapely.wkt
-from hyp3lib.fetch import download_file
-from osgeo import gdal
-from shapely.geometry import LinearRing, Polygon
+from hyp3lib.util import GDALConfigManager
+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'
 
+EARTH_APPROX_CIRCUMFERENCE_KM = 40075017.0
+EARTH_RADIUS_KM = EARTH_APPROX_CIRCUMFERENCE_KM / (2 * np.pi)
+DEM_MARGIN_KM = 200
 
-def check_antimeridean(poly: Polygon) -> list[Polygon]:
+
+def margin_km_to_deg(margin_in_km: float) -> float:
+    """Converts a margin value from kilometers to degrees."""
+    km_to_deg_at_equator = 1000.0 / (EARTH_APPROX_CIRCUMFERENCE_KM / 360.0)
+    margin_in_deg = margin_in_km * km_to_deg_at_equator
+    return margin_in_deg
+
+
+def margin_km_to_longitude_deg(margin_in_km: float, lat: float) -> float:
+    """Converts a margin value from kilometers to degrees as a function of latitude."""
+    delta_lon = 180 * 1000 * margin_in_km / (np.pi * EARTH_RADIUS_KM * np.cos(np.pi * lat / 180))
+    return delta_lon
+
+
+def polygon_from_bounds(bounds: tuple[float, float, float, float]) -> Polygon:
+    """Create a polygon (EPSG:4326) from the lat/lon coordinates corresponding to a provided bounding box."""
+    lon_min, lat_min, lon_max, lat_max = bounds
+    # note we can also use the center lat here
+    lat_worst_case = max([lat_min, lat_max])
+    lat_margin = margin_km_to_deg(DEM_MARGIN_KM)
+    lon_margin = margin_km_to_longitude_deg(DEM_MARGIN_KM, lat=lat_worst_case)
+    # Check if the bbox crosses the antimeridian and apply the margin accordingly
+    # so that any resultant DEM is split properly by check_dateline
+    if lon_max - lon_min > 180:
+        lon_min, lon_max = lon_max, lon_min
+
+    poly = box(
+        lon_min - lon_margin, max([lat_min - lat_margin, -90]), lon_max + lon_margin, min([lat_max + lat_margin, 90])
+    )
+    return poly
+
+
+def split_antimeridian(poly: Polygon) -> list[Polygon]:
+    """Check if the provided polygon crosses the antimeridian and split it if it does."""
     x_min, _, x_max, _ = poly.bounds
 
-    # Check anitmeridean crossing
+    # Check anitmeridian crossing
     if (x_max - x_min > 180.0) or (x_min <= 180.0 <= x_max):
         dateline = shapely.wkt.loads('LINESTRING( 180.0 -90.0, 180.0 90.0)')
 
@@ -54,38 +88,93 @@ def check_antimeridean(poly: Polygon) -> list[Polygon]:
     return polys
 
 
-def get_dem_granule_url(lat: int, lon: int) -> str:
-    lat_tens = np.floor_divide(lat, 10) * 10
-    lat_cardinal = 'S' if lat_tens < 0 else 'N'
-
-    lon_tens = np.floor_divide(lon, 20) * 20
-    lon_cardinal = 'W' if lon_tens < 0 else 'E'
-
-    prefix = f'{lat_cardinal}{np.abs(lat_tens):02d}_{lon_cardinal}{np.abs(lon_tens):03d}'
-    filename = f'DEM_{lat_cardinal}{np.abs(lat):02d}_00_{lon_cardinal}{np.abs(lon):03d}_00.tif'
-    file_url = f'{URL}/{prefix}/{filename}'
-    return file_url
-
-
-def get_latlon_pairs(polygon: Polygon) -> list:
-    minx, miny, maxx, maxy = polygon.bounds
-    lats = np.arange(np.floor(miny), np.floor(maxy) + 1).astype(int)
-    lons = np.arange(np.floor(minx), np.floor(maxx) + 1).astype(int)
-    return list(product(lats, lons))
-
-
-def download_opera_dem_for_footprint(output_path: Path, footprint: Polygon) -> None:
-    footprints = check_antimeridean(footprint)
-    latlon_pairs = []
-    for footprint in footprints:
-        latlon_pairs += get_latlon_pairs(footprint)
-    urls = [get_dem_granule_url(lat, lon) for lat, lon in latlon_pairs]
-
-    with TemporaryDirectory() as tmpdir_str:
-        tmpdir = Path(tmpdir_str)
-        with ThreadPoolExecutor(max_workers=4) as executor:
-            executor.map(lambda url: download_file(url, str(tmpdir)), urls)
-        vrt_filepath = str(tmpdir / 'dem.vrt')
-        input_files = [str(file) for file in tmpdir.glob('*.tif')]
-        gdal.BuildVRT(vrt_filepath, input_files)
-        gdal.Translate(str(output_path), vrt_filepath, format='GTiff')
+def snap_coord(val: float, snap: float, offset: float, round_func: Callable) -> float:
+    return round_func(float(val - offset) / snap) * snap + offset
+
+
+def translate_dem(vrt_filename: str, output_path: str, bounds: tuple[float, float, float, float]) -> None:
+    """Write a local subset of the OPERA DEM for a region matching the provided bounds.
+
+    Params:
+        vrt_filename: Path to the input VRT file
+        output_path: Path to the translated output GTiff file
+        bounds: Bounding box in the form of (lon_min, lat_min, lon_max, lat_max)
+    """
+    ds = gdal.Open(vrt_filename, gdal.GA_ReadOnly)
+
+    # update cropping coordinates to not exceed the input DEM bounding box
+    input_x_min, xres, _, input_y_max, _, yres = ds.GetGeoTransform()
+    length = ds.GetRasterBand(1).YSize
+    width = ds.GetRasterBand(1).XSize
+
+    # Snap edge coordinates using the DEM pixel spacing
+    # (xres and yres) and starting coordinates (input_x_min and
+    # input_x_max). Maximum values are rounded using np.ceil
+    # and minimum values are rounded using np.floor
+    x_min, y_min, x_max, y_max = bounds
+    snapped_x_min = snap_coord(x_min, xres, input_x_min, np.floor)
+    snapped_x_max = snap_coord(x_max, xres, input_x_min, np.ceil)
+    snapped_y_min = snap_coord(y_min, yres, input_y_max, np.floor)
+    snapped_y_max = snap_coord(y_max, yres, input_y_max, np.ceil)
+
+    input_y_min = input_y_max + length * yres
+    input_x_max = input_x_min + width * xres
+
+    adjusted_x_min = max(snapped_x_min, input_x_min)
+    adjusted_x_max = min(snapped_x_max, input_x_max)
+    adjusted_y_min = max(snapped_y_min, input_y_min)
+    adjusted_y_max = min(snapped_y_max, input_y_max)
+
+    try:
+        gdal.Translate(
+            output_path, ds, format='GTiff', projWin=[adjusted_x_min, adjusted_y_max, adjusted_x_max, adjusted_y_min]
+        )
+    except RuntimeError as err:
+        if 'negative width and/or height' in str(err):
+            gdal.Translate(output_path, ds, format='GTiff', projWin=[x_min, y_max, x_max, y_min])
+        else:
+            raise
+
+    # stage_dem.py takes a bbox as an input. The longitude coordinates
+    # of this bbox are unwrapped i.e., range in [0, 360] deg. If the
+    # bbox crosses the anti-meridian, the script divides it in two
+    # bboxes neighboring the anti-meridian. Here, x_min and x_max
+    # represent the min and max longitude coordinates of one of these
+    # bboxes. We Add 360 deg if the min longitude of the downloaded DEM
+    # tile is < 180 deg i.e., there is a dateline crossing.
+    # This ensures that the mosaicked DEM VRT will span a min
+    # range of longitudes rather than the full [-180, 180] deg
+    sr = osr.SpatialReference(ds.GetProjection())
+    epsg_str = sr.GetAttrValue('AUTHORITY', 1)
+
+    if x_min <= -180.0 and epsg_str == '4326':
+        ds = gdal.Open(output_path, gdal.GA_Update)
+        geotransform = list(ds.GetGeoTransform())
+        geotransform[0] += 360.0
+        ds.SetGeoTransform(tuple(geotransform))
+
+
+def download_opera_dem_for_footprint(outfile: Path, bounds: tuple[float, float, float, float]) -> None:
+    """Download a DEM from the specified S3 bucket.
+
+    Params:
+        outfile: Path to the where the output DEM file is to be staged.
+        bounds: Bounding box in the form of (lon_min, lat_min, lon_max, lat_max).
+    """
+    poly = polygon_from_bounds(bounds)
+    polys = split_antimeridian(poly)
+    dem_list = []
+
+    with NamedTemporaryFile(suffix='.txt') as cookie_file:
+        with GDALConfigManager(
+            GDAL_HTTP_COOKIEJAR=cookie_file.name,
+            GDAL_HTTP_COOKIEFILE=cookie_file.name,
+            GDAL_DISABLE_READDIR_ON_OPEN='EMPTY_DIR',
+        ):
+            vrt_filename = '/vsicurl/https://nisar.asf.earthdatacloud.nasa.gov/STATIC/DEM/v1.1/EPSG4326/EPSG4326.vrt'
+            for idx, poly in enumerate(polys):
+                output_path = str(outfile.parent / f'{outfile.stem}_{idx}.tif')
+                dem_list.append(output_path)
+                translate_dem(vrt_filename, output_path, poly.bounds)
+
+            gdal.BuildVRT(str(outfile), dem_list)