Normalized difference calculation fixes

openet/ssebop/image.py

@@ -602,7 +602,11 @@ def tcorr(self):
                 .set({'tcorr_source': f'custom_{self._tcorr_source}'})
             )
         elif 'FANO' == self._tcorr_source.upper():
-            return ee.Image(self.tcorr_FANO).select(['tcorr']).set({'tcorr_source': 'FANO'})
+            return (
+                ee.Image(self.tcorr_FANO).select(['tcorr'])
+                .updateMask(1)
+                .set({'tcorr_source': 'FANO'})
+            )
         else:
             raise ValueError(f'Unsupported tcorr_source: {self._tcorr_source}\n')
 
@@ -849,7 +853,7 @@ def tcorr_FANO(self):
         m_pixels = 65535
 
         lst = ee.Image(self.lst)
-        ndvi = ee.Image(self.ndvi).clamp(-1.0, 1.0)
+        ndvi = ee.Image(self.ndvi)
         tmax = ee.Image(self.tmax)
         dt = ee.Image(self.dt)
 
@@ -858,7 +862,7 @@ def tcorr_FANO(self):
         qa_watermask = ee.Image(self.qa_water_mask)
         ndvi = ndvi.where(qa_watermask.eq(1).And(ndvi.gt(0)), ndvi.multiply(-1))
 
-        # Mask with not_water pixels set to 1 and water pixels set to 0
+        # Mask with not_water pixels set to 1 and other (likely water) pixels set to 0
         not_water_mask = self.tcorr_not_water_mask
 
         # Count not-water pixels and the total number of pixels

openet/ssebop/landsat.py

@@ -111,20 +111,53 @@ def lst(landsat_image):
     return lst.rename(['lst'])
 
 
-def ndvi(landsat_image):
+def ndvi(landsat_image, gsw_extent_flag=True):
     """Normalized difference vegetation index
 
     Parameters
     ----------
     landsat_image : ee.Image
         "Prepped" Landsat image with standardized band names.
+    gsw_extent_flag : boolean
+        If True, apply the global surface water extent mask to the QA_PIXEL water mask
+        The default is True.
 
     Returns
     -------
     ee.Image
 
     """
-    return ee.Image(landsat_image).normalizedDifference(['nir', 'red']).rename(['ndvi'])
+    # Force the input values to be at greater than or equal to zero
+    #   since C02 surface reflectance values can be negative
+    #   but the normalizedDifference function will return nodata
+    ndvi_img = landsat_image.max(0).normalizedDifference(['nir', 'red'])
+
+    b1 = landsat_image.select(['nir'])
+    b2 = landsat_image.select(['red'])
+
+    # Assume that very high reflectance values are unreliable for computing the index
+    #   and set the output value to 0
+    # Threshold value could be set lower, but for now only trying to catch saturated pixels
+    ndvi_img = ndvi_img.where(b1.gte(1).Or(b2.gte(1)), 0)
+
+    # Including the global surface water maximum extent to help remove shadows that
+    #   are misclassified as water
+    # The flag is needed so that the image can be bypassed during testing with constant images
+    qa_water_mask = landsat_c2_qa_water_mask(landsat_image)
+    if gsw_extent_flag:
+        gsw_mask = ee.Image('JRC/GSW1_4/GlobalSurfaceWater').select(['max_extent']).gte(1)
+        qa_water_mask = qa_water_mask.And(gsw_mask)
+
+    # Assume that low reflectance values are unreliable for computing the index
+    # If both reflectance values are below the threshold,
+    #   and if the pixel is flagged as water, set the output to -0.1 (should this be -1?)
+    #   otherwise set the output to 0
+    ndvi_img = ndvi_img.where(b1.lt(0.01).And(b2.lt(0.01)), 0)
+    ndvi_img = ndvi_img.where(b1.lt(0.01).And(b2.lt(0.01)).And(qa_water_mask), -0.1)
+    # Should there be an additional check for if either value was negative?
+    # ndvi_img = ndvi_img.where(b1.lt(0).Or(b2.lt(0)), 0)
+
+    return ndvi_img.clamp(-1.0, 1.0).rename(['ndvi'])
 
 
 def ndwi(landsat_image):
