Repository: ASFHyP3/hyp3-autorift

Issue:
Get associated autoRIFT files from parameter shapefile
Currently, we hardcode the set of needed files for processing: 

https://github.com/ASFHyP3/hyp3-autorift/blob/develop/hyp3_autorift/io.py#L47-L68

These however, at detailed in the new parameter shapefile and could be pulled from there instead to be more flexible to future change (see #47 )


Discussion:

Patch:
diff --git a/hyp3_autorift/geometry.py b/hyp3_autorift/geometry.py
index 5b66d76..a2c6a13 100644
--- a/hyp3_autorift/geometry.py
+++ b/hyp3_autorift/geometry.py
@@ -2,19 +2,18 @@
 
 import logging
 import os
+from typing import Tuple
 
 import isce  # noqa: F401
 import isceobj
 import numpy as np
 from contrib.demUtils import createDemStitcher
 from contrib.geo_autoRIFT.geogrid import Geogrid
-from hyp3lib import DemError
 from isceobj.Orbit.Orbit import Orbit
 from isceobj.Sensor.TOPS.Sentinel1 import Sentinel1
 from osgeo import gdal
 from osgeo import ogr
-
-from hyp3_autorift.io import AUTORIFT_PREFIX, ITS_LIVE_BUCKET
+from osgeo import osr
 
 log = logging.getLogger(__name__)
 
@@ -89,7 +88,7 @@ def bounding_box(safe, priority='reference', polarization='hh', orbits='Orbits',
     return lat_limits, lon_limits
 
 
-def polygon_from_bbox(lat_limits, lon_limits):
+def polygon_from_bbox(lat_limits: Tuple[float, float], lon_limits: Tuple[float, float]) -> ogr.Geometry:
     ring = ogr.Geometry(ogr.wkbLinearRing)
     ring.AddPoint(lon_limits[0], lat_limits[0])
     ring.AddPoint(lon_limits[1], lat_limits[0])
@@ -101,20 +100,24 @@ def polygon_from_bbox(lat_limits, lon_limits):
     return polygon
 
 
-def find_jpl_dem(lat_limits, lon_limits):
-    shape_file = f'/vsicurl/http://{ITS_LIVE_BUCKET}.s3.amazonaws.com/{AUTORIFT_PREFIX}/autorift_parameters.shp'
-    driver = ogr.GetDriverByName('ESRI Shapefile')
-    shapes = driver.Open(shape_file, gdal.GA_ReadOnly)
+def poly_bounds_in_proj(polygon: ogr.Geometry, in_epsg: int, out_epsg: int):
+    in_srs = osr.SpatialReference()
+    in_srs.ImportFromEPSG(in_epsg)
 
-    centroid = polygon_from_bbox(lat_limits, lon_limits).Centroid()
-    for feature in shapes.GetLayer(0):
-        if feature.geometry().Contains(centroid):
-            return f'{feature["name"]}_0240m'
+    out_srs = osr.SpatialReference()
+    out_srs.ImportFromEPSG(out_epsg)
 
-    raise DemError('Could not determine appropriate DEM for:\n'
-                   f'    lat (min, max): {lat_limits}'
-                   f'    lon (min, max): {lon_limits}'
-                   f'    using: {shape_file}')
+    transformation = osr.CoordinateTransformation(in_srs, out_srs)
+    ring = ogr.Geometry(ogr.wkbLinearRing)
+    for point in polygon.GetBoundary().GetPoints():
+        try:
+            lon, lat = point
+        except ValueError:
+            # buffered polygons only have (X, Y) while unbuffered have (X, Y, Z)
+            lon, lat, _ = point
+        ring.AddPoint(*transformation.TransformPoint(lat, lon))
+
+    return ring.GetEnvelope()
 
 
 def prep_isce_dem(input_dem, lat_limits, lon_limits, isce_dem=None):
diff --git a/hyp3_autorift/io.py b/hyp3_autorift/io.py
index 226f3d6..b38fa23 100644
--- a/hyp3_autorift/io.py
+++ b/hyp3_autorift/io.py
@@ -3,24 +3,24 @@
 import argparse
 import logging
 import os
-import shutil
 import textwrap
 from pathlib import Path
 from typing import Union
 
 import boto3
-from boto3.s3.transfer import TransferConfig
-from botocore import UNSIGNED
-from botocore.config import Config
+from hyp3lib import DemError
 from isce.applications.topsApp import TopsInSAR
+from osgeo import gdal
+from osgeo import ogr
 from scipy.io import savemat
 
+from hyp3_autorift.geometry import poly_bounds_in_proj
+
 log = logging.getLogger(__name__)
 
 ITS_LIVE_BUCKET = 'its-live-data.jpl.nasa.gov'
 AUTORIFT_PREFIX = 'autorift_parameters/v001'
 
-_s3_client_unsigned = boto3.client('s3', config=Config(signature_version=UNSIGNED))
 _s3_client = boto3.client('s3')
 
 
@@ -31,53 +31,65 @@ def download_s3_file_requester_pays(target_path: Union[str, Path], bucket: str,
     return filename
 
 
-def _download_s3_files(target_dir, bucket, keys, chunk_size=50*1024*1024):
-    transfer_config = TransferConfig(multipart_threshold=chunk_size, multipart_chunksize=chunk_size)
-    file_list = []
-    for key in keys:
-        filename = os.path.join(target_dir, os.path.basename(key))
-        if os.path.exists(filename):
-            continue
-        file_list.append(filename)
-        log.info(f'Downloading s3://{bucket}/{key} to {filename}')
-        _s3_client_unsigned.download_file(Bucket=bucket, Key=key, Filename=filename, Config=transfer_config)
-    return file_list
-
-
-def _get_s3_keys_for_dem(prefix=AUTORIFT_PREFIX, dem='GRE240m'):
-    tags = [
-        'h',
-        'StableSurface',
-        'dhdx',
-        'dhdy',
-        'dhdxs',
-        'dhdys',
-        'vx0',
-        'vy0',
-        'vxSearchRange',
-        'vySearchRange',
-        'xMinChipSize',
-        'yMinChipSize',
-        'xMaxChipSize',
-        'yMaxChipSize',
-        # FIXME: Was renamed from masks to sp by JPL; change hasn't been propagated to autoRIFT
-        #        keep last so we can easily rename the file after downloading
-        'sp',
-    ]
-    keys = [f'{prefix}/{dem}_{tag}.tif' for tag in tags]
-    return keys
-
-
-def fetch_jpl_tifs(dem='GRE240m', target_dir='DEM', bucket=ITS_LIVE_BUCKET, prefix=AUTORIFT_PREFIX):
-    # FIXME: gdalwarp  needed subset instead?
-    log.info(f"Downloading {dem} tifs from JPL's AWS bucket")
-
-    for logger in ('botocore', 's3transfer'):
-        logging.getLogger(logger).setLevel(logging.WARNING)
-
-    keys = _get_s3_keys_for_dem(prefix, dem)
-    tifs = _download_s3_files(target_dir, bucket, keys)
-    shutil.move(tifs[-1], tifs[-1].replace('_sp.tif', '_masks.tif'))
+def find_jpl_dem(polygon: ogr.Geometry) -> dict:
+    shape_file = f'/vsicurl/http://{ITS_LIVE_BUCKET}.s3.amazonaws.com/{AUTORIFT_PREFIX}/autorift_parameters.shp'
+    driver = ogr.GetDriverByName('ESRI Shapefile')
+    shapes = driver.Open(shape_file, gdal.GA_ReadOnly)
+
+    centroid = polygon.Centroid()
+    for feature in shapes.GetLayer(0):
+        if feature.geometry().Contains(centroid):
+            dem_info = {
+                'name': f'{feature["name"]}_0240m',
+                'epsg': feature['epsg'],
+                'tifs': {
+                    'h': f"/vsicurl/{feature['h']}",
+                    'StableSurface': f"/vsicurl/{feature['StableSurfa']}",
+                    'dhdx': f"/vsicurl/{feature['dhdx']}",
+                    'dhdy': f"/vsicurl/{feature['dhdy']}",
+                    'dhdxs': f"/vsicurl/{feature['dhdxs']}",
+                    'dhdys': f"/vsicurl/{feature['dhdys']}",
+                    'vx0': f"/vsicurl/{feature['vx0']}",
+                    'vy0': f"/vsicurl/{feature['vy0']}",
+                    'vxSearchRange': f"/vsicurl/{feature['vxSearchRan']}",
+                    'vySearchRange': f"/vsicurl/{feature['vySearchRan']}",
+                    'xMinChipSize': f"/vsicurl/{feature['xMinChipSiz']}",
+                    'yMinChipSize': f"/vsicurl/{feature['yMinChipSiz']}",
+                    'xMaxChipSize': f"/vsicurl/{feature['xMaxChipSiz']}",
+                    'yMaxChipSize': f"/vsicurl/{feature['yMaxChipSiz']}",
+                    'sp': f"/vsicurl/{feature['sp']}",
+                },
+            }
+            return dem_info
+
+    raise DemError('Could not determine appropriate DEM for:\n'
+                   f'    centroid: {centroid}'
+                   f'    using: {shape_file}')
+
+
+def subset_jpl_tifs(polygon: ogr.Geometry, buffer: float = 0.15, target_dir: Union[str, Path] = '.'):
+    dem_info = find_jpl_dem(polygon)
+    log.info(f'Subsetting {dem_info["name"]} tifs from s3://{ITS_LIVE_BUCKET}/{AUTORIFT_PREFIX}/')
+
+    min_x, max_x, min_y, max_y = poly_bounds_in_proj(polygon.Buffer(buffer), in_epsg=4326, out_epsg=dem_info['epsg'])
+    output_bounds = (min_x, min_y, max_x, max_y)
+    log.debug(f'Subset bounds: {output_bounds}')
+
+    subset_tifs = {}
+    for key, tif in dem_info['tifs'].items():
+        out_path = os.path.join(target_dir, os.path.basename(tif))
+
+        # FIXME: shouldn't need to do after next autoRIFT upgrade
+        if out_path.endswith('_sp.tif'):
+            out_path = out_path.replace('_sp.tif', '_masks.tif')
+
+        subset_tifs[key] = out_path
+
+        gdal.Warp(
+            out_path, tif, outputBounds=output_bounds, xRes=240, yRes=240, targetAlignedPixels=True, multithread=True,
+        )
+
+    return subset_tifs
 
 
 def format_tops_xml(reference, secondary, polarization, dem, orbits, xml_file='topsApp.xml'):
diff --git a/hyp3_autorift/process.py b/hyp3_autorift/process.py
index 401bd34..d3d3936 100644
--- a/hyp3_autorift/process.py
+++ b/hyp3_autorift/process.py
@@ -196,25 +196,22 @@ def process(reference: str, secondary: str, polarization: str = 'hh', band: str
         lat_limits = (bbox[1], bbox[3])
         lon_limits = (bbox[0], bbox[2])
 
-    dem = geometry.find_jpl_dem(lat_limits, lon_limits)
-    dem_dir = os.path.join(os.getcwd(), 'DEM')
-    mkdir_p(dem_dir)
-    io.fetch_jpl_tifs(dem=dem, target_dir=dem_dir)
-    dem_prefix = os.path.join(dem_dir, dem)
-
-    geogrid_parameters = f'-d {dem_prefix}_h.tif -ssm {dem_prefix}_StableSurface.tif ' \
-                         f'-sx {dem_prefix}_dhdx.tif -sy {dem_prefix}_dhdy.tif ' \
-                         f'-vx {dem_prefix}_vx0.tif -vy {dem_prefix}_vy0.tif ' \
-                         f'-srx {dem_prefix}_vxSearchRange.tif -sry {dem_prefix}_vySearchRange.tif ' \
-                         f'-csminx {dem_prefix}_xMinChipSize.tif -csminy {dem_prefix}_yMinChipSize.tif ' \
-                         f'-csmaxx {dem_prefix}_xMaxChipSize.tif -csmaxy {dem_prefix}_yMaxChipSize.tif'
+    scene_poly = geometry.polygon_from_bbox(lat_limits, lon_limits)
+    tifs = io.subset_jpl_tifs(scene_poly, target_dir=Path.cwd())
+
+    geogrid_parameters = f'-d {tifs["h"]} -ssm {tifs["StableSurface"]} ' \
+                         f'-sx {tifs["dhdx"]} -sy {tifs["dhdy"]} ' \
+                         f'-vx {tifs["vx0"]} -vy {tifs["vy0"]} ' \
+                         f'-srx {tifs["vxSearchRange"]} -sry {tifs["vySearchRange"]} ' \
+                         f'-csminx {tifs["xMinChipSize"]} -csminy {tifs["yMinChipSize"]} ' \
+                         f'-csmaxx {tifs["xMaxChipSize"]} -csmaxy {tifs["yMaxChipSize"]}'
     autorift_parameters = '-g window_location.tif -o window_offset.tif -sr window_search_range.tif ' \
                           '-csmin window_chip_size_min.tif -csmax window_chip_size_max.tif ' \
                           '-vx window_rdr_off2vel_x_vec.tif -vy window_rdr_off2vel_y_vec.tif ' \
                           '-ssm window_stable_surface_mask.tif'
 
     if platform == 'S1':
-        isce_dem = geometry.prep_isce_dem(f'{dem_prefix}_h.tif', lat_limits, lon_limits)
+        isce_dem = geometry.prep_isce_dem(tifs["h"], lat_limits, lon_limits)
 
         io.format_tops_xml(reference, secondary, polarization, isce_dem, orbits)
 
diff --git a/setup.py b/setup.py
index 54a17a4..96de391 100644
--- a/setup.py
+++ b/setup.py
@@ -62,7 +62,6 @@ setup(
         'testGeogrid_ISCE.py = hyp3_autorift.vend.testGeogrid_ISCE:main',
         'testautoRIFT.py = hyp3_autorift.vend.testautoRIFT:main',
         'testGeogridOptical.py = hyp3_autorift.vend.testGeogridOptical:main',
-        # FIXME: Only needed for testautoRIFT.py and testautoRIFT_ISCE.py
         'topsinsar_filename.py = hyp3_autorift.io:topsinsar_mat',
         ]
     },


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index b4399a0..6b44082 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,7 +1,7 @@
 import pytest
 from botocore.stub import Stubber
 
-from hyp3_autorift.io import _s3_client, _s3_client_unsigned
+from hyp3_autorift.io import _s3_client
 
 
 @pytest.fixture
@@ -9,10 +9,3 @@ def s3_stub():
     with Stubber(_s3_client) as stubber:
         yield stubber
         stubber.assert_no_pending_responses()
-
-
-@pytest.fixture
-def s3_unsigned_stub():
-    with Stubber(_s3_client_unsigned) as stubber:
-        yield stubber
-        stubber.assert_no_pending_responses()
diff --git a/tests/test_geometry.py b/tests/test_geometry.py
index 264f07c..54d444e 100644
--- a/tests/test_geometry.py
+++ b/tests/test_geometry.py
@@ -1,5 +1,4 @@
-import pytest
-from hyp3lib import DemError
+import numpy as np
 
 from hyp3_autorift import geometry
 
@@ -11,47 +10,27 @@ def test_polygon_from_bbox():
            == 'POLYGON ((3 1 0,4 1 0,4 2 0,3 2 0,3 1 0))'
 
 
-def test_find_jpl_dem():
-    lat_limits = (55, 56)
-    lon_limits = (40, 41)
-    assert geometry.find_jpl_dem(lat_limits, lon_limits) == 'NPS_0240m'
-
-    lat_limits = (54, 55)
-    lon_limits = (40, 41)
-    assert geometry.find_jpl_dem(lat_limits, lon_limits) == 'N37_0240m'
-
-    lat_limits = (54, 55)
-    lon_limits = (-40, -41)
-    assert geometry.find_jpl_dem(lat_limits, lon_limits) == 'N24_0240m'
-
-    lat_limits = (-54, -55)
-    lon_limits = (-40, -41)
-    assert geometry.find_jpl_dem(lat_limits, lon_limits) == 'S24_0240m'
-
-    lat_limits = (-55, -56)
-    lon_limits = (40, 41)
-    assert geometry.find_jpl_dem(lat_limits, lon_limits) == 'S37_0240m'
-
-    lat_limits = (-56, -57)
-    lon_limits = (40, 41)
-    assert geometry.find_jpl_dem(lat_limits, lon_limits) == 'SPS_0240m'
-
-    lat_limits = (-90, -91)
-    lon_limits = (40, 41)
-    with pytest.raises(DemError):
-        geometry.find_jpl_dem(lat_limits, lon_limits)
-
-    lat_limits = (90, 91)
-    lon_limits = (40, 41)
-    with pytest.raises(DemError):
-        geometry.find_jpl_dem(lat_limits, lon_limits)
-
-    lat_limits = (55, 56)
-    lon_limits = (180, 181)
-    with pytest.raises(DemError):
-        geometry.find_jpl_dem(lat_limits, lon_limits)
-
-    lat_limits = (55, 56)
-    lon_limits = (-180, -181)
-    with pytest.raises(DemError):
-        geometry.find_jpl_dem(lat_limits, lon_limits)
+def test_pol_bounds_in_proj():
+    polygon = geometry.polygon_from_bbox(lat_limits=(55, 56), lon_limits=(40, 41))
+    assert np.allclose(
+        geometry.poly_bounds_in_proj(polygon, in_epsg=4326, out_epsg=3413),  # NPS
+        (3776365.5697414433, 3899644.3388010086, -340706.3423259673, -264432.19003121805)
+    )
+
+    polygon = geometry.polygon_from_bbox(lat_limits=(-58, -57), lon_limits=(40, 41))
+    assert np.allclose(
+        geometry.poly_bounds_in_proj(polygon, in_epsg=4326, out_epsg=3031),  # SPS
+        (2292512.6214329866, 2416952.768825992, 2691684.1770189586, 2822144.2827928355)
+    )
+
+    polygon = geometry.polygon_from_bbox(lat_limits=(22, 23), lon_limits=(40, 41))
+    assert np.allclose(
+        geometry.poly_bounds_in_proj(polygon, in_epsg=4326, out_epsg=32637),
+        (602485.1663686256, 706472.0593133729, 2433164.428653589, 2544918.1043369616)
+    )
+
+    polygon = geometry.polygon_from_bbox(lat_limits=(-23, -22), lon_limits=(40, 41))
+    assert np.allclose(
+        geometry.poly_bounds_in_proj(polygon, in_epsg=4326, out_epsg=32737),
+        (602485.1663686256, 706472.0593133729, 7455081.895663038, 7566835.5713464115)
+    )
diff --git a/tests/test_io.py b/tests/test_io.py
index 5d75942..cb4d334 100644
--- a/tests/test_io.py
+++ b/tests/test_io.py
@@ -1,6 +1,9 @@
 from io import BytesIO
 
-from hyp3_autorift import io
+import pytest
+from hyp3lib import DemError
+
+from hyp3_autorift import geometry, io
 
 
 def test_download_s3_file_requester_pays(tmp_path, s3_stub):
@@ -21,52 +24,43 @@ def test_download_s3_file_requester_pays(tmp_path, s3_stub):
     assert tmp_path / 'foobar.txt' == file
 
 
-def test_get_s3_keys_for_dem():
-    expected = [
-        'Prefix/GRE240m_h.tif',
-        'Prefix/GRE240m_StableSurface.tif',
-        'Prefix/GRE240m_dhdx.tif',
-        'Prefix/GRE240m_dhdy.tif',
-        'Prefix/GRE240m_dhdxs.tif',
-        'Prefix/GRE240m_dhdys.tif',
-        'Prefix/GRE240m_vx0.tif',
-        'Prefix/GRE240m_vy0.tif',
-        'Prefix/GRE240m_vxSearchRange.tif',
-        'Prefix/GRE240m_vySearchRange.tif',
-        'Prefix/GRE240m_xMinChipSize.tif',
-        'Prefix/GRE240m_yMinChipSize.tif',
-        'Prefix/GRE240m_xMaxChipSize.tif',
-        'Prefix/GRE240m_yMaxChipSize.tif',
-        'Prefix/GRE240m_sp.tif',
-    ]
-    assert sorted(io._get_s3_keys_for_dem('Prefix', 'GRE240m')) == sorted(expected)
+def test_find_jpl_dem():
+    polygon = geometry.polygon_from_bbox(lat_limits=(55, 56), lon_limits=(40, 41))
+    dem_info = io.find_jpl_dem(polygon)
+    assert dem_info['name'] == 'NPS_0240m'
+
+    polygon = geometry.polygon_from_bbox(lat_limits=(54, 55), lon_limits=(40, 41))
+    dem_info = io.find_jpl_dem(polygon)
+    assert dem_info['name'] == 'N37_0240m'
+
+    polygon = geometry.polygon_from_bbox(lat_limits=(54, 55), lon_limits=(-40, -41))
+    dem_info = io.find_jpl_dem(polygon)
+    assert dem_info['name'] == 'N24_0240m'
+
+    polygon = geometry.polygon_from_bbox(lat_limits=(-54, -55), lon_limits=(-40, -41))
+    dem_info = io.find_jpl_dem(polygon)
+    assert dem_info['name'] == 'S24_0240m'
+
+    polygon = geometry.polygon_from_bbox(lat_limits=(-55, -56), lon_limits=(40, 41))
+    dem_info = io.find_jpl_dem(polygon)
+    assert dem_info['name'] == 'S37_0240m'
+
+    polygon = geometry.polygon_from_bbox(lat_limits=(-56, -57), lon_limits=(40, 41))
+    dem_info = io.find_jpl_dem(polygon)
+    assert dem_info['name'] == 'SPS_0240m'
+
+    polygon = geometry.polygon_from_bbox(lat_limits=(-90, -91), lon_limits=(40, 41))
+    with pytest.raises(DemError):
+        io.find_jpl_dem(polygon)
+
+    polygon = geometry.polygon_from_bbox(lat_limits=(90, 91), lon_limits=(40, 41))
+    with pytest.raises(DemError):
+        io.find_jpl_dem(polygon)
 
+    polygon = geometry.polygon_from_bbox(lat_limits=(55, 56), lon_limits=(180, 181))
+    with pytest.raises(DemError):
+        io.find_jpl_dem(polygon)
 
-def test_download_s3_files(tmp_path, s3_unsigned_stub):
-    keys = ['foo', 'bar']
-    for key in keys:
-        s3_unsigned_stub.add_response(
-            'head_object',
-            expected_params={
-                'Bucket': 'myBucket',
-                'Key': key,
-            },
-            service_response={
-                'ContentLength': 3,
-            },
-        )
-        s3_unsigned_stub.add_response(
-            'get_object',
-            expected_params={
-                'Bucket': 'myBucket',
-                'Key': key,
-            },
-            service_response={
-                'Body': BytesIO(b'123'),
-            },
-        )
-    downloaded_files = io._download_s3_files(tmp_path, 'myBucket', keys)
-    for key, downloaded_file in zip(keys, downloaded_files):
-        assert (tmp_path / key).exists()
-        assert (tmp_path / key).read_text() == '123'
-        assert str(tmp_path / key) == downloaded_file
+    polygon = geometry.polygon_from_bbox(lat_limits=(55, 56), lon_limits=(-180, -181))
+    with pytest.raises(DemError):
+        io.find_jpl_dem(polygon)