Merge pull request #70 from Open-ET/water-mask

Moved Tcorr water masking code to a dedicated function

openet/ssebop/image.py

@@ -475,6 +475,32 @@ def quality(self):
         """Set quality to 1 for all active pixels (for now)"""
         return self.mask.rename(['quality']).set(self._properties)
 
+    @lazy_property
+    def tcorr_not_water_mask(self):
+        """Mask of pixels that have a high confidence of not being water
+
+        The purpose for this mask is to ensure that water pixels are not used in
+            the Tcorr FANO calculation.
+
+        Output image will be 1 for pixels that are not-water and 0 otherwise
+
+        NDWI in landsat.py is defined as "green - swir1", which is "flipped"
+            compared to NDVI and other NDWI calculations,
+            so water will be positive and land will be negative
+
+        """
+        ndwi_threshold = -0.15
+
+        # TODO: Check if .multiply() is the same as .And() here
+        #   The .And() seems more readable
+        not_water_mask = (
+            ee.Image(self.ndwi).lt(ndwi_threshold)
+            .multiply(self.qa_water_mask.eq(0))
+            # .And(self.qa_water_mask.eq(0))
+        )
+
+        return not_water_mask.rename(['tcorr_not_water']).set(self._properties).uint8()
+
     @lazy_property
     def time(self):
         """Return an image of the 0 UTC time (in milliseconds)"""
@@ -817,7 +843,6 @@ def tcorr_FANO(self):
         coarse_transform = [1000, 0, 15, 0, -1000, 15]
         coarse_transform100 = [100000, 0, 15, 0, -100000, 15]
         dt_coeff = 0.125
-        ndwi_threshold = -0.15
         high_ndvi_threshold = 0.9
         water_pct = 50
         # max pixels argument for .reduceResolution()
@@ -827,24 +852,25 @@ def tcorr_FANO(self):
         ndvi = ee.Image(self.ndvi).clamp(-1.0, 1.0)
         tmax = ee.Image(self.tmax)
         dt = ee.Image(self.dt)
-        ndwi = ee.Image(self.ndwi)
-        qa_watermask = ee.Image(self.qa_water_mask)
 
-        # setting NDVI to negative values where Landsat QA Pixel detects water.
+        # Setting NDVI to negative values where Landsat QA Pixel detects water.
+        # TODO: We may want to switch "qa_watermask" to "not_water_mask.eq(0)"
+        qa_watermask = ee.Image(self.qa_water_mask)
         ndvi = ndvi.where(qa_watermask.eq(1).And(ndvi.gt(0)), ndvi.multiply(-1))
 
-        watermask = ndwi.lt(ndwi_threshold)
-        # combining NDWI mask with QA Pixel watermask.
-        watermask = watermask.multiply(qa_watermask.eq(0))
-        # returns qa_watermask layer masked by combined watermask to get a count of valid pixels
-        watermask_for_coarse = qa_watermask.updateMask(watermask)
+        # Mask with not_water pixels set to 1 and water pixels set to 0
+        not_water_mask = self.tcorr_not_water_mask
 
+        # Count not-water pixels and the total number of pixels
+        # TODO: Rename "watermask_coarse_count" here to "not_water_pixels_count"
         watermask_coarse_count = (
-            watermask_for_coarse
+            self.qa_water_mask.updateMask(not_water_mask)
             .reduceResolution(ee.Reducer.count(), False, m_pixels)
             .reproject(self.crs, coarse_transform)
             .updateMask(1).select([0], ['count'])
         )
+
+        # TODO: Maybe chance ndvi to self.qa_water_mask?
         total_pixels_count = (
             ndvi
             .reduceResolution(ee.Reducer.count(), False, m_pixels)
@@ -865,13 +891,13 @@ def tcorr_FANO(self):
 
         ndvi_avg_masked = (
             ndvi
-            .updateMask(watermask)
+            .updateMask(not_water_mask)
             .reduceResolution(ee.Reducer.mean(), False, m_pixels)
             .reproject(self.crs, coarse_transform)
         )
         ndvi_avg_masked100 = (
             ndvi
-            .updateMask(watermask)
+            .updateMask(not_water_mask)
             .reduceResolution(ee.Reducer.mean(), True, m_pixels)
             .reproject(self.crs, coarse_transform100)
         )
@@ -883,13 +909,13 @@ def tcorr_FANO(self):
         )
         lst_avg_masked = (
             lst
-            .updateMask(watermask)
+            .updateMask(not_water_mask)
             .reduceResolution(ee.Reducer.mean(), False, m_pixels)
             .reproject(self.crs, coarse_transform)
         )
         lst_avg_masked100 = (
             lst
-            .updateMask(watermask)
+            .updateMask(not_water_mask)
             .reduceResolution(ee.Reducer.mean(), True, m_pixels)
             .reproject(self.crs, coarse_transform100)
         )

