diff --git a/openet/ssebop/image.py b/openet/ssebop/image.py index b18e67b..f068970 100644 --- a/openet/ssebop/image.py +++ b/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) ) diff --git a/openet/ssebop/tests/test_c_image.py b/openet/ssebop/tests/test_c_image.py index 1dc4faf..1a32513 100644 --- a/openet/ssebop/tests/test_c_image.py +++ b/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,32 +979,6 @@ 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', [ @@ -911,7 +986,7 @@ def test_Image_time_properties(): ['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 @@ -919,15 +994,19 @@ def test_Image_qa_water_mask(image_id, xy, expected): 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():