@@ -140,7 +173,25 @@ def ndwi(landsat_image):
     ee.Image
 
     """
-    return ee.Image(landsat_image).normalizedDifference(['green', 'swir1']).rename(['ndwi'])
+    # Force the input values to be at greater than or equal to zero
+    #   since C02 surface reflectance values can be negative
+    #   but the normalizedDifference function will return nodata
+    ndwi_img = landsat_image.max(0).normalizedDifference(['green', 'swir1'])
+
+    b1 = landsat_image.select(['green'])
+    b2 = landsat_image.select(['swir1'])
+
+    # Assume that very high reflectance values are unreliable for computing the index
+    #   and set the output value to 0
+    # Threshold value could be set lower, but for now only trying to catch saturated pixels
+    ndwi_img = ndwi_img.where(b1.gte(1).Or(b2.gte(1)), 0)
+
+    # Assume that low reflectance values are unreliable for computing the index
+    # If both reflectance values are below the threshold set the output to 0
+    # May want to check the QA water mask here also, similar to the NDVI calculation
+    ndwi_img = ndwi_img.where(b1.lt(0.01).And(b2.lt(0.01)), 0)
+
+    return ndwi_img.clamp(-1.0, 1.0).rename(['ndwi'])
 
 
 def landsat_c2_qa_water_mask(landsat_image):
@@ -149,7 +200,7 @@ def landsat_c2_qa_water_mask(landsat_image):
     Parameters
     ----------
     landsat_image : ee.Image
-        Image from a Landsat C02 image collection with a QA_PIXEL band.
+        Landsat C02 image with a QA_PIXEL band.
 
     Returns
     -------

openet/ssebop/tests/test_b_landsat.py

@@ -10,13 +10,11 @@
 # import openet.core.utils as utils
 
 
-# Should these be test fixtures instead?
-# I'm not sure how to make them fixtures and allow input parameters
-def toa_image(red=0.1, nir=0.9, bt=305):
-    """Construct a fake Landsat 8 TOA image with renamed bands"""
+def sr_image(blue=0.1, green=0.1, red=0.1, nir=0.9, swir1=0.1, swir2=0.1, bt=305, qa=1):
+    """Construct a fake Landsat 8 SR image with renamed bands"""
     return (
-        ee.Image.constant([red, nir, bt])
-        .rename(['red', 'nir', 'tir'])
+        ee.Image.constant([blue, green, red, nir, swir1, swir2, bt, qa])
+        .rename(['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'tir', 'QA_PIXEL'])
         .set({
             # 'system:time_start': ee.Date(SCENE_DATE).millis(),
             'k1_constant': ee.Number(607.76),
@@ -25,81 +23,173 @@ def toa_image(red=0.1, nir=0.9, bt=305):
     )
 
 
-def sr_image(red=1000, nir=9000, bt=305):
-    """Construct a fake Landsat 8 TOA image with renamed bands"""
-    return (
-        ee.Image.constant([red, nir, bt])
-        .rename(['red', 'nir', 'tir'])
-        .set({
-            # 'system:time_start': ee.Date(SCENE_DATE).millis(),
-            'k1_constant': ee.Number(607.76),
-            'k2_constant': ee.Number(1260.56),
-        })
-    )
+def test_ndvi_band_name():
+    output = utils.getinfo(landsat.ndvi(sr_image()))
+    assert output['bands'][0]['id'] == 'ndvi'
 
 
 @pytest.mark.parametrize(
     'red, nir, expected',
     [
-        [0.2, 9.0 / 55, -0.1],
-        [0.2, 0.2,  0.0],
-        [0.1, 11.0 / 90,  0.1],
-        [0.2, 0.3, 0.2],
-        [0.1, 13.0 / 70, 0.3],
-        [0.3, 0.7, 0.4],
-        [0.2, 0.6, 0.5],
-        [0.2, 0.8, 0.6],
-        [0.1, 17.0 / 30, 0.7],
-        [0.1, 0.9, 0.8],
+        [0.02, 0.9 / 55, -0.1],
+        [0.02, 0.02,  0.0],
+        [0.01, 0.11 / 9, 0.1],
+        [0.02, 0.03, 0.2],
+        [0.01, 0.13 / 7, 0.3],
+        [0.03, 0.07, 0.4],
+        [0.02, 0.06, 0.5],
+        [0.02, 0.08, 0.6],
+        [0.01, 0.17 / 3, 0.7],
+        [0.01, 0.09, 0.8],
     ]
 )
 def test_ndvi_calculation(red, nir, expected, tol=0.000001):
-    toa = toa_image(red=red, nir=nir)
-    output = utils.constant_image_value(landsat.ndvi(toa))
+    output = utils.constant_image_value(landsat.ndvi(sr_image(red=red, nir=nir)))
     assert abs(output['ndvi'] - expected) <= tol
 
 
-def test_ndvi_band_name():
-    output = utils.getinfo(landsat.ndvi(toa_image()))
-    assert output['bands'][0]['id'] == 'ndvi'
+@pytest.mark.parametrize(
+    'red, nir, expected',
+    [
+        [1.0, 0.4, 0.0],
+        [0.4, 1.0, 0.0],
+        [1.0, 1.0, 0.0],
+    ]
+)
+def test_ndvi_saturated_reflectance(red, nir, expected, tol=0.000001):
+    # Check that saturated reflectance values return 0
+    output = utils.constant_image_value(landsat.ndvi(sr_image(red=red, nir=nir)))
+    assert abs(output['ndvi'] - expected) <= tol
 
 
 @pytest.mark.parametrize(
     'red, nir, expected',
     [
-        [0.2, 9.0 / 55, 0.985],      # -0.1
-        [0.2, 0.2,  0.977],          # 0.0
-        [0.1, 11.0 / 90,  0.977],    # 0.1
-        [0.2, 0.2999, 0.977],        # 0.3- (0.3 NIR isn't exactly an NDVI of 0.2)
-        [0.2, 0.3001, 0.986335],     # 0.3+
-        [0.1, 13.0 / 70, 0.986742],  # 0.3
-        [0.3, 0.7, 0.987964],        # 0.4
-        [0.2, 0.6, 0.99],            # 0.5
-        [0.2, 0.8, 0.99],            # 0.6
-        [0.1, 17.0 / 30, 0.99],      # 0.7
+        [-0.1, -0.1, 0.0],
+        [0.0, 0.0, 0.0],
+        [0.009, 0.009, 0.0],
+        [0.009, -0.01, 0.0],
+        [-0.01, 0.009, 0.0],
+        # Check that calculation works correctly if one value is above threshold
+        [-0.01, 0.1, 1.0],
+        [0.1, -0.01, -1.0],
     ]
 )
-def test_emissivity_calculation(red, nir, expected, tol=0.000001):
-    output = utils.constant_image_value(landsat.emissivity(toa_image(red=red, nir=nir)))
-    assert abs(output['emissivity'] - expected) <= tol
+def test_ndvi_negative_non_water(red, nir, expected, tol=0.000001):
+    # Check that non-water pixels with very low or negative reflectance values are set to 0.0
+    output = utils.constant_image_value(landsat.ndvi(sr_image(red=red, nir=nir, qa=1)))
+    assert abs(output['ndvi'] - expected) <= tol
+
+
+@pytest.mark.parametrize(
+    'red, nir, expected',
+    [
+        [-0.1, -0.1, -0.1],
+        [0.0, 0.0, -0.1],
+        [0.009, 0.009, -0.1],
+        [0.009, -0.01, -0.1],
+        [-0.01, 0.009, -0.1],
+    ]
+)
+def test_ndvi_negative_water(red, nir, expected, tol=0.000001):
+    # Check that water pixels with very low or negative reflectance values are set to -0.1
+    output = utils.constant_image_value(landsat.ndvi(
+        sr_image(red=red, nir=nir, qa=128), gsw_extent_flag=False
+    ))
+    assert abs(output['ndvi'] - expected) <= tol
+
+
+def test_ndwi_band_name():
+    output = utils.getinfo(landsat.ndwi(sr_image()))
+    assert output['bands'][0]['id'] == 'ndwi'
+
+
+@pytest.mark.parametrize(
+    'green, swir1, expected',
+    [
+        [0.01, 0.07, -0.75],
+        [0.01, 0.03, -0.5],
+        [0.9 / 55, 0.02, -0.1],
+        [0.2, 0.2, 0.0],
+        [0.11 / 9, 0.01, 0.1],
+        [0.07, 0.03, 0.4],
+        [0.09, 0.01, 0.8],
+    ]
+)
+def test_ndwi_calculation(green, swir1, expected, tol=0.000001):
+    output = utils.constant_image_value(landsat.ndwi(sr_image(green=green, swir1=swir1)))
+    assert abs(output['ndwi'] - expected) <= tol
+
+
+@pytest.mark.parametrize(
+    'green, swir1, expected',
+    [
+        [-0.1, -0.1, 0.0],
+        [0.0, 0.0, 0.0],
+        [0.009, 0.009, 0.0],
+        [0.009, -0.01, 0.0],
+        [-0.01, 0.009, 0.0],
+        # Check that calculation works correctly if one value is above threshold
+        [-0.01, 0.1, -1.0],
+        [0.1, -0.01, 1.0],
+    ]
+)
+def test_ndwi_negative_reflectance(green, swir1, expected, tol=0.000001):
+    # Check that very low or negative reflectance values are set to 0
+    output = utils.constant_image_value(landsat.ndwi(sr_image(green=green, swir1=swir1)))
+    assert abs(output['ndwi'] - expected) <= tol
+
+
+@pytest.mark.parametrize(
+    'green, swir1, expected',
+    [
+        [1.0, 0.4, 0.0],
+        [0.4, 1.0, 0.0],
+        [1.0, 1.0, 0.0],
+    ]
+)
+def test_ndwi_saturated_reflectance(green, swir1, expected, tol=0.000001):
+    # Check that saturated reflectance values return 0
+    output = utils.constant_image_value(landsat.ndwi(sr_image(green=green, swir1=swir1)))
+    assert abs(output['ndwi'] - expected) <= tol
 
 
 def test_emissivity_band_name():
-    output = utils.getinfo(landsat.emissivity(toa_image()))
+    output = utils.getinfo(landsat.emissivity(sr_image()))
     assert output['bands'][0]['id'] == 'emissivity'
 
 
 @pytest.mark.parametrize(
-    'red, nir, bt, expected',
+    'red, nir, expected',
     [
-        [0.2, 0.7, 300, 303.471031],
+        [0.02, 0.9 / 55, 0.985],     # -0.1
+        [0.02, 0.02,  0.977],        # 0.0
+        [0.01, 0.11 / 9,  0.977],    # 0.1
+        [0.02, 0.02999, 0.977],      # 0.3- (0.3 NIR isn't exactly an NDVI of 0.2)
+        [0.02, 0.03001, 0.986335],   # 0.3+
+        [0.01, 0.13 / 7, 0.986742],  # 0.3
+        [0.03, 0.07, 0.987964],      # 0.4
+        [0.02, 0.06, 0.99],          # 0.5
+        [0.02, 0.08, 0.99],          # 0.6
+        [0.01, 0.17 / 3, 0.99],      # 0.7
     ]
 )
-def test_lst_calculation(red, nir, bt, expected, tol=0.000001):
-    output = utils.constant_image_value(landsat.lst(toa_image(red=red, nir=nir, bt=bt)))
-    assert abs(output['lst'] - expected) <= tol
+def test_emissivity_calculation(red, nir, expected, tol=0.000001):
+    output = utils.constant_image_value(landsat.emissivity(sr_image(red=red, nir=nir)))
+    assert abs(output['emissivity'] - expected) <= tol
 
 
 def test_lst_band_name():
-    output = utils.getinfo(landsat.lst(toa_image()))
+    output = utils.getinfo(landsat.lst(sr_image()))
     assert output['bands'][0]['id'] == 'lst'
+
+
+@pytest.mark.parametrize(
+    'red, nir, bt, expected',
+    [
+        [0.02, 0.07, 300, 303.471031],
+    ]
+)
+def test_lst_calculation(red, nir, bt, expected, tol=0.000001):
+    output = utils.constant_image_value(landsat.lst(sr_image(red=red, nir=nir, bt=bt)))
+    assert abs(output['lst'] - expected) <= tol