From a185093e28cc8c8871178e0157f5fdd78c85079e Mon Sep 17 00:00:00 2001 From: kmarkert Date: Thu, 14 Apr 2022 11:26:03 -0600 Subject: [PATCH 1/7] adding wip mods to package --- hydrafloods/MODIS_DNNS.py | 149 ++++++++++--------- hydrafloods/VIIRS_DNNS.py | 151 ++++++++++--------- hydrafloods/__init__.py | 2 + hydrafloods/defaults.py | 35 +++++ hydrafloods/defaults.yml | 21 +++ hydrafloods/fusion.py | 117 ++++++++++++++- hydrafloods/tests/test_datasets.py | 65 ++++++--- hydrafloods/workflows/dswfp.py | 22 +-- hydrafloods/workflows/hydrosar.py | 227 +++++++++++++++++++++++++++++ 9 files changed, 614 insertions(+), 175 deletions(-) create mode 100644 hydrafloods/defaults.py create mode 100644 hydrafloods/defaults.yml create mode 100644 hydrafloods/workflows/hydrosar.py diff --git a/hydrafloods/MODIS_DNNS.py b/hydrafloods/MODIS_DNNS.py index b304d34..abb3df8 100644 --- a/hydrafloods/MODIS_DNNS.py +++ b/hydrafloods/MODIS_DNNS.py @@ -5,106 +5,119 @@ import math import matplotlib.pyplot as pl + def perm_water_mask(): - return ee.Image("MODIS/MOD44W/MOD44W_005_2000_02_24").select(['water_mask'], ['b1']) + return ee.Image("MODIS/MOD44W/MOD44W_005_2000_02_24").select(["water_mask"], ["b1"]) + def DEM(): - return ee.Image('CGIAR/SRTM90_V4') -#ALOS DEM -#ee.Image('JAXA/ALOS/AW3D30_V1_1') -#SRTM DEM -#ee.Image('USGS/SRTMGL1_003') - - -def GEE_classifier(img, name_classifier, arguments={}): - arguments = { - 'training_image' : perm_water_mask(), - 'training_band' : "b1", - 'training_region' : img - } - c_arguments = { - 'subsampling' : 0.1, - 'max_classification': 2, - 'classifier_mode' : 'classification', - 'name_classifier' : name_classifier - } - arguments.update(c_arguments) - classes1 = img.trainClassifier(**arguments) - classes2 = classify(classes1).select(['classification'], ['b1']) - return classes2; + return ee.Image("CGIAR/SRTM90_V4") + + +# ALOS DEM +# ee.Image('JAXA/ALOS/AW3D30_V1_1') +# SRTM DEM +# ee.Image('USGS/SRTMGL1_003') + + +# def GEE_classifier(img, name_classifier, arguments={}): +# arguments = { +# "training_image": perm_water_mask(), +# "training_band": "b1", +# "training_region": img, +# } +# c_arguments = { +# "subsampling": 0.1, +# "max_classification": 2, +# "classifier_mode": "classification", +# "name_classifier": name_classifier, +# } +# arguments.update(c_arguments) +# classes1 = img.trainClassifier(**arguments) +# classes2 = classify(classes1).select(["classification"], ["b1"]) +# return classes2 def dnns(img): - + kernel_scale = 40 - VIIRS_ave_scale = 500 - - k1 = ee.Kernel.square(kernel_scale, 'pixels', False) - - composite = img['b1'].addBands(img['b2']).addBands(img['b6']) - - classes = GEE_classifier(img, 'Pegasos', {'classifier_mode' : 'probability'}) + VIIRS_ave_scale = 500 + + k1 = ee.Kernel.square(kernel_scale, "pixels", False) + + composite = img["b1"].addBands(img["b2"]).addBands(img["b6"]) + + # classes = GEE_classifier(img, "Pegasos", {"classifier_mode": "probability"}) + classes = ee.Image([0.5, 0.5]) water = classes.gte(0.9) - land = classes.lte(0.1) + land = classes.lte(0.1) mix = water.Not().And(land.Not()) ave_water = water.mask(water).multiply(composite) - ave_water_img = ee.Image([ave_water['b1'], ave_water['b2'], ave_water['b6']]) - + ave_water_img = ee.Image([ave_water["b1"], ave_water["b2"], ave_water["b6"]]) + N_nmin_water = 1 N_water = water.convolve(k1) - - water_rf = water.multiply(composite).convolve(k1).multiply(N_water.gte(N_nmin_water)).divide(N_water) + + water_rf = ( + water.multiply(composite) + .convolve(k1) + .multiply(N_water.gte(N_nmin_water)) + .divide(N_water) + ) water_rf = water_rf.add(ave_water_img.multiply(water_rf.Not())) ave_land = composite.mask(land) - - ave_land_img = ee.Image([ave_land['b1'], ave_land['b2'], ave_land['b6']]) - - # Constraints - R1 = img['b1'].divide(img['b6']) - R2 = img['b2'].divide(img['b6']) - R3 = R1.subtract(water_rf.select('b1').divide(img['b6']) ) - R4 = R2.subtract(water_rf.select('b2').divide(img['b6']) ) - - NR1 = R1.neighborhoodToBands(k1) + + ave_land_img = ee.Image([ave_land["b1"], ave_land["b2"], ave_land["b6"]]) + + # Constraints + R1 = img["b1"].divide(img["b6"]) + R2 = img["b2"].divide(img["b6"]) + R3 = R1.subtract(water_rf.select("b1").divide(img["b6"])) + R4 = R2.subtract(water_rf.select("b2").divide(img["b6"])) + + NR1 = R1.neighborhoodToBands(k1) NR2 = R2.neighborhoodToBands(k1) - NI1 = img['b1'].neighborhoodToBands(k1) - NI2 = img['b2'].neighborhoodToBands(k1) - NI3 = img['b6'].neighborhoodToBands(k1) - - + NI1 = img["b1"].neighborhoodToBands(k1) + NI2 = img["b2"].neighborhoodToBands(k1) + NI3 = img["b6"].neighborhoodToBands(k1) + M1 = (NR1.gt(R3)).And(NR1.lt(R1)) M2 = (NR2.gt(R4)).And(NR2.lt(R2)) nLP = M1.And(M2) - + NnLP = nLP.reduce(ee.Reducer.sum()) - ave_nI1 = NI1.multiply(nLP).reduce(ee.Reducer.sum()).divide(numnLP) - ave_nI2 = NI2.multiply(nLP).reduce(ee.Reducer.sum()).divide(numnLP) + ave_nI1 = NI1.multiply(nLP).reduce(ee.Reducer.sum()).divide(NnLP) + ave_nI2 = NI2.multiply(nLP).reduce(ee.Reducer.sum()).divide(NnLP) ave_nI3 = NI3.multiply(nLP).reduce(ee.Reducer.sum()).divide(NnLP) N_nmin_land = 1 ave_land = ave_nI1.addBands(ave_nI2).addBands(ave_nI3) - ave_land = ave_land.multiply(NnLP.gte(N_nmin_land)).add(ave_land_img.multiply(NnLP.lt(N_nmin_land)) ) + ave_land = ave_land.multiply(NnLP.gte(N_nmin_land)).add( + ave_land_img.multiply(NnLP.lt(N_nmin_land)) + ) - ave_landI3 = ave_land.select('b6') - f_water = (ave_landI3.subtract(img['b6'])).divide(ave_landI3.subtract(water_rf.select('b6'))).clamp(0, 1) + ave_landI3 = ave_land.select("b6") + f_water = ( + (ave_landI3.subtract(img["b6"])) + .divide(ave_landI3.subtract(water_rf.select("b6"))) + .clamp(0, 1) + ) f_water = f_water.add(water).subtract(land).clamp(0, 1) - + return f_water def DEM_downscale(img, f_water): - + MODIS_Pixel_meters = 500 DEM_d = img.DEM() - + water_present = f_water.gt(0.0) - h_min = DEM_d.mask(f_water).focal_min(MODIS_Pixel_meters, 'square', 'meters') - h_max = DEM_d.mask(f_water).focal_max(MODIS_Pixel_meters, 'square', 'meters') - - - water_high = h_min.add(h_max.subtract(h_min).multiply(f_water)) - water_high = water_high.multiply(f_water.lt(1.0)).multiply(f_water.gt(0.0)) - return f_water.eq(1.0).select(['elevation'], ['b1']) + h_min = DEM_d.mask(f_water).focal_min(MODIS_Pixel_meters, "square", "meters") + h_max = DEM_d.mask(f_water).focal_max(MODIS_Pixel_meters, "square", "meters") + water_high = h_min.add(h_max.subtract(h_min).multiply(f_water)) + water_high = water_high.multiply(f_water.lt(1.0)).multiply(f_water.gt(0.0)) + return f_water.eq(1.0).select(["elevation"], ["b1"]) diff --git a/hydrafloods/VIIRS_DNNS.py b/hydrafloods/VIIRS_DNNS.py index 165dffe..cfefe7c 100644 --- a/hydrafloods/VIIRS_DNNS.py +++ b/hydrafloods/VIIRS_DNNS.py @@ -1,5 +1,4 @@ - -# coding: utf-8 VIIRS DNNS +# coding: utf-8 VIIRS DNNS import ee import numpy as np @@ -8,107 +7,117 @@ def perm_water_mask(): - return ee.Image("MODIS/MOD44W/MOD44W_005_2000_02_24").select(['water_mask'], ['b1']) + return ee.Image("MODIS/MOD44W/MOD44W_005_2000_02_24").select(["water_mask"], ["b1"]) def DEM(): - return ee.Image('CGIAR/SRTM90_V4') -#ALOS DEM -#ee.Image('JAXA/ALOS/AW3D30_V1_1') -#SRTM DEM -#ee.Image('USGS/SRTMGL1_003') - - -def GEE_classifier(img, name_classifier, arguments={}): - arguments = { - 'training_image' : perm_water_mask(), - 'training_band' : "b1", - 'training_region' : img - } - c_arguments = { - 'subsampling' : 0.1, - 'max_classification': 2, - 'classifier_mode' : 'classification', - 'name_classifier' : name_classifier - } - arguments.update(c_arguments) - classes1 = img.trainClassifier(**arguments) - classes2 = classify(classes1).select(['classification'], ['b1']) - return classes2; + return ee.Image("CGIAR/SRTM90_V4") + + +# ALOS DEM +# ee.Image('JAXA/ALOS/AW3D30_V1_1') +# SRTM DEM +# ee.Image('USGS/SRTMGL1_003') +# def GEE_classifier(img, name_classifier, arguments={}): +# arguments = { +# "training_image": perm_water_mask(), +# "training_band": "b1", +# "training_region": img, +# } +# c_arguments = { +# "subsampling": 0.1, +# "max_classification": 2, +# "classifier_mode": "classification", +# "name_classifier": name_classifier, +# } +# arguments.update(c_arguments) +# classes1 = img.trainClassifier(**arguments) +# classes2 = classify(classes1).select(["classification"], ["b1"]) +# return classes2 + def dnns(img): - + kernel_scale = 40 - VIIRS_ave_scale = 500 - - k1 = ee.Kernel.square(kernel_scale, 'pixels', False) - - composite = img['I1'].addBands(img['I2']).addBands(img['I3']) - - classes = GEE_classifier(img, 'Pegasos', {'classifier_mode' : 'probability'}) + VIIRS_ave_scale = 500 + + k1 = ee.Kernel.square(kernel_scale, "pixels", False) + + composite = img["I1"].addBands(img["I2"]).addBands(img["I3"]) + + # classes = GEE_classifier(img, "Pegasos", {"classifier_mode": "probability"}) + classes = ee.Image([0.5, 0.5]) water = classes.gte(0.9) - land = classes.lte(0.1) + land = classes.lte(0.1) mix = water.Not().And(land.Not()) ave_water = water.mask(water).multiply(composite) - ave_water_img = ee.Image([ave_water['I1'], ave_water['I2'], ave_water['I3']]) - + ave_water_img = ee.Image([ave_water["I1"], ave_water["I2"], ave_water["I3"]]) + N_nmin_water = 1 N_water = water.convolve(k1) - - water_rf = water.multiply(composite).convolve(k1).multiply(N_water.gte(N_nmin_water)).divide(N_water) + + water_rf = ( + water.multiply(composite) + .convolve(k1) + .multiply(N_water.gte(N_nmin_water)) + .divide(N_water) + ) water_rf = water_rf.add(ave_water_img.multiply(water_rf.Not())) ave_land = composite.mask(land) - - ave_land_img = ee.Image([ave_land['I1'], ave_land['I2'], ave_land['I3']]) - - # Constraints - R1 = img['I1'].divide(img['I3']) - R2 = img['I2'].divide(img['I3']) - R3 = R1.subtract(water_rf.select('I1').divide(img['I3']) ) - R4 = R2.subtract(water_rf.select('I2').divide(img['I3']) ) - - NR1 = R1.neighborhoodToBands(k1) + + ave_land_img = ee.Image([ave_land["I1"], ave_land["I2"], ave_land["I3"]]) + + # Constraints + R1 = img["I1"].divide(img["I3"]) + R2 = img["I2"].divide(img["I3"]) + R3 = R1.subtract(water_rf.select("I1").divide(img["I3"])) + R4 = R2.subtract(water_rf.select("I2").divide(img["I3"])) + + NR1 = R1.neighborhoodToBands(k1) NR2 = R2.neighborhoodToBands(k1) - NI1 = img['I1'].neighborhoodToBands(k1) - NI2 = img['I2'].neighborhoodToBands(k1) - NI3 = img['I3'].neighborhoodToBands(k1) - - + NI1 = img["I1"].neighborhoodToBands(k1) + NI2 = img["I2"].neighborhoodToBands(k1) + NI3 = img["I3"].neighborhoodToBands(k1) + M1 = (NR1.gt(R3)).And(NR1.lt(R1)) M2 = (NR2.gt(R4)).And(NR2.lt(R2)) nLP = M1.And(M2) - + NnLP = nLP.reduce(ee.Reducer.sum()) - ave_nI1 = NI1.multiply(nLP).reduce(ee.Reducer.sum()).divide(numnLP) - ave_nI2 = NI2.multiply(nLP).reduce(ee.Reducer.sum()).divide(numnLP) + ave_nI1 = NI1.multiply(nLP).reduce(ee.Reducer.sum()).divide(NnLP) + ave_nI2 = NI2.multiply(nLP).reduce(ee.Reducer.sum()).divide(NnLP) ave_nI3 = NI3.multiply(nLP).reduce(ee.Reducer.sum()).divide(NnLP) N_nmin_land = 1 ave_land = ave_nI1.addBands(ave_nI2).addBands(ave_nI3) - ave_land = ave_land.multiply(NnLP.gte(N_nmin_land)).add(ave_land_img.multiply(NnLP.lt(N_nmin_land)) ) - - ave_landI3 = ave_land.select('I3') - f_water = (ave_landI3.subtract(img['I3'])).divide(ave_landI3.subtract(water_rf.select('I3'))).clamp(0, 1) + ave_land = ave_land.multiply(NnLP.gte(N_nmin_land)).add( + ave_land_img.multiply(NnLP.lt(N_nmin_land)) + ) + + ave_landI3 = ave_land.select("I3") + f_water = ( + (ave_landI3.subtract(img["I3"])) + .divide(ave_landI3.subtract(water_rf.select("I3"))) + .clamp(0, 1) + ) f_water = f_water.add(water).subtract(land).clamp(0, 1) - + return f_water def DEM_downscale(img, f_water): - + VIIRS_Pixel_meters = 500 DEM_d = img.DEM() - + water_present = f_water.gt(0.0) - h_min = DEM_d.mask(f_water).focal_min(VIIRS_Pixel_meters, 'square', 'meters') - h_max = DEM_d.mask(f_water).focal_max(VIIRS_Pixel_meters, 'square', 'meters') - - - water_high = h_min.add(h_max.subtract(h_min).multiply(f_water)) - water_high = water_high.multiply(f_water.lt(1.0)).multiply(f_water.gt(0.0)) - return f_water.eq(1.0).select(['elevation'], ['I1']) + h_min = DEM_d.mask(f_water).focal_min(VIIRS_Pixel_meters, "square", "meters") + h_max = DEM_d.mask(f_water).focal_max(VIIRS_Pixel_meters, "square", "meters") + water_high = h_min.add(h_max.subtract(h_min).multiply(f_water)) + water_high = water_high.multiply(f_water.lt(1.0)).multiply(f_water.gt(0.0)) + return f_water.eq(1.0).select(["elevation"], ["I1"]) diff --git a/hydrafloods/__init__.py b/hydrafloods/__init__.py index a6ede88..947f83e 100644 --- a/hydrafloods/__init__.py +++ b/hydrafloods/__init__.py @@ -9,6 +9,8 @@ from hydrafloods.floods import * from hydrafloods import fetch, utils +from hydrafloods import fusion + # from hydrafloods import * __version__ = "2021.11.10" diff --git a/hydrafloods/defaults.py b/hydrafloods/defaults.py new file mode 100644 index 0000000..ed49914 --- /dev/null +++ b/hydrafloods/defaults.py @@ -0,0 +1,35 @@ +s1 = { + "edge_otsu": { + "initial_threshold": -16, + "edge_buffer": 300, + "scale": 150, + "band": "VV", + }, + "bmax_otsu": { + "initial_threshold": -16, + "grid_size": 0.25, + "scale": 150, + "band": "VV", + }, +} + +lc8 = { + "edge_otsu": { + "initial_threshold": 0.0, + "edge_buffer": 300, + "scale": 150, + "invert": True, + "band": "mndwi", + }, + "bmax_otsu": {"initial_threshold": 0.0, "scale": 150, "invert": True}, +} + +fusion = { + "edge_otsu": { + "initial_threshold": 0.0, + "edge_buffer": 300, + "scale": 150, + "invert": True, + }, + "bmax_otsu": {"initial_threshold": 0.0, "scale": 150, "invert": True}, +} diff --git a/hydrafloods/defaults.yml b/hydrafloods/defaults.yml new file mode 100644 index 0000000..7324869 --- /dev/null +++ b/hydrafloods/defaults.yml @@ -0,0 +1,21 @@ +sentinel1: + edge_otsu: + initial_threshold: -16 + edge_buffer: 300 + scale: 150 + + bmax_otsu: + initial_threshold: -16 + scale: 150 + +landsat: + edge_otsu: + initial_threshold: 0.0 + edge_buffer: 300 + scale: 150 + invert: True + + bmax_otsu: + initial_threshold: -16 + scale: 150, + invert: True \ No newline at end of file diff --git a/hydrafloods/fusion.py b/hydrafloods/fusion.py index 2303068..bd7c43c 100644 --- a/hydrafloods/fusion.py +++ b/hydrafloods/fusion.py @@ -1,6 +1,6 @@ from __future__ import print_function, division import ee -from hydrafloods import decorators +from hydrafloods import decorators, ml, thresholding def starfm( @@ -143,3 +143,118 @@ def minimizeDepth(d): ) return final + + +@decorators.keep_attrs +def pca_fusion( + img, + img_stats=None, + match_stats=None, + eigen_vecs=None, + scaling_dict=None, + img_band=None, + match_band=None, + inverse_relationship=True, + reverse_pca_bands=False, +): + pca_bands = [img_band, match_band] + + if reverse_pca_bands: + pca_bands = pca_bands[::-1] + + img_z = ( + img.select(img_band) + .subtract(img_stats.select(img_band + "_mean")) + .divide(img_stats.select(img_band + "_stdDev")) + ) + + if inverse_relationship: + img_z = img_z.multiply(-1) + + match_synth = ( + match_stats.select(match_band + "_mean") + .add(img_z.multiply((match_stats.select(match_band + "_stdDev")))) + .rename(match_band) + ) + + if reverse_pca_bands: + input_img = match_synth.addBands(img) + else: + input_img = img.addBands(match_synth) + + if scaling_dict is not None: + input_img = ml.standard_image_scaling(input_img, scaling_dict, pca_bands) + + pcs = ml.apply_image_pca(input_img, eigen_vecs, pca_bands) + + df_index = pcs.select([0]).clamp(-3.75, 3.75).rename("df_index") + + return df_index + +@decorators.keep_attrs +def dnns(img, region=None): + """Dynamic Nearest Neighbor Search algorithm for estimating + + """ + pimg = img.select(["red","nir","swir1"]) + + kernel_size = 40 #pixels? + scale = pimg.projection().nominalScale() + + kernel = ee.Kernel.square(kernel_size, "pixels", False) + kernel_norm = ee.Kernel.square(kernel_size, 'pixels', True) + + water = img.select("mndwi").gt(0.15)#hf.bmax_otsu(img, band="mndwi", initial_threshold=0.1, region=region).Not() + land = img.select("mndwi").lt(-0.15) + + mix = water.Not().And(land.Not()) + ave_water = water.multiply(pimg).reduceRegion(ee.Reducer.mean(),region,scale,maxPixels=1e4,bestEffort=True).toImage() + + N_nmin_water = 1 + N_water = water.convolve(kernel) + + water_rf = ( + water.multiply(pimg) + .convolve(kernel) + .multiply(N_water.gte(N_nmin_water)) + .divide(N_water) + ) + water_rf = water_rf.add(ave_water.multiply(water_rf.Not())) + + ave_land = pimg.multiply(land).reduceRegion(ee.Reducer.mean(),region,scale,maxPixels=1e4,bestEffort=True).toImage() + + R1 = pimg.select("red").divide(pimg.select("swir1")) + R2 = pimg.select("nir").divide(pimg.select("swir1")) + R3 = R1.subtract(water_rf.select("red").divide(pimg.select("swir1"))) + R4 = R2.subtract(water_rf.select("red").divide(pimg.select("swir1"))) + + NR1 = R1.neighborhoodToBands(kernel) + NR2 = R2.neighborhoodToBands(kernel) + NI1 = img.select("red").neighborhoodToBands(kernel) + NI2 = img.select("nir").neighborhoodToBands(kernel) + NI3 = img.select("swir1").neighborhoodToBands(kernel) + + M1 = (NR1.gt(R3)).And(NR1.lt(R1)) + M2 = (NR2.gt(R4)).And(NR2.lt(R2)) + nLP = M1.And(M2) + + NnLP = nLP.reduce(ee.Reducer.sum()) + ave_nI1 = NI1.multiply(nLP).reduce(ee.Reducer.sum()).divide(NnLP) + ave_nI2 = NI2.multiply(nLP).reduce(ee.Reducer.sum()).divide(NnLP) + ave_nI3 = NI3.multiply(nLP).reduce(ee.Reducer.sum()).divide(NnLP) + + N_nmin_land = 1 + ave_pureland = ee.Image.cat([ave_nI1,ave_nI2,ave_nI3]) + ave_pureland = ave_pureland.multiply(NnLP.gte(N_nmin_land)).add( + ave_land.multiply(NnLP.lt(N_nmin_land)) + ) + + ave_landI3 = ave_land.select([2]) + f_water_all = ( + (ave_pureland.subtract(pimg)) + .divide(ave_pureland.subtract(water_rf)) + .clamp(0, 1) + ) + f_water = f_water_all.add(water).subtract(land).clamp(0, 1) + + return f_water.reduce(ee.Reducer.mean()).rename("water_fraction").toFloat() diff --git a/hydrafloods/tests/test_datasets.py b/hydrafloods/tests/test_datasets.py index d69a4e7..34d2bef 100644 --- a/hydrafloods/tests/test_datasets.py +++ b/hydrafloods/tests/test_datasets.py @@ -6,33 +6,50 @@ ee.Initialize() TEST_REGION = ee.Geometry.Rectangle([102.3335, 10.4045, 107.6277, 14.6900]) -TEST_START_TIME = "2019-02-07" -TEST_END_TIME = "2019-02-14" +# define two time periods to test for newer/older data +TEST_START_TIME_E2 = "2022-02-07" +TEST_END_TIME_E2 = "2022-02-14" + +TEST_START_TIME_E1 = "2009-02-07" +TEST_END_TIME_E1 = "2009-02-14" + class TestDatasets: def test_sentinel1(self): - s1 = hf.Sentinel1(TEST_REGION, TEST_START_TIME, TEST_END_TIME) - assert s1.n_images == 24 + s1 = hf.Sentinel1(TEST_REGION, TEST_START_TIME_E2, TEST_END_TIME_E2) + assert s1.n_images == 20 def test_sentinel2(self): - s2 = hf.Sentinel2(TEST_REGION, TEST_START_TIME, TEST_END_TIME) - assert s2.n_images == 101 + s2 = hf.Sentinel2(TEST_REGION, TEST_START_TIME_E2, TEST_END_TIME_E2) + assert s2.n_images == 100 + + def test_landsat9(self): + lc9 = hf.Landsat9(TEST_REGION, TEST_START_TIME_E2, TEST_END_TIME_E2) + assert lc9.n_images == 11 def test_landsat8(self): - lc8 = hf.Landsat8(TEST_REGION, TEST_START_TIME, TEST_END_TIME) - assert lc8.n_images == 11 + lc8 = hf.Landsat8(TEST_REGION, TEST_START_TIME_E2, TEST_END_TIME_E2) + assert lc8.n_images == 8 def test_landsat7(self): - lc7 = hf.Landsat7(TEST_REGION, TEST_START_TIME, TEST_END_TIME) - assert lc7.n_images == 7 + le7 = hf.Landsat7(TEST_REGION, TEST_START_TIME_E1, TEST_END_TIME_E1) + assert le7.n_images == 5 - def test_modis(self): - modis = hf.Modis(TEST_REGION, TEST_START_TIME, TEST_END_TIME) - assert modis.n_images == 7 + def test_landsat5(self): + lt5 = hf.Landsat5(TEST_REGION, TEST_START_TIME_E1, TEST_END_TIME_E1) + assert lt5.n_images == 7 + + def test_modis_terra(self): + modis_t = hf.Modis(TEST_REGION, TEST_START_TIME_E1, TEST_END_TIME_E1) + assert modis_t.n_images == 7 + + def test_modis_aqua(self): + modis_a = hf.Modis(TEST_REGION, TEST_START_TIME_E1, TEST_END_TIME_E1, asset_id="MODIS/006/MYD09GA") + assert modis_a.n_images == 7 def test_viirs(self): - viirs = hf.Viirs(TEST_REGION, TEST_START_TIME, TEST_END_TIME) + viirs = hf.Viirs(TEST_REGION, TEST_START_TIME_E2, TEST_END_TIME_E2) assert viirs.n_images == 7 def test_from_imgcollection(self): @@ -49,30 +66,31 @@ def test_from_imgcollection(self): assert planet.n_images == 33 def test_merge(self): - lc8 = hf.Landsat8(TEST_REGION, TEST_START_TIME, TEST_END_TIME) - s2 = hf.Sentinel2(TEST_REGION, TEST_START_TIME, TEST_END_TIME) + lc8 = hf.Landsat8(TEST_REGION, TEST_START_TIME_E2, TEST_END_TIME_E2) + viirs = hf.Viirs(TEST_REGION, TEST_START_TIME_E2, TEST_END_TIME_E2) - merged = lc8.merge(s2) - assert merged.n_images == 112 + merged = lc8.merge(viirs) + assert merged.n_images == 15 def test_join(self): - lc8 = hf.Landsat8(TEST_REGION, TEST_START_TIME, TEST_END_TIME) - s1 = hf.Sentinel1(TEST_REGION, TEST_START_TIME, TEST_END_TIME) + TEST_REGION_JOIN = ee.Geometry.Rectangle([0,0,120,30],"EPSG:4326", False) + lc8 = hf.Landsat8(TEST_REGION_JOIN, TEST_START_TIME_E2, TEST_END_TIME_E2) + s1 = hf.Sentinel1(TEST_REGION_JOIN, TEST_START_TIME_E2, TEST_END_TIME_E2) joined = lc8.join(s1) first = joined.collection.first() - assert (joined.n_images == 8) and ( + assert (joined.n_images == 110) and ( first.bandNames().getInfo() == ["blue", "green", "red", "nir", "swir1", "swir2", "VV", "VH", "angle"] ) def test_temporalagg(self): - terra = hf.Modis(TEST_REGION, TEST_START_TIME, TEST_END_TIME) + terra = hf.Modis(TEST_REGION, TEST_START_TIME_E1, TEST_END_TIME_E1) # get the aqua MODIS dataset # note calling the asset_id explicitly aqua = hf.Modis( - TEST_REGION, TEST_START_TIME, TEST_END_TIME, asset_id="MODIS/006/MYD09GA" + TEST_REGION, TEST_START_TIME_E1, TEST_END_TIME_E1, asset_id="MODIS/006/MYD09GA" ) merged = terra.merge(aqua) @@ -95,4 +113,3 @@ def test_temporalagg(self): "2018-01-01 00:00:00.000", "2019-01-01 00:00:00.000", ] - diff --git a/hydrafloods/workflows/dswfp.py b/hydrafloods/workflows/dswfp.py index 0e5e86f..beefb89 100644 --- a/hydrafloods/workflows/dswfp.py +++ b/hydrafloods/workflows/dswfp.py @@ -206,9 +206,9 @@ def export_fusion_samples( """ - dem = ee.Image("NASA/NASADEM_HGT/001").select("elevation") + dem = ee.Image("NASA/NASADEM_HGT/001").select("elevation").unmask(0) - optical_water_indices = ["mndwi", "nwi", "aewish", "aewinsh"] + optical_water_indices = ["mndwi", "nwi", "aewish", "aewinsh", "mbwi"] ds_kwargs = dict( region=region, start_time=start_time, end_time=end_time, rescale=True @@ -216,8 +216,8 @@ def export_fusion_samples( dsa_kwargs = {**ds_kwargs, **{"apply_band_adjustment": True}} lc8 = datasets.Landsat8(**ds_kwargs) - le7 = datasets.Landsat7(**dsa_kwargs) - s2 = datasets.Sentinel2(**dsa_kwargs) + le7 = datasets.Landsat7(**ds_kwargs) + s2 = datasets.Sentinel2(**ds_kwargs) _ = ds_kwargs.pop("rescale") @@ -235,7 +235,7 @@ def export_fusion_samples( corrections.slope_correction, dict( elevation=dem, - buffer=50, + buffer=100, ), ), hf.gamma_map, @@ -245,15 +245,15 @@ def export_fusion_samples( s1a.pipe(sar_proc, inplace=True) s1d.pipe(sar_proc, inplace=True) - s1a_anomalies = _calc_sar_anomalies(years, s1a) - s1d_anomalies = _calc_sar_anomalies(years, s1d) - - s1a.collection = s1a_anomalies - s1d.collection = s1d_anomalies + # s1a_anomalies = _calc_sar_anomalies(years, s1a) + # s1d_anomalies = _calc_sar_anomalies(years, s1d) + # + # s1a.collection = s1a_anomalies + # s1d.collection = s1d_anomalies s1 = s1a.merge(s1d) - optical = lc8.merge(s2).merge(le7) + optical = lc8.merge(le7) # .merge(s2) optical = optical.apply_func(geeutils.add_indices, indices=optical_water_indices) diff --git a/hydrafloods/workflows/hydrosar.py b/hydrafloods/workflows/hydrosar.py new file mode 100644 index 0000000..6b32851 --- /dev/null +++ b/hydrafloods/workflows/hydrosar.py @@ -0,0 +1,227 @@ +import os +import ee +import time +from ee.ee_exception import EEException +import gcsfs +import logging +import datetime +from functools import partial +import warnings +import hydrafloods as hf +from hydrafloods import ( + datasets, + timeseries, + filtering, + ml, + utils, + geeutils, + indices, + thresholding, + decorators, + corrections, + fuzzy, +) + + +# set global constant variables in module that are used throughout +MERIT = ee.Image("MERIT/Hydro/v1_0_1") +DEM = MERIT.select("elv") +HAND = MERIT.select("hnd") +HAND_THRESH = 15 # hand threshold +HAND_MASK = HAND.lt(HAND_THRESH) +SLOPE = ee.Terrain.slope(HAND.unmask(0)) +SAND = ee.Image("OpenLandMap/SOL/SOL_SAND-WFRACTION_USDA-3A1A1A_M/v02") + + +def hydro30( + region, + target_date, + speckle_filter="gamma_map", + grid_size=0.1, + scale=90, + combine_scenes=True, +): + + if not isinstance(target_date, ee.Date): + target_date = ee.Date(target_date) + + s1 = datasets.Sentinel1(region, target_date, target_date.advance(1, "day")) + + filter_opts = dict( + lee_sigma=filtering.lee_sigma, + gamma_map=filtering.gamma_map, + refined_lee=filtering.refined_lee, + ) + filter_func = filter_opts[speckle_filter] + + proc = ( + (corrections.slope_correction, dict(elevation=DEM, buffer=60)), + filter_func, + (_water_classification, dict(grid_size=grid_size, scale=scale)), + ) + + processed = s1.pipe(proc) + + if combine_scenes: + output = processed.collection.mosaic() + else: + output = processed.collection + + return output + + +@decorators.keep_attrs +def _water_classification(img, grid_size=0.1, scale=90): + def _per_band_water(img, tiles, bandname): + # !!! this function assumes one band image !!! + def _calc_t(tile): + histogram = band.updateMask(band.lt(0)).reduceRegion( + ee.Reducer.histogram(max_buckets, min_bucket_width, max_raw) + .combine("mean", None, True) + .combine("variance", None, True), + tile.geometry(), + scale=scale, + bestEffort=True, + tileScale=16, + ) + + threshold = hf.otsu(histogram.get(bandname.cat("_histogram"))) + + return tile.set("t", threshold) + + max_buckets = 255 + min_bucket_width = 0.001 + max_raw = 1e6 + bandname = ee.String(bandname) + band = img.select(bandname) + + thresh_tiles = tiles.map(_calc_t) + db_t = ee.Number(thresh_tiles.aggregate_array("t").reduce(ee.Reducer.mean())) + + water_init = band.lt(db_t) + + # db_stats = band.updateMask(band.lt(0)).reduceRegion( + # reducer=ee.Reducer.median(), + # geometry = region, + # scale=scale, + # bestEffort=True, + # tileScale=8 + # ) + # + # db_lower = db_stats.get(bandname) + # db_mid = db_t.add(db_lower).divide(2) + + fuzzy_db = fuzzy.fuzzy_large(band.select(bandname), db_t, 5) + + connected_objects = ( + water_init.select(bandname) + .selfMask() + .connectedComponents(connectedness=ee.Kernel.square(1.5), maxSize=1024) + ) + + connected_n = connected_objects.select("labels").connectedPixelCount( + maxSize=1024 + ) + + fuzzy_connected = fuzzy.fuzzy_zmf(connected_n.unmask(0), 10, 3) + + fuzziness = [ + fuzzy_slope, + fuzzy_hand, + fuzzy_db.select(bandname), + fuzzy_connected, + fuzzy_sand, + ] + + weights = [1 / len(fuzziness) for _ in range(len(fuzziness))] + + fuzzy_combo = hf.fuzzy_weighted(fuzziness, weights).updateMask( + band.select([0]).mask() + ) + + fuzzy_mask = fuzzy_combo.gt(0.60) + + water_final = water_init.And(fuzzy_mask) + + return water_final.rename(bandname.cat("_water")).uint8() + + def tile_stats(tile): + stats = img.reduceRegion( + reducer=stat_reducer, + geometry=tile.geometry(), + scale=scale, + bestEffort=True, + tileScale=16, + ) + + vv_denom = ee.Number(stats.get("VV_stdDev")) + vh_denom = ee.Number(stats.get("VH_stdDev")) + + vv_denom = ee.Algorithms.If(vv_denom, vv_denom, ee.Number(0.00001)) + vh_denom = ee.Algorithms.If(vh_denom, vh_denom, ee.Number(0.00001)) + + vv_ratio = ee.Number(img_stats.get("VV_p50")).divide(vv_denom) + vh_ratio = ee.Number(img_stats.get("VH_p50")).divide(vh_denom) + + stats = stats.set("VV_ratio", vv_ratio).set("VH_ratio", vh_ratio) + + return tile.set(stats) + + # + + db_bands = img.bandNames().removeAll(["angle", "local_inc_angle"]) + hd_bands = ["hnd", "hand_mask"] + img = img.select(db_bands).addBands(HAND_MASK).addBands(HAND) + + region = img.geometry() + + # this only creates tiles that are within the image so some tiles along edges will drop + tiles = hf.tile_region(region, contain_geom=region, grid_size=grid_size) + + stat_reducer = ee.Reducer.stdDev().combine(ee.Reducer.mean(), None, True) + + img_stats = img.reduceRegion( + reducer=ee.Reducer.percentile([50, 90]), + geometry=region, + scale=scale, + bestEffort=True, + tileScale=16, + ) + + stat_tiles = tiles.map(tile_stats).filter(ee.Filter.gte("hnd_mean", 0.8)) + + vv_sel_tiles = stat_tiles.limit(5, "VV_ratio", False) + vh_sel_tiles = stat_tiles.limit(5, "VH_ratio", False) + + # calculate the fuzzy membership of HAND Slope + fuzzy_slope = fuzzy.fuzzy_zmf(SLOPE, 0, 15) + + hand_stats = HAND.updateMask( + HAND.neq(0).And(HAND.lt(ee.Number(img_stats.get("hnd_p90")))) + ).reduceRegion( + reducer=ee.Reducer.median().combine(ee.Reducer.stdDev(), None, True), + geometry=region, + scale=scale, + bestEffort=True, + tileScale=16, + ) + + hand_lower = ee.Number(hand_stats.get("hnd_median")) + hand_std = ee.Number(hand_stats.get("hnd_stdDev")) + hand_upper = hand_lower.add(ee.Number(3.0).multiply(hand_std)).add(5.0) + + fuzzy_hand = fuzzy.fuzzy_zmf(HAND, hand_lower, hand_upper).unmask(1) + + fuzzy_sand = fuzzy.fuzzy_small(SAND.select([0]).unmask(0), 37.5, 8) + + vv_water = _per_band_water(img, vv_sel_tiles, "VV") + vh_water = _per_band_water(img, vh_sel_tiles, "VH") + + overlap = vv_water.Or(vh_water) + + # post processing opening of the segmented water + opened_closed_water = filtering.close_binary(filtering.open_binary(overlap)) + + out_water = opened_closed_water.rename("water").updateMask(img.select([0]).mask()) + + return out_water From 87eb9757c1ff9e3f3f23a08e1a0f56c8bc75ca27 Mon Sep 17 00:00:00 2001 From: Kel Markert Date: Mon, 20 Mar 2023 23:46:53 -0500 Subject: [PATCH 2/7] changing s2 collection to harmonized version and adding method to deduplicate images in collection due to repeats from processing versions --- hydrafloods/datasets.py | 45 +++++++++++++++++++++++++++++- hydrafloods/tests/test_datasets.py | 5 ++++ 2 files changed, 49 insertions(+), 1 deletion(-) diff --git a/hydrafloods/datasets.py b/hydrafloods/datasets.py index 7865b80..5951eb1 100644 --- a/hydrafloods/datasets.py +++ b/hydrafloods/datasets.py @@ -974,7 +974,7 @@ class Sentinel2(Dataset): def __init__( self, *args, - asset_id="COPERNICUS/S2_SR", + asset_id="COPERNICUS/S2_SR_HARMONIZED", use_qa=True, apply_band_adjustment=False, **kwargs, @@ -1013,6 +1013,49 @@ def __init__( return + + def deduplicate(self, inplace=False): + """Deduplicate Sentinel2 images by removing images with the same data take identifier""" + def _dedupe(id): + """Helper function to remove duplicate images""" + id_collection = ( + self.collection + .filter( + ee.Filter.stringStartsWith("DATATAKE_IDENTIFIER", ee.String(id)) + ) + .sort('PROCESSING_BASELINE') + ) + image = ( + id_collection + .mosaic() + ) + + imgbounds = id_collection.map(lambda x: x.geometry(1e3)) + + start_time = id_collection.aggregate_array('system:time_start').reduce(ee.Reducer.min()) + end_time = id_collection.aggregate_array('system:time_start').reduce(ee.Reducer.max()) + + return image.clipToCollection(imgbounds).set( + { + 'system:time_start': start_time, + 'system:time_end':end_time + } + ) + + datatakes = ( + self.collection + .aggregate_histogram("DATATAKE_IDENTIFIER") + .keys() + .map(lambda x: ee.String(x).slice(0,-7)) + ) + + images = ee.ImageCollection.fromImages( + datatakes.map(_dedupe) + ) + + return self._inplace_wrapper(images, inplace) + + @decorators.keep_attrs def qa(self, img): """Custom QA masking method for Sentinel2 surface reflectance dataset""" diff --git a/hydrafloods/tests/test_datasets.py b/hydrafloods/tests/test_datasets.py index 34d2bef..e303d64 100644 --- a/hydrafloods/tests/test_datasets.py +++ b/hydrafloods/tests/test_datasets.py @@ -24,6 +24,11 @@ def test_sentinel2(self): s2 = hf.Sentinel2(TEST_REGION, TEST_START_TIME_E2, TEST_END_TIME_E2) assert s2.n_images == 100 + def test_sentinel2_dedupe(self): + s2 = hf.Sentinel2(TEST_REGION, TEST_START_TIME_E2, TEST_END_TIME_E2) + s2_dedupe = s2.deduplicate() + assert s2_dedupe.n_images == 5 + def test_landsat9(self): lc9 = hf.Landsat9(TEST_REGION, TEST_START_TIME_E2, TEST_END_TIME_E2) assert lc9.n_images == 11 From 66e426f5a1dc3300388e4365c6dba0592efc490c Mon Sep 17 00:00:00 2001 From: Kel Markert Date: Mon, 20 Mar 2023 23:52:56 -0500 Subject: [PATCH 3/7] updating expected result for join --- hydrafloods/decorators.py | 6 +----- hydrafloods/tests/test_datasets.py | 2 +- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/hydrafloods/decorators.py b/hydrafloods/decorators.py index a292fb8..73817dc 100644 --- a/hydrafloods/decorators.py +++ b/hydrafloods/decorators.py @@ -26,11 +26,7 @@ def wrapper(*args, **kwargs): # this assumption is true for 99% of functions used for ee.ImageCollection.map() result = ee.Image(func(*args, **kwargs)) img = [i for i in args if isinstance(i, ee.Image)][0] - return ee.Image( - result.copyProperties(img).set( - "system:time_start", img.get("system:time_start") - ) - ) + return img.select().addBands(result) return wrapper diff --git a/hydrafloods/tests/test_datasets.py b/hydrafloods/tests/test_datasets.py index e303d64..206b8aa 100644 --- a/hydrafloods/tests/test_datasets.py +++ b/hydrafloods/tests/test_datasets.py @@ -85,7 +85,7 @@ def test_join(self): joined = lc8.join(s1) first = joined.collection.first() - assert (joined.n_images == 110) and ( + assert (joined.n_images == 105) and ( first.bandNames().getInfo() == ["blue", "green", "red", "nir", "swir1", "swir2", "VV", "VH", "angle"] ) From d4964e18b6b85099812b4e0189c389f9580e85f5 Mon Sep 17 00:00:00 2001 From: Kel Markert Date: Mon, 20 Mar 2023 23:58:14 -0500 Subject: [PATCH 4/7] setting system:index as well when deduping s2 collection --- hydrafloods/datasets.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/hydrafloods/datasets.py b/hydrafloods/datasets.py index 5951eb1..8d19617 100644 --- a/hydrafloods/datasets.py +++ b/hydrafloods/datasets.py @@ -1038,7 +1038,8 @@ def _dedupe(id): return image.clipToCollection(imgbounds).set( { 'system:time_start': start_time, - 'system:time_end':end_time + 'system:time_end':end_time, + 'system:index':id } ) From db29dca3aa46ce62b05b1cb6225e95e29034c682 Mon Sep 17 00:00:00 2001 From: Kel Markert Date: Tue, 21 Mar 2023 00:51:19 -0500 Subject: [PATCH 5/7] updating qa function to include flag for poor masked data --- hydrafloods/datasets.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/hydrafloods/datasets.py b/hydrafloods/datasets.py index 8d19617..486ab2e 100644 --- a/hydrafloods/datasets.py +++ b/hydrafloods/datasets.py @@ -997,6 +997,9 @@ def __init__( self.BANDREMAP.get("sen2"), self.BANDREMAP.get("new") ) + if use_qa: + coll = coll.filter(ee.Filter.eq('cloud_mask_success', True)) + if apply_band_adjustment: # band bass adjustment coefficients taken HLS project https://hls.gsfc.nasa.gov/algorithms/bandpass-adjustment/ # slope coefficients @@ -1069,10 +1072,14 @@ def qa(self, img): # Get s2cloudless image, subset the probability band. cld_prb = ee.Image( ee.ImageCollection("COPERNICUS/S2_CLOUD_PROBABILITY") + .filterDate(self.start_time, self.end_time) + .filterBounds(self.region) .filter(ee.Filter.eq("system:index", img.get("system:index"))) .first() ).select("probability") + worked = cld_prb.bandNames().length().gt(0) + # Condition s2cloudless by the probability threshold value. is_cloud = cld_prb.gt(CLD_PRB_THRESH) @@ -1111,4 +1118,4 @@ def qa(self, img): ) # Subset reflectance bands and update their masks, return the result. - return geeutils.rescale(img).select("B.*").updateMask(is_cld_shdw.Not()) + return geeutils.rescale(img).select("B.*").updateMask(is_cld_shdw.Not()).set({'cloud_mask_success':worked}) From 4bf64cedcc09bb4d3c9af3352a0ec1a4cc6e2a91 Mon Sep 17 00:00:00 2001 From: Kel Markert Date: Tue, 21 Mar 2023 00:52:40 -0500 Subject: [PATCH 6/7] updating qa function to include flag for poor masked data --- hydrafloods/datasets.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hydrafloods/datasets.py b/hydrafloods/datasets.py index 486ab2e..f8dc7d9 100644 --- a/hydrafloods/datasets.py +++ b/hydrafloods/datasets.py @@ -998,7 +998,7 @@ def __init__( ) if use_qa: - coll = coll.filter(ee.Filter.eq('cloud_mask_success', True)) + coll = coll.filter(ee.Filter.eq('cloud_mask_success', 1)) if apply_band_adjustment: # band bass adjustment coefficients taken HLS project https://hls.gsfc.nasa.gov/algorithms/bandpass-adjustment/ From 996f24514af9564e9b2f61a39974e269f5a62606 Mon Sep 17 00:00:00 2001 From: Kel Markert Date: Sat, 14 Oct 2023 15:09:53 -0600 Subject: [PATCH 7/7] version bump --> 2023.10.14 --- hydrafloods/__init__.py | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/hydrafloods/__init__.py b/hydrafloods/__init__.py index 947f83e..a7f4844 100644 --- a/hydrafloods/__init__.py +++ b/hydrafloods/__init__.py @@ -13,4 +13,4 @@ # from hydrafloods import * -__version__ = "2021.11.10" +__version__ = "2023.10.14" diff --git a/setup.py b/setup.py index 32a3134..9cde2d6 100644 --- a/setup.py +++ b/setup.py @@ -6,7 +6,7 @@ setup( name="hydrafloods", - version="2021.11.10", + version="2023.10.14", description="HYDrologic Remote sensing Analysis for Floods", long_description=long_description, long_description_content_type="text/markdown",