Servir-Mekong / hydra-floods

HYDrologic Remote sensing Analysis for Floods Python package
https://servir-mekong.github.io/hydra-floods/
GNU General Public License v3.0
165 stars 49 forks source link

Bug: New version (2023.10.14) does not handle Landsat8 dataset correctly. #48

Open mn5hk opened 10 months ago

mn5hk commented 10 months ago

Hi team, the new version appears to have problems when taking reducers with the Landsat 8 data.

Sharing a code snippet here:

#Acquire optical data using hf.dataset
optical = hf.Landsat8(tonlesap, start_date, end_date)
elv = ee.Image("JAXA/ALOS/AW3D30/V2_2").select("AVE_DSM")

opt_corrected = optical.apply_func(corrections.illumination_correction, elevation = elv)

opt_reduced = opt_corrected.collection.median()

Map = geemap.Map(center = (12.7, 105), zoom = 8)

Map.addLayer(tonlesap, {}, 'Region of Interest')
Map.addLayer(opt_reduced, opt_params, 'Optical Imagery')
Map.addLayerControl()
Map

Also attached are the snippets of the error msg: image image

KMarkert commented 10 months ago

Can you please share a full working example for me to replicate and investigate? What are the variables tonlesap, start_date, end_date?

mn5hk commented 10 months ago

Hi @KMarkert, thanks for a quick response. Please find an extended code snippet with roi and timeslice definitions:

#Define the region of interest
tonlesap = ee.Geometry.Polygon([
        [106.40316, 11.7778],
        [106.40316, 13.6104],
        [103.56869, 13.6104],
        [103.56869, 11.7778]
    ])

#Define time period of interest
start_date = "2021-09-26"
end_date = "2021-10-15"

opt_params = {
    "min": 50,
    "max": 5500,
    "bands": "red,green,blue",
    "gamma":1.5,
    "dimensions": 1024
}

#Acquire optical data using hf.dataset
optical = hf.Landsat8(tonlesap, start_date, end_date)
elv = ee.Image("JAXA/ALOS/AW3D30/V2_2").select("AVE_DSM")

opt_corrected = optical.apply_func(corrections.illumination_correction, elevation = elv)

opt_reduced = opt_corrected.collection.median()

Map = geemap.Map(center = (12.7, 105), zoom = 8)

Map.addLayer(tonlesap, {}, 'Region of Interest')
Map.addLayer(opt_reduced, opt_params, 'Optical Imagery')
Map.addLayerControl()
Map
KMarkert commented 10 months ago

Found the source of the error. There is a lookup for sun angle properties in the illumination_correction function (see here).

The latest update upgraded the Landsat collections to use Collection 2 (C1 is deprecated) and there were changes to property names including the ones that are looked up by the function which causes it to fail.

here is a quick fix that worked for me:

#Acquire optical data using hf.dataset
optical = hf.Landsat8(tonlesap, start_date, end_date)

# function to rename the sun angle properties to what is expected by
# the illumination_correction function
def rename_sunangle_properties(image):
    sz = ee.Number(90).subtract(image.get('SUN_ELEVATION'))
    sa = image.get('SUN_AZIMUTH')
    return image.set({"SOLAR_ZENITH_ANGLE": sz, "SOLAR_AZIMUTH_ANGLE": sa})

# apply the function to rename properties
optical = optical.apply_func(rename_sunangle_properties)

elv = ee.Image("JAXA/ALOS/AW3D30/V2_2").select("AVE_DSM").unmask(0)

opt_corrected = optical.apply_func(hf.corrections.illumination_correction, elevation = elv, scale=180)

Longer term fix will be to update the illumination_correction function lookup section to look for the correct properties. That should be a relatively straightforward fix.

mn5hk commented 10 months ago

Hi Kel, the quick fix worked for me. Thanks a lot for your prompt support and a detailed explanation to the issue.

KMarkert commented 10 months ago

No problem!

I am reopening this issue so we have in in the backlog for things to fix.