Merge pull request #164 from ASFHyP3/develop

Release 0.9.1

CHANGELOG.md

@@ -6,8 +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).
 
-## [0.9.0]
+## [0.9.1]
+### Fixed
+* Water Masking now pulls the global shapefile and clips it, rather than pulling a partitioned parque, due to issues around the partition boundaries.
 
+## [0.9.0]
 ### Changed
 * Upgraded `hyp3lib` dependency to version `2.x.x`.
 * As of [HyP3-lib v2.0.0](https://github.com/ASFHyP3/hyp3-lib/releases/tag/v2.0.0), the [Copernicus Data Space Ecosystem (CDSE)](https://dataspace.copernicus.eu/) will now be used for downloading Sentinel-1 orbit files from ESA.

environment.yml

@@ -15,9 +15,6 @@ dependencies:
   - geopandas
   - gdal
   - opencv
-  - s3fs
-  - pyarrow
-  - hyp3lib>=2,<3
   # For packaging, and testing
   - flake8
   - flake8-import-order
@@ -30,3 +27,5 @@ dependencies:
   - pytest-console-scripts
   - pytest-cov
   - pytest-mock
+  # For running
+  - hyp3lib>=2,<3

pyproject.toml

@@ -30,8 +30,6 @@ dependencies = [
   "asf_search>=6.4.0",
   "geopandas",
   "gdal",
-  "s3fs",
-  "pyarrow",
   "hyp3lib>=2,<3",
   # Conda-forge only dependencies are listed below
   # "opencv",

src/hyp3_isce2/insar_tops_burst.py

@@ -19,6 +19,7 @@
 from isceobj.TopsProc.runMergeBursts import multilook
 from lxml import etree
 from osgeo import gdal
+from pyproj import CRS
 
 import hyp3_isce2
 import hyp3_isce2.metadata.util
@@ -397,6 +398,43 @@ def get_pixel_size(looks: str) -> float:
     return {'20x4': 80.0, '10x2': 40.0, '5x1': 20.0}[looks]
 
 
+def convert_raster_from_isce2_gdal(input_image, ref_image, output_image):
+    """Convert the water mask in WGS84 to be the same projection and extent of the output product.
+
+    Args:
+        input_image: dem file name
+        ref_image: output geotiff file name
+        output_image: water mask file name
+    """
+
+    ref_ds = gdal.Open(ref_image)
+
+    gt = ref_ds.GetGeoTransform()
+
+    pixel_size = gt[1]
+
+    minx = gt[0]
+    maxx = gt[0] + gt[1] * ref_ds.RasterXSize
+    maxy = gt[3]
+    miny = gt[3] + gt[5] * ref_ds.RasterYSize
+
+    crs = ref_ds.GetSpatialRef()
+    epsg = CRS.from_wkt(crs.ExportToWkt()).to_epsg()
+
+    del ref_ds
+
+    gdal.Warp(
+        output_image,
+        input_image,
+        dstSRS=f'epsg:{epsg}',
+        creationOptions=['TILED=YES', 'COMPRESS=LZW', 'NUM_THREADS=ALL_CPUS'],
+        outputBounds=[minx, miny, maxx, maxy],
+        xRes=pixel_size,
+        yRes=pixel_size,
+        targetAlignedPixels=True
+        )
+
+
 def main():
     """HyP3 entrypoint for the burst TOPS workflow"""
     parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.ArgumentDefaultsHelpFormatter)
@@ -457,22 +495,20 @@ def main():
     translate_outputs(isce_output_dir, product_name, pixel_size=pixel_size)
 
     unwrapped_phase = f'{product_name}/{product_name}_unw_phase.tif'
-    wrapped_phase = f'{product_name}/{product_name}_wrapped_phase.tif'
     water_mask = f'{product_name}/{product_name}_water_mask.tif'
-    create_water_mask(wrapped_phase, water_mask)
 
     if apply_water_mask:
-        for geotiff in [wrapped_phase, unwrapped_phase]:
-            cmd = (
-                'gdal_calc.py '
-                f'--outfile {geotiff} '
-                f'-A {geotiff} -B {water_mask} '
-                '--calc A*B '
-                '--overwrite '
-                '--NoDataValue 0 '
-                '--creation-option TILED=YES --creation-option COMPRESS=LZW --creation-option NUM_THREADS=ALL_CPUS'
-            )
-            subprocess.run(cmd.split(' '), check=True)
+        convert_raster_from_isce2_gdal('water_mask.wgs84', unwrapped_phase, water_mask)
+        cmd = (
+            'gdal_calc.py '
+            f'--outfile {unwrapped_phase} '
+            f'-A {unwrapped_phase} -B {water_mask} '
+            '--calc A*B '
+            '--overwrite '
+            '--NoDataValue 0 '
+            '--creation-option TILED=YES --creation-option COMPRESS=LZW --creation-option NUM_THREADS=ALL_CPUS'
+        )
+        subprocess.run(cmd.split(' '), check=True)
 
     make_browse_image(unwrapped_phase, f'{product_name}/{product_name}_unw_phase.png')
 

src/hyp3_isce2/water_mask.py

@@ -5,68 +5,44 @@
 from tempfile import TemporaryDirectory
 
 import geopandas as gpd
-import numpy as np
-import pandas as pd
-import s3fs
-import shapely
 from osgeo import gdal
-from shapely import geometry, wkt
+from pyproj import CRS
+from shapely import geometry
 
 from hyp3_isce2.utils import GDALConfigManager
 
 gdal.UseExceptions()
 
 
-def get_geo_partition(coordinate, partition_size=90):
-    """Get the geo-partition for a given coordinate (i.e., the lat/lon box it falls in)
-
-    Args:
-        coordinate: A coordinate tuple (lon, lat)
-        partition_size: The partition size (in degrees) to use for the geo-partition
-            using a value of 90 will result in geo-partitions of 90x90 degrees
-
-    Returns:
-        A string representing the geo-partition for the given coordinate and partition size
-    """
-    x, y = coordinate
-    x_rounded = int(np.floor(x / partition_size)) * partition_size
-    y_rounded = int(np.floor(y / partition_size)) * partition_size
-    x_fill = str(x_rounded).zfill(4)
-    y_fill = str(y_rounded).zfill(4)
-    partition = f'{y_fill}_{x_fill}'
-    return partition
-
-
 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_water_mask_gdf(extent: geometry.Polygon) -> gpd.GeoDataFrame:
-    """Get a GeoDataFrame of the water mask for a given extent
-
+def get_envelope_wgs84(input_image: str):
+    """Get the envelope around a GeoTIFF.
     Args:
-        extent: The extent to get the water mask for
-
+        input_image: The path to the desired GeoTIFF, as a string.
     Returns:
-        GeoDataFrame of the water mask for the given extent
+        envelope_wgs84_gdf: The WGS84 envelope around the GeoTIFF, as a GeoDataFrame.
     """
-    mask_location = 'asf-dem-west/WATER_MASK/GSHHG/hyp3_water_mask_20220912'
-    corrected_extent = split_geometry_on_antimeridian(json.loads(shapely.to_geojson(extent)))
-
-    filters = list(set([('lat_lon', '=', get_geo_partition(coord)) for coord in extent.exterior.coords]))
-    s3_fs = s3fs.S3FileSystem(anon=True, default_block_size=5 * (2**20))
+    info = gdal.Info(input_image, format='json')
+    prj = CRS.from_wkt(info["coordinateSystem"]["wkt"])
+    epsg = prj.to_epsg()
+    extent = info['wgs84Extent']
+    poly = geometry.shape(extent)
+    poly_gdf = gpd.GeoDataFrame(index=[0], geometry=[poly], crs='EPSG:4326')
+    envelope_gdf = poly_gdf.to_crs(epsg).envelope.to_crs(4326)
+    envelope_poly = envelope_gdf.geometry[0]
+    envelope = geometry.mapping(envelope_poly)
 
-    # TODO the conversion from pd -> gpd can be removed when gpd adds the filter param for read_parquet
-    df = pd.read_parquet(mask_location, filesystem=s3_fs, filters=filters)
-    df['geometry'] = df['geometry'].apply(wkt.loads)
-    df['lat_lon'] = df['lat_lon'].astype(str)
-    gdf = gpd.GeoDataFrame(df, crs='EPSG:4326')
+    correct_extent = split_geometry_on_antimeridian(envelope)
+    envelope_wgs84 = geometry.shape(correct_extent)
+    envelope_wgs84_gdf = gpd.GeoDataFrame(index=[0], geometry=[envelope_wgs84], crs='EPSG:4326')
 
-    mask = gpd.clip(gdf, geometry.shape(corrected_extent))
-    return mask
+    return envelope_wgs84_gdf
 
 
 def create_water_mask(input_image: str, output_image: str, gdal_format='GTiff'):
@@ -79,7 +55,7 @@ def create_water_mask(input_image: str, output_image: str, gdal_format='GTiff'):
     land in the pixel.
 
     Args:
-        input_imge: Path for the input GDAL-compatible image
+        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
     """
@@ -101,8 +77,13 @@ def create_water_mask(input_image: str, output_image: str, gdal_format='GTiff'):
     dst_ds.SetProjection(src_ds.GetProjection())
     dst_ds.SetMetadataItem('AREA_OR_POINT', src_ds.GetMetadataItem('AREA_OR_POINT'))
 
-    extent = gdal.Info(input_image, format='json')['wgs84Extent']
-    mask = get_water_mask_gdf(geometry.shape(extent))
+    envelope_wgs84_gdf = get_envelope_wgs84(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_wgs84_gdf)
+
+    mask = mask.clip(envelope_wgs84_gdf)
 
     with TemporaryDirectory() as temp_dir:
         temp_file = str(Path(temp_dir) / 'mask.shp')

tests/test_water_mask.py

@@ -1,36 +1,29 @@
 import numpy as np
-from osgeo import gdal
+from osgeo import gdal, osr
 
 from hyp3_isce2 import water_mask
 
 gdal.UseExceptions()
 
 
-def test_get_geo_partition():
-    coordinates = [
-        (0, 0),
-        (0, 90),
-        (0, -45),
-        (-45, 0),
-        (45, 90),
-    ]
-    expected_90partition = [
-        '0000_0000',
-        '0090_0000',
-        '-090_0000',
-        '0000_-090',
-        '0090_0000',
-    ]
-    expected_45partition = [
-        '0000_0000',
-        '0090_0000',
-        '-045_0000',
-        '0000_-045',
-        '0090_0045',
-    ]
-    for coord, expected90, expected45 in zip(coordinates, expected_90partition, expected_45partition):
-        assert water_mask.get_geo_partition(coord, partition_size=90) == expected90
-        assert water_mask.get_geo_partition(coord, partition_size=45) == expected45
+def get_envelope_with_args(out_path, xres, yres, bounds):
+    xmin, ymin, xmax, ymax = bounds
+    xsize = abs(int((xmax - xmin) / xres))
+    ysize = abs(int((ymax - ymin) / yres))
+
+    out_img_path = str(out_path)
+    driver = gdal.GetDriverByName('GTiff')
+    spatref = osr.SpatialReference()
+    spatref.ImportFromEPSG(4326)
+    wkt = spatref.ExportToWkt()
+    ds = driver.Create(out_img_path, xsize, ysize, options=['COMPRESS=LZW', 'TILED=YES'])
+    ds.SetProjection(wkt)
+    ds.SetGeoTransform([xmin, xres, 0, ymin, 0, yres])
+    ds.GetRasterBand(1).Fill(0)
+    ds.FlushCache()
+    ds = None
+
+    return water_mask.get_envelope_wgs84(out_img_path)
 
 
 def test_split_geometry_on_antimeridian():
@@ -79,6 +72,15 @@ def test_split_geometry_on_antimeridian():
     assert result == geometry
 
 
+def test_get_envelope_wgs84(tmp_path):
+    out_path_1 = str(tmp_path / 'envelope1.tif')
+    out_path_2 = str(tmp_path / 'envelope2.tif')
+    envelope1 = get_envelope_with_args(out_path_1, 0.01, 0.01, [0, 40, 10, 50])
+    envelope2 = get_envelope_with_args(out_path_2, 0.01, 0.01, [-177, 40, 178, 50])
+    assert np.all(envelope1.bounds.values == np.asarray([[0.0, 40.0, 10.0, 50.0]]))
+    assert np.all(envelope2.bounds.values == np.asarray([[-180.0, 40.0, 180.0, 50.0]]))
+
+
 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')