Merge pull request #182 from ASFHyP3/develop

Release v0.5.2
ASFHyP3 · Mar 15, 2023 · 329d8e1 · 329d8e1
2 parents 6313bd7 + c616243
commit 329d8e1
Show file tree

Hide file tree

Showing 6 changed files with 65 additions and 26 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -7,6 +7,16 @@ 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.5.2]
+
+### Added
+* Updated tests to use a water_map/flood_map job that was created using a known input SLC
+
+### Fixed
+* Patched issue with gdalcompare.py's handling of nan values by allowing one differences between two rasters
+  that contain nan values. This patch can be remove once the upstream fix is released within GDAL (likely v3.7.0)
+* Fixed incorrect datatype being set for `flood_mask` GeoTIFFs leading to missing nodata.
+
 ## [0.5.1]
 
 ### Changed

diff --git a/src/asf_tools/flood_map.py b/src/asf_tools/flood_map.py
@@ -214,13 +214,15 @@ def make_flood_map(out_raster: Union[str, Path],  vv_raster: Union[str, Path],
     flood_depth[flood_depth < 0] = 0
 
     nodata = -1
-    floodmask_nodata = np.iinfo(np.uint8).max
     flood_depth[padding_mask] = nodata
-    flood_mask[padding_mask] = floodmask_nodata
+
+    floodmask_nodata = np.iinfo(np.uint8).max
+    flood_mask_byte = flood_mask.astype(np.uint8)
+    flood_mask_byte[padding_mask] = floodmask_nodata
 
     write_cog(str(out_raster).replace('.tif', f'_{estimator}_WaterDepth.tif'), flood_depth, transform=geotransform,
               epsg_code=epsg, dtype=gdal.GDT_Float64, nodata_value=nodata)
-    write_cog(str(out_raster).replace('.tif', f'_{estimator}_FloodMask.tif'), flood_mask, transform=geotransform,
+    write_cog(str(out_raster).replace('.tif', f'_{estimator}_FloodMask.tif'), flood_mask_byte, transform=geotransform,
               epsg_code=epsg, dtype=gdal.GDT_Byte, nodata_value=floodmask_nodata)
 
     flood_mask[known_water_mask] = False

diff --git a/src/asf_tools/hand/calculate.py b/src/asf_tools/hand/calculate.py
@@ -149,7 +149,7 @@ def calculate_hand_for_basins(out_raster:  Union[str, Path], geometries: Geometr
     """
     with rasterio.open(dem_file) as src:
         basin_mask, basin_affine_tf, basin_window = rasterio.mask.raster_geometry_mask(
-            src, geometries, all_touched=True, crop=True, pad=True, pad_width=1
+            src, geometries.geoms, all_touched=True, crop=True, pad=True, pad_width=1
         )
         basin_array = src.read(1, window=basin_window)
 

diff --git a/tests/test_flood_map.py b/tests/test_flood_map.py
@@ -9,25 +9,28 @@
 
 
 def test_get_coordinates():
-    water_raster = '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/flood_map/watermap.tif'
+    water_raster = '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/' \
+                   'S1A_IW_20230228T120437_DVR_RTC30/flood_map/watermap.tif'
     info = gdal.Info(water_raster, format='json')
 
     west, south, east, north = flood_map.get_coordinates(info)
 
-    assert west == 75390.0
-    assert south == 3265560.0
-    assert east == 357030.0
-    assert north == 3473460.0
+    assert west == 101460.0
+    assert south == 2457570.0
+    assert east == 386160.0
+    assert north == 2681970.0
 
 
 @pytest.mark.integration
 def test_get_waterbody():
-    water_raster = '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/flood_map/watermap.tif'
+    water_raster = '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/' \
+                   'S1A_IW_20230228T120437_DVR_RTC30/flood_map/watermap.tif'
     info = gdal.Info(water_raster, format='json')
 
     known_water_mask = flood_map.get_waterbody(info, threshold=30)
 
-    test_mask = '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/flood_map/known_water_mask.tif'
+    test_mask = '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/S1A_IW_20230228T120437_DVR_RTC30/' \
+                'flood_map/known_water_mask.tif'
     test_mask_array = gdal.Open(test_mask, gdal.GA_ReadOnly).ReadAsArray()
 
     assert np.all(known_water_mask == test_mask_array)
@@ -40,13 +43,15 @@ def test_logstat():
     assert np.isclose(logstd, 25.95455351947008)
 
 
+"""
 @pytest.mark.integration
 def test_estimate_flood_depths_iterative(flood_window, hand_window):
     water_height = flood_map.estimate_flood_depth(1, hand_window, flood_window, estimator='iterative',
                                                   water_level_sigma=3,
                                                   iterative_bounds=(0, 25))
     # FIXME: Basin-hopping appears to be non-deterministic. Return values vary *wildly*.
     # assert np.isclose(water_height, 7.394713346252969)
+"""
 
 
 @pytest.mark.integration
@@ -74,18 +79,22 @@ def test_estimate_flood_depths_numpy(flood_window, hand_window):
 
 @pytest.mark.integration
 def test_make_flood_map(tmp_path):
-    water_raster = '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/flood_map/watermap.tif'
-    vv_raster = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/flood_map/RTC_VV.tif'
-    hand_raster = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/flood_map/watermap_HAND.tif'
+    water_raster = '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/S1A_IW_20230228T120437_DVR_RTC30/' \
+                   'flood_map/watermap.tif'
+    vv_raster = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/S1A_IW_20230228T120437_DVR_RTC30/' \
+                'flood_map/RTC_VV.tif'
+    hand_raster = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/S1A_IW_20230228T120437_DVR_RTC30/' \
+                  'flood_map/watermap_HAND.tif'
 
     out_flood_map = tmp_path / 'flood_map.tif'
-    flood_map.make_flood_map(out_flood_map, vv_raster, water_raster, hand_raster)
-    out_flood_map = out_flood_map.parent / f'{out_flood_map.stem}_iterative_FloodDepth.tif'
+    flood_map.make_flood_map(out_flood_map, vv_raster, water_raster, hand_raster, estimator='nmad')
+    out_flood_map = out_flood_map.parent / f'{out_flood_map.stem}_nmad_FloodDepth.tif'
 
     assert out_flood_map.exists()
 
-    golden_flood_map = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/flood_map/' \
-                       'flood_map_iterative.tif'
+    golden_flood_map = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/' \
+                       'asf-tools/S1A_IW_20230228T120437_DVR_RTC30/' \
+                       'flood_map/flood_map_nmad.tif'
 
     diffs = find_diff(golden_flood_map, str(out_flood_map))
     assert diffs == 0
diff --git a/tests/test_hand.py b/tests/test_hand.py
@@ -3,17 +3,29 @@
 import pytest
 from osgeo import gdal, ogr
 from osgeo_utils.gdalcompare import find_diff
+import numpy as np
 
 from asf_tools import hand
 
+
 HAND_BASINS = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/' \
-              'asf-tools/hand/hybas_af_lev12_v1c_firstpoly.geojson'
+              'asf-tools/S1A_IW_20230228T120437_DVR_RTC30/hand/hybas_af_lev12_v1c_firstpoly.geojson'
 GOLDEN_HAND = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/' \
-              'asf-tools/hand/hybas_af_lev12_v1c_firstpoly.tif'
+              'asf-tools/S1A_IW_20230228T120437_DVR_RTC30/hand/hybas_af_lev12_v1c_firstpoly.tif'
 
 gdal.UseExceptions()
 
 
+def nodata_equal_nan(GOLDEN_HAND, out_hand):
+    ds_golden = gdal.Open(str(GOLDEN_HAND))
+    ds_out = gdal.Open(str(out_hand))
+    nodata_golden = ds_golden.GetRasterBand(1).GetNoDataValue()
+    nodata_out = ds_out.GetRasterBand(1).GetNoDataValue()
+    if nodata_golden and nodata_out:
+        return np.isnan(nodata_golden) and np.isnan(nodata_out)
+    else:
+        return False
+
 @pytest.mark.integration
 def test_make_copernicus_hand(tmp_path):
 
@@ -23,7 +35,10 @@ def test_make_copernicus_hand(tmp_path):
     assert out_hand.exists()
 
     diffs = find_diff(str(GOLDEN_HAND), str(out_hand))
-    assert diffs == 0
+    if nodata_equal_nan(str(GOLDEN_HAND), str(out_hand)):
+        assert diffs == 1
+    else:
+        assert diffs == 0
 
 
 def test_prepare_hand_vrt_no_coverage():

diff --git a/tests/test_water_map.py b/tests/test_water_map.py
@@ -38,16 +38,19 @@ def test_select_backscatter_tiles(hand_candidates):
 
 @pytest.mark.integration
 def test_make_water_map(tmp_path):
-    vv_geotif = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/water-map/20200603_VV.tif'
-    vh_geotif = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/water-map/20200603_VH.tif'
-    hand_geotif = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/water-map/20200603_HAND.tif'
+    vv_geotif = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' \
+                'S1A_IW_20230228T120437_DVR_RTC30/water_map/RTC_VV.tif'
+    vh_geotif = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' \
+                'S1A_IW_20230228T120437_DVR_RTC30/water_map/RTC_VH.tif'
+    hand_geotif = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' \
+                  'S1A_IW_20230228T120437_DVR_RTC30/water_map/HAND.tif'
 
     out_water_map = tmp_path / 'water_map.tif'
     water_map.make_water_map(out_water_map, vv_geotif, vh_geotif, hand_geotif)
 
     assert out_water_map.exists()
 
-    golden_water_map = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/water-map/' \
-                       'fuzzy-water-map.tif'
+    golden_water_map = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' \
+                       'S1A_IW_20230228T120437_DVR_RTC30/water_map/fuzzy_water_map.tif'
     diffs = find_diff(golden_water_map, str(out_water_map))
     assert diffs == 0