Release 6.4.1

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

@@ -54,7 +54,7 @@ jobs:
           aws-region: ${{ env.AWS_REGION }}
 
       - name: Login to Amazon ECR
-        uses: aws-actions/amazon-ecr-login@v1
+        uses: aws-actions/amazon-ecr-login@v2
 
       - name: Fetch GAMMA
         run: |

CHANGELOG.md

@@ -6,6 +6,11 @@ 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).
 
+## [6.4.1]
+
+### Fixed
+- Incorrect / blank water masks in some areas (primarily Europe), by clipping the water mask to the envelope of the product before rasterization.
+
 ## [6.4.0]
 
 ### Added

environment.yml

@@ -23,6 +23,7 @@ dependencies:
   - numpy>=1.21,<1.22
   - pillow
   - python-dateutil
+  - rtree
   # Addresses https://github.com/ASFHyP3/hyp3-gamma/issues/421
   # - hyp3lib>=1.7.0,<2
   # - lxml

hyp3_gamma/water_mask.py

@@ -1,24 +1,36 @@
 """Create and apply a water body mask"""
-import json
-import subprocess
+from pathlib import Path
 from tempfile import TemporaryDirectory
 
-import geopandas
+import geopandas as gpd
 from osgeo import gdal
+from pyproj import CRS
+from shapely import geometry
+
 
 from hyp3_gamma.util import GDALConfigManager
 
 gdal.UseExceptions()
 
 
-def split_geometry_on_antimeridian(geometry: dict):
-    geometry_as_bytes = json.dumps(geometry).encode()
-    cmd = ['ogr2ogr', '-wrapdateline', '-datelineoffset', '20', '-f', 'GeoJSON', '/vsistdout/', '/vsistdin/']
-    geojson_str = subprocess.run(cmd, input=geometry_as_bytes, stdout=subprocess.PIPE, check=True).stdout
-    return json.loads(geojson_str)['features'][0]['geometry']
+def get_envelope(input_image: str):
+    """ Get the envelope of the input_image
+
+    Args:
+        input_image: Path for the input GDAL-compatible image
+
+    Returns:
+        (envelope, epsg): The envelope and epsg code of the GeoTIFF.
+    """
+    info = gdal.Info(input_image, format='json')
+    prj = CRS.from_wkt(info["coordinateSystem"]["wkt"])
+    epsg = prj.to_epsg()
+    extent = info['wgs84Extent']
+    extent_gdf = gpd.GeoDataFrame(index=[0], geometry=[geometry.shape(extent)], crs='EPSG:4326').to_crs(epsg)
+    return extent_gdf.envelope, epsg
 
 
-def create_water_mask(input_tif: str, output_tif: str):
+def create_water_mask(input_image: str, output_image: str, gdal_format='GTiff'):
     """Create a water mask GeoTIFF with the same geometry as a given input GeoTIFF
 
     The water mask is assembled from GSHHG v2.3.7 Levels 1, 2, 3, and 5 at full resolution. To learn more, visit
@@ -28,25 +40,40 @@ def create_water_mask(input_tif: str, output_tif: str):
     land in the pixel.
 
     Args:
-        input_tif: Path for the input GeoTIFF
-        output_tif: Path for the output GeoTIFF
+        input_image: Path for the input GDAL-compatible image
+        output_image: Path for the output image
+        gdal_format: GDAL format name to create output image as
     """
-    mask_location = '/vsicurl/https://asf-dem-west.s3.amazonaws.com/WATER_MASK/GSHHG/hyp3_water_mask_20220912.shp'
+    src_ds = gdal.Open(input_image)
 
-    src_ds = gdal.Open(input_tif)
+    driver_options = []
+    if gdal_format == 'GTiff':
+        driver_options = ['COMPRESS=LZW', 'TILED=YES', 'NUM_THREADS=ALL_CPUS']
 
-    dst_ds = gdal.GetDriverByName('GTiff').Create(output_tif, src_ds.RasterXSize, src_ds.RasterYSize, 1, gdal.GDT_Byte)
+    dst_ds = gdal.GetDriverByName(gdal_format).Create(
+        output_image,
+        src_ds.RasterXSize,
+        src_ds.RasterYSize,
+        1,
+        gdal.GDT_Byte,
+        driver_options,
+    )
     dst_ds.SetGeoTransform(src_ds.GetGeoTransform())
     dst_ds.SetProjection(src_ds.GetProjection())
     dst_ds.SetMetadataItem('AREA_OR_POINT', src_ds.GetMetadataItem('AREA_OR_POINT'))
 
-    extent = gdal.Info(input_tif, format='json')['wgs84Extent']
-    extent = split_geometry_on_antimeridian(extent)
+    envelope, epsg = get_envelope(input_image)
+
+    mask_location = '/vsicurl/https://asf-dem-west.s3.amazonaws.com/WATER_MASK/GSHHG/hyp3_water_mask_20220912.shp'
+
+    mask = gpd.read_file(mask_location, mask=envelope).to_crs(epsg)
+
+    mask = gpd.clip(mask, envelope)
 
-    mask = geopandas.read_file(mask_location, mask=extent)
-    with TemporaryDirectory() as temp_shapefile:
-        mask.to_file(temp_shapefile, driver='ESRI Shapefile')
+    with TemporaryDirectory() as temp_dir:
+        temp_file = str(Path(temp_dir) / 'mask.shp')
+        mask.to_file(temp_file, driver='ESRI Shapefile')
         with GDALConfigManager(OGR_ENABLE_PARTIAL_REPROJECTION='YES'):
-            gdal.Rasterize(dst_ds, temp_shapefile, allTouched=True, burnValues=[1])
+            gdal.Rasterize(dst_ds, temp_file, allTouched=True, burnValues=[1])
 
     del src_ds, dst_ds

setup.py

@@ -41,6 +41,7 @@
         'numpy>=1.21,<1.22',
         'pillow',
         'python-dateutil',
+        'rtree'
     ],
 
     extras_require={

tests/test_water_mask.py

@@ -6,52 +6,6 @@
 gdal.UseExceptions()
 
 
-def test_split_geometry_on_antimeridian():
-    geometry = {
-        'type': 'Polygon',
-        'coordinates': [[
-            [170, 50],
-            [175, 55],
-            [-170, 55],
-            [-175, 50],
-            [170, 50],
-        ]],
-    }
-    result = water_mask.split_geometry_on_antimeridian(geometry)
-    assert result == {
-        'type': 'MultiPolygon',
-        'coordinates': [
-            [[
-                [175.0, 55.0],
-                [180.0, 55.0],
-                [180.0, 50.0],
-                [170.0, 50.0],
-                [175.0, 55.0],
-            ]],
-            [[
-                [-170.0, 55.0],
-                [-175.0, 50.0],
-                [-180.0, 50.0],
-                [-180.0, 55.0],
-                [-170.0, 55.0],
-            ]],
-        ],
-    }
-
-    geometry = {
-        'type': 'Polygon',
-        'coordinates': [[
-            [150, 50],
-            [155, 55],
-            [-150, 55],
-            [-155, 50],
-            [150, 50],
-        ]],
-    }
-    result = water_mask.split_geometry_on_antimeridian(geometry)
-    assert result == geometry
-
-
 def test_create_water_mask_with_no_water(tmp_path, test_data_dir):
     input_tif = str(test_data_dir / 'test_geotiff.tif')
     output_tif = str(tmp_path / 'water_mask.tif')