From 810466cd94367d2dbf1e8a01bbf4b298543a823d Mon Sep 17 00:00:00 2001 From: Charles Morton Date: Wed, 24 Jul 2024 13:47:33 -0700 Subject: [PATCH 1/4] Testing out putting Tcorr water mask code in a dedicated function Renamed "water_mask" to "not_water_mask" so that pixels with a value of 1 are not-water (or land) and pixels with a value of 0 are water. --- openet/ssebop/image.py | 52 +++++--- openet/ssebop/tests/test_c_image.py | 199 +++++++++++++++++++--------- pyproject.toml | 2 +- 3 files changed, 177 insertions(+), 76 deletions(-) diff --git a/openet/ssebop/image.py b/openet/ssebop/image.py index f80ac43..b1878d5 100644 --- a/openet/ssebop/image.py +++ b/openet/ssebop/image.py @@ -475,6 +475,29 @@ 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 used for FANO Tcorr calculation indicating pixels that are not water + + The purpose for this mask is to get water pixels out of the FANO calculation, + but then put the water pixels back in with the original LST value at the end + + Output image will be 0 for water and 1 for not-water (land?) + + """ + + 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 +840,6 @@ def tcorr_FANO(self): coarse_transform = [5000, 0, 15, 0, -5000, 15] coarse_transform100 = [100000, 0, 15, 0, -100000, 15] dt_coeff = 0.125 - ndwi_threshold = -0.15 high_ndvi_threshold = 0.9 water_pct = 10 # max pixels argument for .reduceResolution() @@ -827,26 +849,26 @@ 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']) ) + total_pixels_count = ( - ndvi + self.qa_water_mask .reduceResolution(ee.Reducer.count(), False, m_pixels) .reproject(self.crs, coarse_transform) .updateMask(1).select([0], ['count']) @@ -865,13 +887,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 +905,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 18ac91a..0509ec6 100644 --- a/openet/ssebop/tests/test_c_image.py +++ b/openet/ssebop/tests/test_c_image.py @@ -1,5 +1,5 @@ import datetime -# import pprint +import pprint import ee import pytest @@ -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,6 +348,7 @@ 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], + # TODO: This image collection is deprecated! ['USGS/NED', [-106.03249, 37.17777], 2364.351], ] ) @@ -254,11 +359,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', [ @@ -506,17 +606,18 @@ def test_Image_from_landsat_c2_sr_lst_source_values(): assert output_img.get('lst_source_id').getInfo().startswith(lst_source) -def test_Image_from_landsat_c2_sr_lst_source_missing(): - """Test that the LST is masked if the scene is not present in lst_source""" - # This image does not currently exist in the source collection, - # but if this test stops working check to see if this image was added - image_id = 'LANDSAT/LC08/C02/T1_L2/LC08_031034_20160702' - xy = (-102.08284, 37.81728) - lst_source = 'projects/openet/assets/lst/landsat/c02' - output_img = ssebop.Image.from_landsat_c2_sr(image_id, lst_source=lst_source).lst - output = utils.point_image_value(output_img, xy) - assert output['lst'] == None - assert output_img.get('lst_source_id').getInfo() == 'None' +# # TODO: Find a new missing LST source image +# def test_Image_from_landsat_c2_sr_lst_source_missing(): +# """Test that the LST is masked if the scene is not present in lst_source""" +# # This image does not currently exist in the source collection, +# # but if this test stops working check to see if this image was added +# image_id = 'LANDSAT/LC08/C02/T1_L2/LC08_031034_20160702' +# xy = (-102.08284, 37.81728) +# lst_source = 'projects/openet/assets/lst/landsat/c02' +# output_img = ssebop.Image.from_landsat_c2_sr(image_id, lst_source=lst_source).lst +# output = utils.point_image_value(output_img, xy) +# assert output['lst'] == None +# assert output_img.get('lst_source_id').getInfo() == 'None' # # DEADBEEF - Keep for now in case approach changes for handling missing scenes in LST source @@ -873,32 +974,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', [ @@ -906,7 +981,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 @@ -914,15 +989,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(): diff --git a/pyproject.toml b/pyproject.toml index 2f9678f..1487df3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "openet-ssebop" -version = "0.4.4" +version = "0.4.5" authors = [ { name = "Gabe Parrish", email = "gparrish@contractor.usgs.gov" }, { name = "Mac Friedrichs", email = "mfriedrichs@contractor.usgs.gov" }, From 9a9c36336dc8a550e81c7105271ca66eec8d43a9 Mon Sep 17 00:00:00 2001 From: Charles Morton Date: Wed, 24 Jul 2024 14:16:29 -0700 Subject: [PATCH 2/4] NED collection in elevation source test is deprecated --- openet/ssebop/tests/test_c_image.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openet/ssebop/tests/test_c_image.py b/openet/ssebop/tests/test_c_image.py index 0509ec6..f2315b8 100644 --- a/openet/ssebop/tests/test_c_image.py +++ b/openet/ssebop/tests/test_c_image.py @@ -348,8 +348,9 @@ def test_Image_elev_band_name(): ['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], - # TODO: This image collection is deprecated! - ['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): From 289b7c18d66bbd627af83cda9ce1a9578b98279c Mon Sep 17 00:00:00 2001 From: Charles Morton Date: Fri, 25 Oct 2024 14:46:52 -0700 Subject: [PATCH 3/4] Updating tcorr_not_water_mask function doc string --- openet/ssebop/image.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/openet/ssebop/image.py b/openet/ssebop/image.py index 9adef93..d284755 100644 --- a/openet/ssebop/image.py +++ b/openet/ssebop/image.py @@ -477,15 +477,18 @@ def quality(self): @lazy_property def tcorr_not_water_mask(self): - """Mask used for FANO Tcorr calculation indicating pixels that are not water + """Mask of pixels that have a high confidence of not being water - The purpose for this mask is to get water pixels out of the FANO calculation, - but then put the water pixels back in with the original LST value at the end + The purpose for this mask is to ensure that water pixels are not used in + the Tcorr FANO calculation. - Output image will be 0 for water and 1 for not-water (land?) + 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 From 0dd183a7c3cd3026196221fb4043d4446929ea34 Mon Sep 17 00:00:00 2001 From: Charles Morton Date: Sat, 26 Oct 2024 10:43:22 -0700 Subject: [PATCH 4/4] Update image.py --- openet/ssebop/image.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openet/ssebop/image.py b/openet/ssebop/image.py index d284755..f068970 100644 --- a/openet/ssebop/image.py +++ b/openet/ssebop/image.py @@ -870,8 +870,9 @@ def tcorr_FANO(self): .updateMask(1).select([0], ['count']) ) + # TODO: Maybe chance ndvi to self.qa_water_mask? total_pixels_count = ( - self.qa_water_mask + ndvi .reduceResolution(ee.Reducer.count(), False, m_pixels) .reproject(self.crs, coarse_transform) .updateMask(1).select([0], ['count'])