New Water Masks

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).
 
+## [0.10.0]
+### Added
+* Support for a new water masking dataset based off of OpenStreetMaps and ESA WorldCover data.
+### Removed
+* Polygon processing functions: `split_geometry_on_antimeridian` and `get_envelope_wgs84` from `water_mask.py`.
 
 ## [0.9.3]
 ### Changed

pyproject.toml

@@ -65,6 +65,7 @@ insar_stripmap = "hyp3_isce2.insar_stripmap:main"
 [tool.pytest.ini_options]
 testpaths = ["tests"]
 script_launch_mode = "subprocess"
+markers = "integration"
 
 [tool.setuptools]
 include-package-data = true

src/hyp3_isce2/insar_tops_burst.py

@@ -136,7 +136,7 @@ def insar_tops_burst(
     if apply_water_mask:
         topsapp.run_topsapp_burst(start='computeBaselines', end='filter', config_xml=config_path)
         water_mask_path = 'water_mask.wgs84'
-        create_water_mask(str(dem_path), water_mask_path, gdal_format='ISCE')
+        create_water_mask(str(dem_path), water_mask_path)
         multilook('merged/lon.rdr.full', outname='merged/lon.rdr', alks=azimuth_looks, rlks=range_looks)
         multilook('merged/lat.rdr.full', outname='merged/lat.rdr', alks=azimuth_looks, rlks=range_looks)
         resample_to_radar_io(water_mask_path, 'merged/lat.rdr', 'merged/lon.rdr', 'merged/water_mask.rdr')

src/hyp3_isce2/metadata/templates/insar_burst/readme.md.txt.j2

@@ -280,13 +280,9 @@ process caused by invalid coherence over water bodies. The water mask will also
 in the output products. This product {{ "has" if apply_water_mask else "has not" }} 
 had the water mask applied.
 
-The water mask is generated using the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG)
-dataset (https://www.ngdc.noaa.gov/mgg/shorelines). To generate the global mask, we combined the full-resolution
-L1 and L5 datasets, and removed the L2 dataset minus the L3 dataset. The L1 dataset is the boundary between land
-and ocean, excluding Antarctica, and the L5 dataset is the boundary between Antarctic ice and ocean. The L2
-dataset is the boundary between lakes and land, and the L3 dataset is the boundary between islands and the
-lakes they are within. The GSHHG dataset was last updated in 2017, so there may be discrepancies where
-shorelines have changed.
+The water mask is generated using data from OpenStreetMaps and/or ESA WorldCover depeding on location. Areas within 
+Canada, Alaska, and Russia are primarily covered by ESA WorldCover data, while the rest of the world is covered
+by OpenStreetMaps data. Water masks were previously generated from the GSHHG dataset.
 
 *************
 # Burst InSAR Processing #

src/hyp3_isce2/water_mask.py

@@ -1,55 +1,76 @@
 """Create and apply a water body mask"""
-import json
-import subprocess
+from os import system
 from pathlib import Path
-from tempfile import TemporaryDirectory
+from typing import Optional
 
-import geopandas as gpd
+import numpy as np
 from osgeo import gdal
-from pyproj import CRS
-from shapely import geometry
-
-from hyp3_isce2.utils import GDALConfigManager
 
 gdal.UseExceptions()
 
+TILE_PATH = '/vsicurl/https://asf-dem-west.s3.amazonaws.com/WATER_MASK/TILES/'
 
-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_corners(filename, tmp_path: Optional[Path]):
+    """Get all four corners of the given image: [upper_left, bottom_left, upper_right, bottom_right].
 
-def get_envelope_wgs84(input_image: str):
-    """Get the envelope around a GeoTIFF.
     Args:
-        input_image: The path to the desired GeoTIFF, as a string.
-    Returns:
-        envelope_wgs84_gdf: The WGS84 envelope around the GeoTIFF, as a GeoDataFrame.
+        filename: The path to the input image.
+        tmp_path: An optional path to a temporary directory for temp files.
     """
-    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)
+    tmp_file = 'tmp.tif' if not tmp_path else str(tmp_path / Path('tmp.tif'))
+    ds = gdal.Warp(tmp_file, filename, dstSRS='EPSG:4326')
+    geotransform = ds.GetGeoTransform()
+    x_min = geotransform[0]
+    x_max = x_min + geotransform[1] * ds.RasterXSize
+    y_max = geotransform[3]
+    y_min = y_max + geotransform[5] * ds.RasterYSize
+    upper_left = [x_min, y_max]
+    bottom_left = [x_min, y_min]
+    upper_right = [x_max, y_max]
+    bottom_right = [x_max, y_min]
+    return [upper_left, bottom_left, upper_right, bottom_right]
+
+
+def coord_to_tile(coord: tuple[float, float]) -> str:
+    """Get the filename of the tile which encloses the inputted coordinate.
 
-    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')
+    Args:
+        coord: The (lon, lat) tuple containing the desired coordinate.
+    """
+    lat_rounded = np.floor(coord[1] / 5) * 5
+    lon_rounded = np.floor(coord[0] / 5) * 5
+    if lat_rounded >= 0:
+        lat_part = 'n' + str(int(lat_rounded)).zfill(2)
+    else:
+        lat_part = 's' + str(int(np.abs(lat_rounded))).zfill(2)
+    if lon_rounded >= 0:
+        lon_part = 'e' + str(int(lon_rounded)).zfill(3)
+    else:
+        lon_part = 'w' + str(int(np.abs(lon_rounded))).zfill(3)
+    return lat_part + lon_part + '.tif'
+
+
+def get_tiles(filename: str, tmp_path: Optional[Path]) -> None:
+    """Get the AWS vsicurl path's to the tiles necessary to cover the inputted file.
 
-    return envelope_wgs84_gdf
+    Args:
+        filename: The path to the input file.
+        tmp_path: An optional path to a temporary directory for temp files.
+    """
+    tiles = []
+    corners = get_corners(filename, tmp_path=tmp_path)
+    for corner in corners:
+        tile = TILE_PATH + coord_to_tile(corner)
+        if tile not in tiles:
+            tiles.append(tile)
+    return tiles
 
 
-def create_water_mask(input_image: str, output_image: str, gdal_format='GTiff'):
+def create_water_mask(input_image: str, output_image: str, gdal_format='ISCE', tmp_path: Optional[Path] = Path('.')):
     """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
-    https://www.soest.hawaii.edu/pwessel/gshhg/
+    The water mask is assembled from OpenStreetMaps data.
 
     Shoreline data is unbuffered and pixel values of 1 indicate land touches the pixel and 0 indicates there is no
     land in the pixel.
@@ -58,37 +79,52 @@ def create_water_mask(input_image: str, output_image: str, gdal_format='GTiff'):
         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
+        tmp_path: An optional path to a temporary directory for temp files.
     """
-    src_ds = gdal.Open(input_image)
-
-    driver_options = []
-    if gdal_format == 'GTiff':
-        driver_options = ['COMPRESS=LZW', 'TILED=YES', 'NUM_THREADS=ALL_CPUS']
-
-    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'))
 
-    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')
-        mask.to_file(temp_file, driver='ESRI Shapefile')
-        with GDALConfigManager(OGR_ENABLE_PARTIAL_REPROJECTION='YES'):
-            gdal.Rasterize(dst_ds, temp_file, allTouched=True, burnValues=[1])
+    tiles = get_tiles(input_image, tmp_path=tmp_path)
+
+    if len(tiles) < 1:
+        raise ValueError(f'No water mask tiles found for {tiles}.')
+
+    tmp_px_size_path = str(tmp_path / 'tmp_px_size.tif')
+    merged_tif_path = str(tmp_path / 'merged.tif')
+    merged_vrt_path = str(tmp_path / 'merged.vrt')
+    merged_warped_path = str(tmp_path / 'merged_warped.tif')
+    shape_path = str(tmp_path / 'tmp.shp')
+
+    pixel_size = gdal.Warp(tmp_px_size_path, input_image, dstSRS='EPSG:4326').GetGeoTransform()[1]
+
+    # This is WAY faster than using gdal_merge, because of course it is.
+    if len(tiles) > 1:
+        build_vrt_command = ' '.join(['gdalbuildvrt', merged_vrt_path] + tiles)
+        system(build_vrt_command)
+        translate_command = ' '.join(['gdal_translate', merged_vrt_path, merged_tif_path])
+        system(translate_command)
+
+    shapefile_command = ' '.join(['gdaltindex', shape_path, input_image])
+    system(shapefile_command)
+
+    warp_filename = merged_tif_path if len(tiles) > 1 else tiles[0]
+
+    gdal.Warp(
+        merged_warped_path,
+        warp_filename,
+        cutlineDSName=shape_path,
+        cropToCutline=True,
+        xRes=pixel_size,
+        yRes=pixel_size,
+        targetAlignedPixels=True,
+        dstSRS='EPSG:4326',
+        format='GTiff'
+    )
 
-    del src_ds, dst_ds
+    flip_values_command = ' '.join([
+        'gdal_calc.py',
+        '-A',
+        merged_warped_path,
+        f'--outfile={output_image}',
+        '--calc="numpy.abs((A.astype(numpy.int16) + 1) - 2)"',  # Change 1's to 0's and 0's to 1's.
+        f'--format={gdal_format}'
+    ])
+    system(flip_values_command)

tests/data/water_mask_output_info.txt

@@ -0,0 +1,31 @@
+Driver: ISCE/ISCE raster
+Files: tests/data/water_mask_output.wgs84
+       tests/data/water_mask_output.wgs84.aux.xml
+       tests/data/water_mask_output.wgs84.xml
+Size is 11, 10
+Coordinate System is:
+GEOGCRS["WGS 84",
+    DATUM["World Geodetic System 1984",
+        ELLIPSOID["WGS 84",6378137,298.257223563,
+            LENGTHUNIT["metre",1]]],
+    PRIMEM["Greenwich",0,
+        ANGLEUNIT["degree",0.0174532925199433]],
+    CS[ellipsoidal,2],
+        AXIS["geodetic latitude (Lat)",north,
+            ORDER[1],
+            ANGLEUNIT["degree",0.0174532925199433]],
+        AXIS["geodetic longitude (Lon)",east,
+            ORDER[2],
+            ANGLEUNIT["degree",0.0174532925199433]],
+    ID["EPSG",4326]]
+Data axis to CRS axis mapping: 2,1
+Origin = (-95.798036030162891,15.873381477294892)
+Pixel Size = (0.000734844751507,-0.000734844751507)
+Corner Coordinates:
+Upper Left  ( -95.7980360,  15.8733815) ( 95d47'52.93"W, 15d52'24.17"N)
+Lower Left  ( -95.7980360,  15.8660330) ( 95d47'52.93"W, 15d51'57.72"N)
+Upper Right ( -95.7899527,  15.8733815) ( 95d47'23.83"W, 15d52'24.17"N)
+Lower Right ( -95.7899527,  15.8660330) ( 95d47'23.83"W, 15d51'57.72"N)
+Center      ( -95.7939944,  15.8697073) ( 95d47'38.38"W, 15d52'10.95"N)
+Band 1 Block=11x1 Type=Byte, ColorInterp=Undefined
+  NoData Value=255