From 680c9c8f30b1688a32fc43088a6f05222208cf4f Mon Sep 17 00:00:00 2001
From: weiqi-tori <zhouwq@umich.edu>
Date: Tue, 3 Dec 2024 18:35:40 +0800
Subject: [PATCH] add HAND layer

---
 city_metrix/layers/__init__.py                |  1 +
 .../layers/height_above_nearest_drainage.py   | 45 +++++++++++++++++++
 tests/test_layers.py                          |  5 +++
 3 files changed, 51 insertions(+)
 create mode 100644 city_metrix/layers/height_above_nearest_drainage.py

diff --git a/city_metrix/layers/__init__.py b/city_metrix/layers/__init__.py
index 421b975..aa16c1f 100644
--- a/city_metrix/layers/__init__.py
+++ b/city_metrix/layers/__init__.py
@@ -26,3 +26,4 @@
 from .glad_lulc import LandCoverHabitatGlad
 from .glad_lulc import LandCoverHabitatChangeGlad
 from .cams import Cams
+from .height_above_nearest_drainage import HeightAboveNearestDrainage
diff --git a/city_metrix/layers/height_above_nearest_drainage.py b/city_metrix/layers/height_above_nearest_drainage.py
new file mode 100644
index 0000000..9f34a82
--- /dev/null
+++ b/city_metrix/layers/height_above_nearest_drainage.py
@@ -0,0 +1,45 @@
+from dask.diagnostics import ProgressBar
+import xarray as xr
+import xee
+import ee
+
+from .layer import Layer, get_utm_zone_epsg, get_image_collection
+
+
+class HeightAboveNearestDrainage(Layer):
+    """
+    Attributes:
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+        river_head: number of river head threshold cells
+    """
+
+    def __init__(self, spatial_resolution=30, river_head=1000, **kwargs):
+        super().__init__(**kwargs)
+        self.spatial_resolution = spatial_resolution
+        self.river_head = river_head
+
+    def get_data(self, bbox):
+        if self.spatial_resolution == 30 and self.river_head == 100:
+            hand = ee.ImageCollection('users/gena/global-hand/hand-100')
+            # smoothen HAND a bit, scale varies a little in the tiles
+            hand = hand.mosaic().focal_mean(0.1)
+        elif self.spatial_resolution == 30:
+            hand = ee.Image(f'users/gena/GlobalHAND/30m/hand-{self.river_head}')
+            # smoothen HAND a bit, scale varies a little in the tiles
+            hand = hand.focal_mean(0.1)
+        elif self.spatial_resolution == 90:
+            hand = ee.Image(f'users/gena/GlobalHAND/90m-global/hand-{self.river_head}')
+            # smoothen HAND a bit, scale varies a little in the tiles
+            hand = hand.focal_mean(0.1)
+
+        # MOD44W.005 Land Water Mask Derived From MODIS and SRTM
+        swbd = ee.Image('MODIS/MOD44W/MOD44W_005_2000_02_24').select('water_mask')
+        swbdMask = swbd.unmask().Not().focal_median(1)
+
+        thresh = 1
+        HANDthresh = hand.lte(thresh).focal_max(1).focal_mode(2, 'circle', 'pixels', 5).mask(swbdMask)
+        HANDthresh = HANDthresh.mask(HANDthresh)
+
+        data = get_image_collection(ee.ImageCollection(HANDthresh), bbox, self.spatial_resolution, "height above nearest drainage")
+
+        return data.b1
diff --git a/tests/test_layers.py b/tests/test_layers.py
index 7875ab9..5aa9eeb 100644
--- a/tests/test_layers.py
+++ b/tests/test_layers.py
@@ -9,6 +9,7 @@
     Era5HottestDay,
     EsaWorldCover,
     EsaWorldCoverClass,
+    HeightAboveNearestDrainage,
     HighLandSurfaceTemperature,
     ImperviousSurface,
     LandCoverGlad,
@@ -68,6 +69,10 @@ def test_esa_world_cover():
     data = EsaWorldCover(land_cover_class=land_cover_class).get_data(BBOX)
     assert np.size(data) > 0
 
+def test_height_above_nearest_drainage():
+    data = HeightAboveNearestDrainage().get_data(BBOX)
+    assert np.size(data) > 0
+
 def test_high_land_surface_temperature():
     data = HighLandSurfaceTemperature().get_data(BBOX)
     assert np.size(data) > 0