diff --git a/.github/workflows/test-and-build.yml b/.github/workflows/test-and-build.yml index 8953e299..2fc867d4 100644 --- a/.github/workflows/test-and-build.yml +++ b/.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: | diff --git a/CHANGELOG.md b/CHANGELOG.md index 4939c836..0b672cab 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/environment.yml b/environment.yml index fc334c40..f79a6aa2 100644 --- a/environment.yml +++ b/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 diff --git a/hyp3_gamma/water_mask.py b/hyp3_gamma/water_mask.py index e29438bf..0ee493da 100755 --- a/hyp3_gamma/water_mask.py +++ b/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 diff --git a/setup.py b/setup.py index 9311d076..8ead1fa0 100644 --- a/setup.py +++ b/setup.py @@ -41,6 +41,7 @@ 'numpy>=1.21,<1.22', 'pillow', 'python-dateutil', + 'rtree' ], extras_require={ diff --git a/tests/test_water_mask.py b/tests/test_water_mask.py index c4587c8a..79850b3b 100644 --- a/tests/test_water_mask.py +++ b/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')