openet/ssebop/tests/test_c_image.py

@@ -47,16 +47,18 @@
 #     #     })
 
 
-def default_image(lst=305, ndvi=0.8, qa_water=0):
+def default_image(lst=305, ndvi=0.8, ndwi=-0.5, qa_water=0):
     # First construct a fake 'prepped' input image
     mask_img = ee.Image(f'{COLL_ID}/{SCENE_ID}').select(['SR_B3']).multiply(0)
-    return ee.Image([mask_img.add(lst), mask_img.add(ndvi), mask_img.add(qa_water)]) \
-        .rename(['lst', 'ndvi', 'qa_water']) \
+    return (
+        ee.Image([mask_img.add(lst), mask_img.add(ndvi), mask_img.add(ndwi), mask_img.add(qa_water)])
+        .rename(['lst', 'ndvi', 'ndwi', 'qa_water'])
         .set({
             'system:index': SCENE_ID,
             'system:time_start': SCENE_TIME,
             'system:id': f'{COLL_ID}/{SCENE_ID}',
         })
+    )
     # return ee.Image.constant([lst, ndvi]).rename(['lst', 'ndvi']) \
     #     .set({
     #         'system:index': SCENE_ID,
@@ -68,7 +70,10 @@ def default_image(lst=305, ndvi=0.8, qa_water=0):
 # Setting et_reference_source and et_reference_band on the default image to
 #   simplify testing but these do not have defaults in the Image class init
 def default_image_args(
-        lst=305, ndvi=0.85,
+        lst=305,
+        ndvi=0.85,
+        ndwi=-0.5,
+        qa_water=0,
         # et_reference_source='IDAHO_EPSCOR/GRIDMET',
         et_reference_source=9.5730,
         et_reference_band='etr',
@@ -87,7 +92,7 @@ def default_image_args(
         tcorr_resample='nearest',
         ):
     return {
-        'image': default_image(lst=lst, ndvi=ndvi),
+        'image': default_image(lst=lst, ndvi=ndvi, ndwi=ndwi, qa_water=qa_water),
         'et_reference_source': et_reference_source,
         'et_reference_band': et_reference_band,
         'et_reference_factor': et_reference_factor,
@@ -107,7 +112,10 @@ def default_image_args(
 
 
 def default_image_obj(
-        lst=305, ndvi=0.85,
+        lst=305,
+        ndvi=0.85,
+        ndwi=-0.5,
+        qa_water=0,
         # et_reference_source='IDAHO_EPSCOR/GRIDMET',
         et_reference_source=9.5730,
         et_reference_band='etr',
@@ -126,7 +134,10 @@ def default_image_obj(
         tcorr_resample='nearest',
         ):
     return ssebop.Image(**default_image_args(
-        lst=lst, ndvi=ndvi,
+        lst=lst,
+        ndvi=ndvi,
+        ndwi=ndwi,
+        qa_water=qa_water,
         et_reference_source=et_reference_source,
         et_reference_band=et_reference_band,
         et_reference_factor=et_reference_factor,
@@ -221,6 +232,94 @@ def test_Image_lst_properties():
     assert output['properties']['image_id'] == f'{COLL_ID}/{SCENE_ID}'
 
 
+def test_Image_mask_properties():
+    """Test if properties are set on the time image"""
+    output = utils.getinfo(default_image_obj().mask)
+    assert output['bands'][0]['id'] == 'mask'
+    assert output['properties']['system:index'] == SCENE_ID
+    assert output['properties']['system:time_start'] == SCENE_TIME
+    assert output['properties']['image_id'] == f'{COLL_ID}/{SCENE_ID}'
+
+
+def test_Image_mask_values():
+    output_img = default_image_obj(
+        ndvi=0.5, lst=308, dt_source=10, elev_source=50,
+        tcorr_source=0.98, tmax_source=310).mask
+    output = utils.point_image_value(output_img, TEST_POINT)
+    assert output['mask'] == 1
+
+
+def test_Image_time_properties():
+    """Test if properties are set on the time image"""
+    output = utils.getinfo(default_image_obj().time)
+    assert output['bands'][0]['id'] == 'time'
+    assert output['properties']['system:index'] == SCENE_ID
+    assert output['properties']['system:time_start'] == SCENE_TIME
+    assert output['properties']['image_id'] == f'{COLL_ID}/{SCENE_ID}'
+
+
+def test_Image_time_values():
+    # The time band should be the 0 UTC datetime, not the system:time_start
+    # The time image is currently being built from the et_fraction image, so all
+    #   the ancillary values must be set.
+    output_img = default_image_obj(
+        ndvi=0.5, lst=308, dt_source=10, elev_source=50,
+        tcorr_source=0.98, tmax_source=310).time
+    output = utils.point_image_value(output_img, TEST_POINT)
+    assert output['time'] == utils.millis(SCENE_DT)
+
+
+@pytest.mark.parametrize(
+    'qa_water, expected',
+    [
+        [0, 0],
+        [1, 1],
+    ]
+)
+def test_Image_qa_water_mask_values(qa_water, expected):
+    """Test QA water mask"""
+    output_img = default_image_obj(qa_water=qa_water)
+    mask_img = output_img.qa_water_mask
+    output = utils.point_image_value(mask_img, SCENE_POINT)
+    assert output['qa_water'] == expected
+
+
+@pytest.mark.parametrize(
+    'ndvi, ndwi, qa_water, expected',
+    [
+        # Not water (Land) Tests
+        # Start with all not water conditions
+        [0.85, -0.5, 0, 1],
+        # NDVI is not currently checked in the Tcorr water/not-water mask function
+        #   so setting NDVI negative will not change mask here
+        #   but NDVI in Tcorr calculation may/will be adjusted
+        [-0.5, -0.5, 0, 1],
+
+        # Water Tests
+        # Any indication of water will set mask to 0
+        # So either qa_water being True or NDWI being > threshold
+        [-0.5, 0.5, 1, 0],
+        [-0.5, -0.5, 1, 0],
+        [-0.5, 0.5, 0, 0],
+        # NDVI is not considered so changing NDVI to positive here will not change mask value
+        [0.5, 0.5, 1, 0],
+        [0.5, -0.5, 1, 0],
+        [0.5, 0.5, 0, 0],
+
+        # Check right around NDWI threshold
+        [0.85, -0.149, 0, 0],
+        [0.85, -0.150, 0, 0],
+        [0.85, -0.151, 0, 1],
+    ]
+)
+def test_Image_tcorr_not_water_mask_values(ndvi, ndwi, qa_water, expected):
+    """Test water mask values"""
+    output_img = default_image_obj(ndvi=ndvi, ndwi=ndwi, qa_water=qa_water)
+    mask_img = output_img.tcorr_not_water_mask
+    output = utils.point_image_value(mask_img, SCENE_POINT)
+    assert output['tcorr_not_water'] == expected
+
+
 @pytest.mark.parametrize(
     'elev_source',
     [None],
@@ -234,6 +333,11 @@ def test_Image_elev_source_exception(elev_source):
         utils.getinfo(default_image_obj(elev_source=elev_source).elev)
 
 
+def test_Image_elev_band_name():
+    output = utils.getinfo(default_image_obj().elev)['bands'][0]['id']
+    assert output == 'elev'
+
+
 @pytest.mark.parametrize(
     'elev_source, xy, expected',
     [
@@ -244,7 +348,9 @@ def test_Image_elev_source_exception(elev_source):
         ['projects/usgs-ssebop/srtm_1km', [-106.03249, 37.17777], 2369.0],
         ['projects/earthengine-legacy/assets/projects/usgs-ssebop/srtm_1km',
          [-106.03249, 37.17777], 2369.0],
-        ['USGS/NED', [-106.03249, 37.17777], 2364.351],
+        ['USGS/3DEP/10m', [-106.03249, 37.17777], 2364.1377],
+        # # NED collection is deprecated
+        # ['USGS/NED', [-106.03249, 37.17777], 2364.351],
     ]
 )
 def test_Image_elev_source(elev_source, xy, expected, tol=0.001):
@@ -254,11 +360,6 @@ def test_Image_elev_source(elev_source, xy, expected, tol=0.001):
     assert abs(output['elev'] - expected) <= tol
 
 
-def test_Image_elev_band_name():
-    output = utils.getinfo(default_image_obj().elev)['bands'][0]['id']
-    assert output == 'elev'
-
-
 @pytest.mark.parametrize(
     'dt_source, doy, xy, expected',
     [
@@ -878,56 +979,34 @@ def test_Image_from_landsat_c2_sr_et():
     assert output['properties']['system:index'] == image_id.split('/')[-1]
 
 
-def test_Image_mask_properties():
-    """Test if properties are set on the time image"""
-    output = utils.getinfo(default_image_obj().mask)
-    assert output['bands'][0]['id'] == 'mask'
-    assert output['properties']['system:index'] == SCENE_ID
-    assert output['properties']['system:time_start'] == SCENE_TIME
-    assert output['properties']['image_id'] == f'{COLL_ID}/{SCENE_ID}'
-
-
-def test_Image_mask_values():
-    output_img = default_image_obj(
-        ndvi=0.5, lst=308, dt_source=10, elev_source=50,
-        tcorr_source=0.98, tmax_source=310).mask
-    output = utils.point_image_value(output_img, TEST_POINT)
-    assert output['mask'] == 1
-
-
-def test_Image_time_properties():
-    """Test if properties are set on the time image"""
-    output = utils.getinfo(default_image_obj().time)
-    assert output['bands'][0]['id'] == 'time'
-    assert output['properties']['system:index'] == SCENE_ID
-    assert output['properties']['system:time_start'] == SCENE_TIME
-    assert output['properties']['image_id'] == f'{COLL_ID}/{SCENE_ID}'
-
-
 @pytest.mark.parametrize(
     'image_id, xy, expected',
     [
         ['LANDSAT/LC08/C02/T1_L2/LC08_044033_20170716', (-122.15, 38.6153), 0],
         ['LANDSAT/LC08/C02/T1_L2/LC08_044033_20170716', (-122.2571, 38.6292), 1],
     ]
 )
-def test_Image_qa_water_mask(image_id, xy, expected):
+def test_Image_qa_water_mask_image_values(image_id, xy, expected):
     """Test if qa pixel waterband exists"""
     output_img = ssebop.Image.from_image_id(image_id)
     mask_img = output_img.qa_water_mask
     output = utils.point_image_value(mask_img, xy)
     assert output['qa_water'] == expected
 
 
-def test_Image_time_values():
-    # The time band should be the 0 UTC datetime, not the system:time_start
-    # The time image is currently being built from the et_fraction image, so all
-    #   the ancillary values must be set.
-    output_img = default_image_obj(
-        ndvi=0.5, lst=308, dt_source=10, elev_source=50,
-        tcorr_source=0.98, tmax_source=310).time
-    output = utils.point_image_value(output_img, TEST_POINT)
-    assert output['time'] == utils.millis(SCENE_DT)
+@pytest.mark.parametrize(
+    'image_id, xy, expected',
+    [
+        ['LANDSAT/LC08/C02/T1_L2/LC08_044033_20170716', (-122.15, 38.6153), 1],
+        ['LANDSAT/LC08/C02/T1_L2/LC08_044033_20170716', (-122.2571, 38.6292), 0],
+    ]
+)
+def test_Image_tcorr_not_water_mask_image_values(image_id, xy, expected):
+    """Test water mask values"""
+    output_img = ssebop.Image.from_image_id(image_id)
+    mask_img = output_img.tcorr_not_water_mask
+    output = utils.point_image_value(mask_img, xy)
+    assert output['tcorr_not_water'] == expected
 
 
 def test_Image_calculate_properties():