Merge pull request #59 from ASFHyP3/develop

Release L8 IndexError (partial) fix

CHANGELOG.md

@@ -6,6 +6,17 @@ 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.4.2](https://github.com/ASFHyP3/hyp3-autorift/compare/v0.4.1...v0.4.2)
+
+### Fixed
+* A *partial* fix was implemented to correct out of index errors when processing
+  optical scenes (typically seen with Landsat-8 pairs) due to calculating different
+  overlapping subset sizes when co-registering the images. Currently, only the
+  smallest subset size is used, so the bounding box may be 1px too small in x
+  and/or y, but there shouldn't be any pixel offsets. Full fix will need to be
+  implemented upstream in [autoRIFT](https://github.com/leiyangleon/autoRIFT).
+
 ## [0.4.1](https://github.com/ASFHyP3/hyp3-autorift/compare/v0.4.0...v0.4.1)
 
 ### Changed

hyp3_autorift/geometry.py

@@ -88,36 +88,49 @@ def bounding_box(safe, priority='reference', polarization='hh', orbits='Orbits',
     return lat_limits, lon_limits
 
 
-def polygon_from_bbox(lat_limits: Tuple[float, float], lon_limits: Tuple[float, float]) -> ogr.Geometry:
+def polygon_from_bbox(x_limits: Tuple[float, float], y_limits: Tuple[float, float],
+                      epsg_code: int = 4326) -> ogr.Geometry:
     ring = ogr.Geometry(ogr.wkbLinearRing)
-    ring.AddPoint(lon_limits[0], lat_limits[0])
-    ring.AddPoint(lon_limits[1], lat_limits[0])
-    ring.AddPoint(lon_limits[1], lat_limits[1])
-    ring.AddPoint(lon_limits[0], lat_limits[1])
-    ring.AddPoint(lon_limits[0], lat_limits[0])
+    ring.AddPoint_2D(x_limits[0], y_limits[1])
+    ring.AddPoint_2D(x_limits[1], y_limits[1])
+    ring.AddPoint_2D(x_limits[1], y_limits[0])
+    ring.AddPoint_2D(x_limits[0], y_limits[0])
+    ring.AddPoint_2D(x_limits[0], y_limits[1])
     polygon = ogr.Geometry(ogr.wkbPolygon)
     polygon.AddGeometry(ring)
+
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(epsg_code)
+    polygon.AssignSpatialReference(srs)
+
     return polygon
 
 
-def poly_bounds_in_proj(polygon: ogr.Geometry, in_epsg: int, out_epsg: int):
-    in_srs = osr.SpatialReference()
-    in_srs.ImportFromEPSG(in_epsg)
+def poly_bounds_in_proj(polygon: ogr.Geometry, out_epsg: int):
+    in_srs = polygon.GetSpatialReference()
+    if in_srs is None:
+        log.warning('Polygon does not have an assigned spatial reference; assuming EPSG:4326.')
+        in_srs = osr.SpatialReference()
+        in_srs.ImportFromEPSG(4326)
 
     out_srs = osr.SpatialReference()
     out_srs.ImportFromEPSG(out_epsg)
 
     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()
+    out_poly = ogr.Geometry(wkb=polygon.ExportToWkb())
+    out_poly.Transform(transformation)
+
+    return out_poly.GetEnvelope()
+
+
+def flip_point_coordinates(point: ogr.Geometry):
+    if not point.GetGeometryName() == 'POINT':
+        raise ValueError('Can only flip POINT geometries')
+
+    flipped = ogr.Geometry(ogr.wkbPoint)
+    flipped.AddPoint_2D(point.GetY(), point.GetX())
+
+    return flipped
 
 
 def prep_isce_dem(input_dem, lat_limits, lon_limits, isce_dem=None):

hyp3_autorift/io.py

@@ -14,7 +14,7 @@
 from osgeo import ogr
 from scipy.io import savemat
 
-from hyp3_autorift.geometry import poly_bounds_in_proj
+from hyp3_autorift.geometry import flip_point_coordinates, poly_bounds_in_proj
 
 log = logging.getLogger(__name__)
 
@@ -36,7 +36,7 @@ def find_jpl_dem(polygon: ogr.Geometry) -> dict:
     driver = ogr.GetDriverByName('ESRI Shapefile')
     shapes = driver.Open(shape_file, gdal.GA_ReadOnly)
 
-    centroid = polygon.Centroid()
+    centroid = flip_point_coordinates(polygon.Centroid())
     for feature in shapes.GetLayer(0):
         if feature.geometry().Contains(centroid):
             dem_info = {
@@ -71,7 +71,7 @@ def subset_jpl_tifs(polygon: ogr.Geometry, buffer: float = 0.15, target_dir: Uni
     dem_info = find_jpl_dem(polygon)
     log.info(f'Subsetting {dem_info["name"]} tifs: {dem_info["tifs"]["h"].replace("_h.tif", "_*")}')
 
-    min_x, max_x, min_y, max_y = poly_bounds_in_proj(polygon.Buffer(buffer), in_epsg=4326, out_epsg=dem_info['epsg'])
+    min_x, max_x, min_y, max_y = poly_bounds_in_proj(polygon.Buffer(buffer), out_epsg=dem_info['epsg'])
     output_bounds = (min_x, min_y, max_x, max_y)
     log.debug(f'Subset bounds: {output_bounds}')
 

hyp3_autorift/process.py

@@ -206,7 +206,7 @@ def process(reference: str, secondary: str, band: str = 'B08') -> Path:
         lat_limits = (bbox[1], bbox[3])
         lon_limits = (bbox[0], bbox[2])
 
-    scene_poly = geometry.polygon_from_bbox(lat_limits, lon_limits)
+    scene_poly = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     tifs = io.subset_jpl_tifs(scene_poly, target_dir=Path.cwd())
 
     geogrid_parameters = f'-d {tifs["h"]} -ssm {tifs["StableSurface"]} ' \

hyp3_autorift/vend/PRE109-PATCH-2.diff

@@ -0,0 +1,14 @@
+diff --git a/hyp3_autorift/vend/testGeogridOptical.py b/hyp3_autorift/vend/testGeogridOptical.py
+index 27793bd..eee302c 100755
+--- a/hyp3_autorift/vend/testGeogridOptical.py
++++ b/hyp3_autorift/vend/testGeogridOptical.py
+@@ -126,6 +126,9 @@ def coregisterLoadMetadata(indir_r, indir_s, urlflag):
+
+     x1a, y1a, xsize1, ysize1, x2a, y2a, xsize2, ysize2, trans = obj.coregister(indir_r, indir_s, urlflag)
+
++    xsize1 = min((xsize1, xsize2))
++    ysize1 = min((ysize1, ysize2))
++
+     if urlflag is 1:
+         DS = gdal.Open('/vsicurl/%s' % (indir_r))
+     else:

hyp3_autorift/vend/README.md

@@ -24,6 +24,12 @@ language.
 
 ## `testGeogrid_ISCE.py` and `testGeogridOptical.py`
 
+---
+*Note: A patch from autoRIFT was applied to these files to prevent out of index
+errors when coregistering optical scenes, which will be included in the next autoRIFT
+release. See `PRE109-PATCH-2.diff` for the changes applied.*
+---
+
 These modules are required for the expected workflow provided to ASF, but are
 not provided in the autoRIFT v1.0.8 release and instead resides in the "sister"
 Geogrid package (https://github.com/leiyangleon/Geogrid). Geogrid and autoRIFT

hyp3_autorift/vend/testGeogridOptical.py

@@ -126,6 +126,9 @@ def coregisterLoadMetadata(indir_r, indir_s, urlflag):
 
     x1a, y1a, xsize1, ysize1, x2a, y2a, xsize2, ysize2, trans = obj.coregister(indir_r, indir_s, urlflag)
 
+    xsize1 = min((xsize1, xsize2))
+    ysize1 = min((ysize1, ysize2))
+
     if urlflag is 1:
         DS = gdal.Open('/vsicurl/%s' % (indir_r))
     else:

tests/test_geometry.py

@@ -1,36 +1,65 @@
 import numpy as np
+from osgeo import ogr
 
 from hyp3_autorift import geometry
 
 
 def test_polygon_from_bbox():
     lat_limits = (1, 2)
     lon_limits = (3, 4)
-    assert geometry.polygon_from_bbox(lat_limits, lon_limits).ExportToWkt() \
-           == 'POLYGON ((3 1 0,4 1 0,4 2 0,3 2 0,3 1 0))'
+
+    polygon = geometry.polygon_from_bbox(lat_limits, lon_limits)
+    assert str(polygon) == 'POLYGON ((1 4,2 4,2 3,1 3,1 4))'
+    srs = polygon.GetSpatialReference()
+    assert srs.GetAttrValue('AUTHORITY', 0) == 'EPSG'
+    assert srs.GetAttrValue('AUTHORITY', 1) == '4326'
+
+    polygon = geometry.polygon_from_bbox(lat_limits, lon_limits, epsg_code=3413)
+    assert str(polygon) == 'POLYGON ((1 4,2 4,2 3,1 3,1 4))'
+    srs = polygon.GetSpatialReference()
+    assert srs.GetAttrValue('AUTHORITY', 0) == 'EPSG'
+    assert srs.GetAttrValue('AUTHORITY', 1) == '3413'
 
 
 def test_pol_bounds_in_proj():
-    polygon = geometry.polygon_from_bbox(lat_limits=(55, 56), lon_limits=(40, 41))
+    lat_limits = (55, 56)
+    lon_limits = (40, 41)
+    polygon = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     assert np.allclose(
-        geometry.poly_bounds_in_proj(polygon, in_epsg=4326, out_epsg=3413),  # NPS
+        geometry.poly_bounds_in_proj(polygon, 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))
+    lat_limits = (-58, -57)
+    lon_limits = (40, 41)
+    polygon = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     assert np.allclose(
-        geometry.poly_bounds_in_proj(polygon, in_epsg=4326, out_epsg=3031),  # SPS
+        geometry.poly_bounds_in_proj(polygon, 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))
+    lat_limits = (22, 23)
+    lon_limits = (40, 41)
+    polygon = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     assert np.allclose(
-        geometry.poly_bounds_in_proj(polygon, in_epsg=4326, out_epsg=32637),
+        geometry.poly_bounds_in_proj(polygon, out_epsg=32637),
         (602485.1663686256, 706472.0593133729, 2433164.428653589, 2544918.1043369616)
     )
 
-    polygon = geometry.polygon_from_bbox(lat_limits=(-23, -22), lon_limits=(40, 41))
+    lat_limits = (-23, -22)
+    lon_limits = (40, 41)
+    polygon = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     assert np.allclose(
-        geometry.poly_bounds_in_proj(polygon, in_epsg=4326, out_epsg=32737),
+        geometry.poly_bounds_in_proj(polygon, out_epsg=32737),
         (602485.1663686256, 706472.0593133729, 7455081.895663038, 7566835.5713464115)
     )
+
+
+def test_flip_point_coordinates():
+    point = ogr.Geometry(ogr.wkbPoint)
+    point.AddPoint_2D(0, 1)
+    assert str(geometry.flip_point_coordinates(point)) == 'POINT (1 0)'
+
+    point = ogr.Geometry(ogr.wkbPoint)
+    point.AddPoint(2, 3, 4)
+    assert str(geometry.flip_point_coordinates(point)) == 'POINT (3 2)'

tests/test_io.py

@@ -25,42 +25,62 @@ def test_download_s3_file_requester_pays(tmp_path, s3_stub):
 
 
 def test_find_jpl_dem():
-    polygon = geometry.polygon_from_bbox(lat_limits=(55, 56), lon_limits=(40, 41))
+    lat_limits = (55, 56)
+    lon_limits = (40, 41)
+    polygon = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     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))
+    lat_limits = (54, 55)
+    lon_limits = (40, 41)
+    polygon = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     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))
+    lat_limits = (54, 55)
+    lon_limits = (-40, -41)
+    polygon = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     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))
+    lat_limits = (-54, -55)
+    lon_limits = (-40, -41)
+    polygon = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     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))
+    lat_limits = (-55, -56)
+    lon_limits = (40, 41)
+    polygon = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     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))
+    lat_limits = (-56, -57)
+    lon_limits = (40, 41)
+    polygon = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     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))
+    lat_limits = (-90, -91)
+    lon_limits = (40, 41)
+    polygon = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     with pytest.raises(DemError):
         io.find_jpl_dem(polygon)
 
-    polygon = geometry.polygon_from_bbox(lat_limits=(90, 91), lon_limits=(40, 41))
+    lat_limits = (90, 91)
+    lon_limits = (40, 41)
+    polygon = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     with pytest.raises(DemError):
         io.find_jpl_dem(polygon)
 
-    polygon = geometry.polygon_from_bbox(lat_limits=(55, 56), lon_limits=(180, 181))
+    lat_limits = (55, 56)
+    lon_limits = (180, 181)
+    polygon = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     with pytest.raises(DemError):
         io.find_jpl_dem(polygon)
 
-    polygon = geometry.polygon_from_bbox(lat_limits=(55, 56), lon_limits=(-180, -181))
+    lat_limits = (55, 56)
+    lon_limits = (-180, -181)
+    polygon = geometry.polygon_from_bbox(x_limits=lat_limits, y_limits=lon_limits)
     with pytest.raises(DemError):
         io.find_jpl_dem(polygon)