Repository: ASFHyP3/asf-tools (language: python)

Issue:
Update workflow
<!--
If this is a pull request for a new release, please use the release template:
   https://github.com/ASFHyP3/asf-tools/compare/main...develop?template=release.md
 
NOTE: Pull requests should only be opened for merges to protected branches (required) and any 
changes which you'd like reviewed. Do not open a pull request to update a feature or personal
branch -- simply merge with `git`.
-->


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/static-analysis.yml b/.github/workflows/static-analysis.yml
index c99dfb38..d53e0804 100644
--- a/.github/workflows/static-analysis.yml
+++ b/.github/workflows/static-analysis.yml
@@ -23,15 +23,8 @@ jobs:
               --application-import-names asf_tools ArcGIS-toolbox/ASF_Tools.pyt asf_tools
 
 
-  gitleaks:
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v2
-        with:
-          fetch-depth: 0
-
-      - name: Scan for secrets with gitleaks
-        uses: zricethezav/gitleaks-action@master
+  call-secrets-analysis-workflow:
+    uses: ASFHyP3/actions/.github/workflows/reusable-secrets-analysis.yml@main
 
 
   trufflehog:
@@ -47,11 +40,11 @@ jobs:
 
       - name: Install dependencies
         run: |
-          python -m pip install trufflehog gitdb2=="2.*"
+          python -m pip install trufflehog gitdb2=='2.*'
 
       - name: Scan for secrets with trufflehog
         run: |
           git fetch origin +refs/tags/*:refs/tags/*
           export LAST_TAG_HASH=$(git show-ref --hash -- $(git describe --abbrev=0))
-          trufflehog --regex --entropy True --since_commit "${LAST_TAG_HASH}" \
-              --exclude_paths .trufflehog.txt file://"${PWD}"
+          trufflehog --regex --entropy True --since_commit '${LAST_TAG_HASH}' \
+              --exclude_paths .trufflehog.txt file://'${PWD}'
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4f60f507..d570ee26 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,16 @@ 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.1](https://github.com/ASFHyP3/asf-tools/compare/v0.4.0...v0.4.1)
+
+### Changed
+* `asf_tools.dem` now uses the updated [2021 release of the Copernicus GLO-30 Public DEM](https://spacedata.copernicus.eu/blogs/-/blogs/copernicus-dem-2021-release-now-available)
+
+### Fixed
+* `asf_tools.hand.calculate` correctly uses [pyshed's `sGrid`](https://github.com/mdbartos/pysheds) for calculating HAND across all hydrobasins simultaneously. 
+* `asf_tools.hand.calculate` will fill NaNs within the hydrobasins in the calculated HAND array with values interpolated from their neighbor's HOND (height of nearest drainage)
+* `asf_tools.flood_map.iteartive` now produces more precise flood depth estimates by averaging water levels from a range of different initial guesses
+
 ## [0.4.0]
 
 ### Added
diff --git a/asf_tools/dem.py b/asf_tools/dem.py
index 9ee9b9ad..8e2a4fbb 100644
--- a/asf_tools/dem.py
+++ b/asf_tools/dem.py
@@ -8,7 +8,7 @@
 from asf_tools import vector
 from asf_tools.util import GDALConfigManager
 
-DEM_GEOJSON = '/vsicurl/https://asf-dem-west.s3.amazonaws.com/v2/cop30.geojson'
+DEM_GEOJSON = '/vsicurl/https://asf-dem-west.s3.amazonaws.com/v2/cop30-2021.geojson'
 
 gdal.UseExceptions()
 ogr.UseExceptions()
diff --git a/asf_tools/flood_map.py b/asf_tools/flood_map.py
index abf969b6..345e1a4b 100644
--- a/asf_tools/flood_map.py
+++ b/asf_tools/flood_map.py
@@ -69,14 +69,15 @@ def __call__(self, **kwargs):
             return tmax and tmin
 
     bounds = MyBounds()
-    x0 = [np.mean(water_levels)]
-    opt_res = optimize.basinhopping(_goal_ts, x0, niter=10000, niter_success=100, accept_test=bounds)
-    if opt_res.message[0] == 'success condition satisfied' \
-            or opt_res.message[0] == 'requested number of basinhopping iterations completed successfully':
-        best_water_level = opt_res.x[0]
-    else:
-        best_water_level = np.inf  # unstable solution.
-    return best_water_level
+    temp_wl = np.zeros(max(water_levels))
+    for i in range(1, max(water_levels)):
+        opt_res = optimize.basinhopping(_goal_ts, i, niter=10000, niter_success=100, accept_test=bounds)
+        if opt_res.message[0] == 'success condition satisfied' \
+                or opt_res.message[0] == 'requested number of basinhopping iterations completed successfully':
+            temp_wl[i] = opt_res.x[0]
+        else:
+            temp_wl[i] = np.inf  # set as inf to mark unstable solution
+    return np.nanmean(temp_wl)
 
 
 def logstat(data: np.ndarray, func: Callable = np.nanstd) -> Union[np.ndarray, float]:
diff --git a/asf_tools/hand/calculate.py b/asf_tools/hand/calculate.py
index eb18fd50..1e81d6a6 100644
--- a/asf_tools/hand/calculate.py
+++ b/asf_tools/hand/calculate.py
@@ -12,7 +12,7 @@
 import numpy as np
 import rasterio.crs
 import rasterio.mask
-from pysheds.pgrid import Grid as Pgrid
+from pysheds.sgrid import sGrid
 from shapely.geometry import GeometryCollection, shape
 
 from asf_tools.composite import write_cog
@@ -30,14 +30,31 @@ def fill_nan(array: np.ndarray) -> np.ndarray:
     kernel = astropy.convolution.Gaussian2DKernel(x_stddev=3)  # kernel x_size=8*stddev
     with warnings.catch_warnings():
         warnings.simplefilter("ignore")
-
-        array = astropy.convolution.interpolate_replace_nans(
-            array, kernel, convolve=astropy.convolution.convolve
-        )
+        while np.any(np.isnan(array)):
+            array = astropy.convolution.interpolate_replace_nans(
+                array, kernel, convolve=astropy.convolution.convolve
+            )
 
     return array
 
 
+def fill_hand(hand: np.ndarray, dem: np.ndarray):
+    """Replace NaNs in a HAND array with values interpolated from their neighbor's HOND
+
+    Replace NaNs in a HAND array with values interpolated from their neighbor's HOND (height of nearest drainage)
+    using a 2D Gaussian kernel. Here, HOND is defined as the DEM value less the HAND value. For the kernel, see:
+    https://docs.astropy.org/en/stable/convolution/#using-astropy-s-convolution-to-replace-bad-data
+    """
+    hond = dem - hand
+    hond = fill_nan(hond)
+
+    hand_mask = np.isnan(hand)
+    hand[hand_mask] = dem[hand_mask] - hond[hand_mask]
+    hand[hand < 0] = 0
+
+    return hand
+
+
 def calculate_hand(dem_array, dem_affine: rasterio.Affine, dem_crs: rasterio.crs.CRS, basin_mask,
                    acc_thresh: Optional[int] = 100):
     """Calculate the Height Above Nearest Drainage (HAND)
@@ -48,6 +65,8 @@ def calculate_hand(dem_array, dem_affine: rasterio.Affine, dem_crs: rasterio.crs
      see: https://www.hydrosheds.org/page/hydrobasins
 
      This involves:
+        * Filling pits (single-cells lower than their surrounding neighbors)
+            in the Digital Elevation Model (DEM)
         * Filling depressions (regions of cells lower than their surrounding neighbors)
             in the Digital Elevation Model (DEM)
         * Resolving un-drainable flats
@@ -56,7 +75,7 @@ def calculate_hand(dem_array, dem_affine: rasterio.Affine, dem_crs: rasterio.crs
         * Create a drainage mask using the accumulation threshold `acc_thresh`
         * Calculating HAND
 
-    In the HAND calculation, NaNs inside the basin filled using `fill_nan`
+    In the HAND calculation, NaNs inside the basin filled using `fill_hand`
 
     Args:
         dem_array: DEM to calculate HAND for
@@ -67,46 +86,51 @@ def calculate_hand(dem_array, dem_affine: rasterio.Affine, dem_crs: rasterio.crs
         acc_thresh: Accumulation threshold for determining the drainage mask.
             If `None`, the mean accumulation value is used
     """
+    nodata_fill_value = np.finfo(float).eps
+    with NamedTemporaryFile() as temp_file:
+        write_cog(temp_file.name, dem_array,
+                  transform=dem_affine.to_gdal(), epsg_code=dem_crs.to_epsg(),
+                  # Prevents PySheds from assuming using zero as the nodata value
+                  nodata_value=nodata_fill_value)
+
+        # From PySheds; see example usage: http://mattbartos.com/pysheds/
+        grid = sGrid.from_raster(str(temp_file.name))
+        dem = grid.read_raster(str(temp_file.name))
 
-    grid = Pgrid()
-    grid.add_gridded_data(dem_array, data_name='dem', affine=dem_affine, crs=dem_crs.to_dict(), mask=~basin_mask)
+    log.info('Fill pits in DEM')
+    pit_filled_dem = grid.fill_pits(dem)
 
     log.info('Filling depressions')
-    grid.fill_depressions('dem', out_name='flooded_dem')
-    if np.isnan(grid.flooded_dem).any():
-        log.debug('NaNs encountered in flooded DEM; filling.')
-        grid.flooded_dem = fill_nan(grid.flooded_dem)
+    flooded_dem = grid.fill_depressions(pit_filled_dem)
+    del pit_filled_dem
 
     log.info('Resolving flats')
-    grid.resolve_flats('flooded_dem', out_name='inflated_dem')
-    if np.isnan(grid.inflated_dem).any():
-        log.debug('NaNs encountered in inflated DEM; replacing NaNs with original DEM values')
-        grid.inflated_dem[np.isnan(grid.inflated_dem)] = dem_array[np.isnan(grid.inflated_dem)]
+    inflated_dem = grid.resolve_flats(flooded_dem)
+    del flooded_dem
 
     log.info('Obtaining flow direction')
-    grid.flowdir(data='inflated_dem', out_name='dir', apply_mask=True)
-    if np.isnan(grid.dir).any():
-        log.debug('NaNs encountered in flow direction; filling.')
-        grid.dir = fill_nan(grid.dir)
+    flow_dir = grid.flowdir(inflated_dem, apply_mask=True)
 
     log.info('Calculating flow accumulation')
-    grid.accumulation(data='dir', out_name='acc')
-    if np.isnan(grid.acc).any():
-        log.debug('NaNs encountered in accumulation; filling.')
-        grid.acc = fill_nan(grid.acc)
+    acc = grid.accumulation(flow_dir)
 
     if acc_thresh is None:
-        acc_thresh = grid.acc.mean()
+        acc_thresh = acc.mean()
 
     log.info(f'Calculating HAND using accumulation threshold of {acc_thresh}')
-    hand = grid.compute_hand('dir', 'inflated_dem', grid.acc > acc_thresh, inplace=False)
+    hand = grid.compute_hand(flow_dir, inflated_dem, acc > acc_thresh, inplace=False)
+
     if np.isnan(hand).any():
-        log.debug('NaNs encountered in HAND; filling.')
-        hand = fill_nan(hand)
+        log.info('Filling NaNs in the HAND')
+        # mask outside of basin with a not-NaN value to prevent NaN-filling outside of basin (optimization)
+        hand[basin_mask] = nodata_fill_value
+        hand = fill_hand(hand, dem_array)
 
-    # ensure non-basin is masked after fill_nan
+    # set pixels outside of basin to nodata
     hand[basin_mask] = np.nan
 
+    # TODO: also mask ocean pixels here?
+
     return hand
 
 
@@ -129,7 +153,9 @@ def calculate_hand_for_basins(out_raster:  Union[str, Path], geometries: Geometr
 
         hand = calculate_hand(basin_array, basin_affine_tf, src.crs, basin_mask)
 
-        write_cog(str(out_raster), hand, transform=basin_affine_tf.to_gdal(), epsg_code=src.crs.to_epsg())
+        write_cog(
+            out_raster, hand, transform=basin_affine_tf.to_gdal(), epsg_code=src.crs.to_epsg(), nodata_value=np.nan,
+        )
 
 
 def make_copernicus_hand(out_raster:  Union[str, Path], vector_file: Union[str, Path]):


Test patch:
diff --git a/tests/test_flood_map.py b/tests/test_flood_map.py
index cd28edbe..14359160 100644
--- a/tests/test_flood_map.py
+++ b/tests/test_flood_map.py
@@ -37,7 +37,7 @@ def test_logstat():
     arr = [10, 100, 1000, 10000, 100000]
     logstd = flood_map.logstat(arr)
 
-    assert logstd == 25.95455351947008
+    assert np.isclose(logstd, 25.95455351947008)
 
 
 @pytest.mark.integration
diff --git a/tests/test_vector.py b/tests/test_vector.py
index 37c32faf..6ced5610 100644
--- a/tests/test_vector.py
+++ b/tests/test_vector.py
@@ -40,7 +40,7 @@ def test_get_intersecting_feature_properties():
     }
     geometry = ogr.CreateGeometryFromJson(json.dumps(geojson))
     assert vector.intersecting_feature_properties(geometry, dem_tile_features, 'file_path') == [
-        '/vsicurl/https://copernicus-dem-30m.s3.amazonaws.com/'
+        '/vsicurl/https://copernicus-dem-30m.s3.amazonaws.com/2021/'
         'Copernicus_DSM_COG_10_S46_00_E169_00_DEM/Copernicus_DSM_COG_10_S46_00_E169_00_DEM.tif'
     ]
 
@@ -50,13 +50,13 @@ def test_get_intersecting_feature_properties():
     }
     geometry = ogr.CreateGeometryFromJson(json.dumps(geojson))
     assert vector.intersecting_feature_properties(geometry, dem_tile_features, 'file_path') == [
-        '/vsicurl/https://copernicus-dem-30m.s3.amazonaws.com/'
-        'Copernicus_DSM_COG_10_S46_00_E169_00_DEM/Copernicus_DSM_COG_10_S46_00_E169_00_DEM.tif',
-        '/vsicurl/https://copernicus-dem-30m.s3.amazonaws.com/'
+        '/vsicurl/https://copernicus-dem-30m.s3.amazonaws.com/2021/'
         'Copernicus_DSM_COG_10_N73_00_W122_00_DEM/Copernicus_DSM_COG_10_N73_00_W122_00_DEM.tif',
+        '/vsicurl/https://copernicus-dem-30m.s3.amazonaws.com/2021/'
+        'Copernicus_DSM_COG_10_S46_00_E169_00_DEM/Copernicus_DSM_COG_10_S46_00_E169_00_DEM.tif',
     ]
 
 
 def test_get_features():
-    for vector_file in (dem.DEM_GEOJSON, hand.prepare.HAND_GEOJSON):
-        assert len(vector.get_features(vector_file)) == 26445
+    assert len(vector.get_features(dem.DEM_GEOJSON)) == 26450
+    assert len(vector.get_features(hand.prepare.HAND_GEOJSON)) == 26445