Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

add water fill confidence margin #15

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open

add water fill confidence margin #15

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0c2ad17
add water fill confidence margin
Jaapel Aug 25, 2022
955c601
formatting using black
Jaapel Aug 25, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
363 changes: 212 additions & 151 deletions eepackages/applications/waterbody_area.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,161 +2,222 @@
from eepackages import assets
from eepackages import utils

def computeSurfaceWaterArea(waterbody, start_filter, start, stop, scale, waterOccurrence, opt_missions):
geom = ee.Feature(waterbody).geometry()

missions = ['L4', 'L5', 'L7', 'L8', 'S2']

if opt_missions:
missions = opt_missions

images = assets.getImages(geom, {
'resample': True,
'filter': ee.Filter.date(start_filter, stop),
'missions': missions,
'scale': scale * 10
})

# print('Image count: ', images.size())

options = {
# 'cloudFrequencyThresholdDelta': -0.15
'scale': scale * 5
}

g = geom.buffer(300, scale)

# images = assets.getMostlyCleanImages(images, g, options).sort('system:time_start')
# images = images.filterDate(start, stop)

# all images, but with quality score
images = (assets.addQualityScore(images, g, options)
.filter(ee.Filter.gt('quality_score', 0))) # sometimes null?!


# print('Image count (clean): ', images.size())

water = images.map(lambda i: computeSurfaceWaterArea_SingleImage(i, waterbody, scale, waterOccurrence))

water = water.filter(ee.Filter.neq('area', 0))

return water

def computeSurfaceWaterArea_SingleImage(i, waterbody, scale, waterOccurrence):
geom = ee.Feature(waterbody).geometry()

fillPercentile = 50 # // we don't trust our prior

ndwiBands = ['green', 'swir']
# var ndwiBands = ['green', 'nir']

waterMaxImage = ee.Image().float().paint(waterbody.buffer(150), 1)

maxArea = waterbody.area(scale)

t = i.get('system:time_start')

i = (i
.updateMask(waterMaxImage)
.updateMask(i.select('swir').min(i.select('nir')).gt(0.001)))

ndwi = i.normalizedDifference(ndwiBands)

# var water = ndwi.gt(0)

th = utils.computeThresholdUsingOtsu(ndwi, scale, geom, 0.5, 0.7, -0.2)
water = ndwi.gt(th)

area = (ee.Image.pixelArea().mask(water)
.reduceRegion(**{
'reducer': ee.Reducer.sum(),
'geometry': geom,
'scale': scale
}).get('area'))

# # fill missing, estimate water probability as the lowest percentile.
# var waterEdge = ee.Algorithms.CannyEdgeDetector(water, 0.1, 0)
waterEdge = ee.Algorithms.CannyEdgeDetector(ndwi, 0.5, 0.7)

# image mask
imageMask = ndwi.mask()

# var imageMask = i.select(ndwiBands).reduce(ee.Reducer.allNonZero())

# imageMask = utils.focalMin(imageMask, ee.Number(scale) * 1.5)
imageMask = imageMask.focal_min(ee.Number(scale).multiply(1.5), 'square', 'meters')
#.multiply(waterMaxImage)

# TODO: exclude non-water/missing boundsry
waterEdge = waterEdge.updateMask(imageMask)

# # clip by clouds
# var bad = ee.Image(1).float().subtract(i.select('weight'))
# bad = utils.focalMax(bad, 90).not()
# waterEdge = waterEdge.updateMask(bad)

# get water probability around edges
# P(D|W) = P(D|W) * P(W) / P(D) ~= P(D|W) * P(W)
p = waterOccurrence.mask(waterEdge).reduceRegion(**{
'reducer': ee.Reducer.percentile([fillPercentile]),
'geometry': geom,
'scale': scale
}).values().get(0)

# TODO: exclude edges belonging to cloud/water or cloud/land

# TODO: use multiple percentiles (confidence margin)

p = ee.Algorithms.If(ee.Algorithms.IsEqual(p, None), 101, p)

waterFill = (waterOccurrence.gt(ee.Image.constant(p))
.updateMask(water.unmask(0, False).Not()))

# exclude false-positive, where we're sure in a non-water
nonWater = ndwi.lt(-0.15).unmask(0, False)
waterFill = waterFill.updateMask(nonWater.Not())

fill = (ee.Image.pixelArea().mask(waterFill)
.reduceRegion(**{
'reducer': ee.Reducer.sum(),
'geometry': geom,
'scale': scale
}).get('area'))

area_filled = ee.Number(area).add(fill)

filled_fraction = ee.Number(fill).divide(area_filled)

return (i
.addBands(waterFill.rename('water_fill'))
.addBands(waterEdge.rename('water_edge'))
.addBands(ndwi.rename('ndwi'))
.addBands(water.rename('water'))
.set({
'p': p,
'area': area,
'area_filled': area_filled,
'filled_fraction': filled_fraction,
'system:time_start': t,
'ndwi_threshold': th,
'waterOccurrenceExpected': waterOccurrence.mask(waterEdge)
}))

def computeSurfaceWaterArea(
waterbody, start_filter, start, stop, scale, waterOccurrence, opt_missions
):
geom = ee.Feature(waterbody).geometry()

missions = ["L4", "L5", "L7", "L8", "S2"]

if opt_missions:
missions = opt_missions

images = assets.getImages(
geom,
{
"resample": True,
"filter": ee.Filter.date(start_filter, stop),
"missions": missions,
"scale": scale * 10,
},
)

# print('Image count: ', images.size())

options = {
# 'cloudFrequencyThresholdDelta': -0.15
"scale": scale
* 5
}

g = geom.buffer(300, scale)

# images = assets.getMostlyCleanImages(images, g, options).sort('system:time_start')
# images = images.filterDate(start, stop)

# all images, but with quality score
images = assets.addQualityScore(images, g, options).filter(
ee.Filter.gt("quality_score", 0)
) # sometimes null?!

# print('Image count (clean): ', images.size())

water = images.map(
lambda i: computeSurfaceWaterArea_SingleImage(
i, waterbody, scale, waterOccurrence
)
)

water = water.filter(ee.Filter.neq("area", 0))

return water


def computeSurfaceWaterArea_SingleImage(
i, waterbody, scale, waterOccurrence, fillPercentiles=None
):
geom = ee.Feature(waterbody).geometry()
if fillPercentiles is None:
fillPercentiles = [5, 50, 95] # we don't trust our prior
fillPercentilesStr = list(map(lambda n: str(n), fillPercentiles))

ndwiBands = ["green", "swir"]
# var ndwiBands = ['green', 'nir']

waterMaxImage = ee.Image().float().paint(waterbody.buffer(150), 1)

# maxArea = waterbody.area(scale)

t = i.get("system:time_start")

i = i.updateMask(waterMaxImage).updateMask(
i.select("swir").min(i.select("nir")).gt(0.001)
)

ndwi = i.normalizedDifference(ndwiBands)

# var water = ndwi.gt(0)

th = utils.computeThresholdUsingOtsu(ndwi, scale, geom, 0.5, 0.7, -0.2)
water = ndwi.gt(th)

area = (
ee.Image.pixelArea()
.mask(water)
.reduceRegion(**{"reducer": ee.Reducer.sum(), "geometry": geom, "scale": scale})
.get("area")
)

# # fill missing, estimate water probability as the lowest percentile.
# var waterEdge = ee.Algorithms.CannyEdgeDetector(water, 0.1, 0)
waterEdge = ee.Algorithms.CannyEdgeDetector(ndwi, 0.5, 0.7)

# image mask
imageMask = ndwi.mask()

# var imageMask = i.select(ndwiBands).reduce(ee.Reducer.allNonZero())

# imageMask = utils.focalMin(imageMask, ee.Number(scale) * 1.5)
imageMask = imageMask.focal_min(ee.Number(scale).multiply(1.5), "square", "meters")
# .multiply(waterMaxImage)

# TODO: exclude non-water/missing boundsry
waterEdge = waterEdge.updateMask(imageMask)

# # clip by clouds
# var bad = ee.Image(1).float().subtract(i.select('weight'))
# bad = utils.focalMax(bad, 90).not()
# waterEdge = waterEdge.updateMask(bad)

# get water probability around edges
# P(D|W) = P(D|W) * P(W) / P(D) ~= P(D|W) * P(W)
pValues = (
waterOccurrence.mask(waterEdge)
.reduceRegion(
**{
"reducer": ee.Reducer.percentile(fillPercentiles),
"geometry": geom,
"scale": scale,
}
)
.values()
)

# TODO: exclude edges belonging to cloud/water or cloud/land

pValues = pValues.map(
lambda p: ee.Algorithms.If(ee.Algorithms.IsEqual(p, None), 101, p)
)
pValuesDict = ee.Dictionary(
ee.List(fillPercentilesStr)
.map(lambda s: ee.String("p").cat(ee.String(s)))
.zip(pValues)
.flatten()
)
i = i.set(pValuesDict)

waterFills = pValues.map(
lambda p: (
waterOccurrence.gt(ee.Image.constant(p)).updateMask(
water.unmask(0, False).Not()
)
)
)

# exclude false-positive, where we're sure in a non-water
nonWater = ndwi.lt(-0.15).unmask(0, False)
# waterFill = waterFill.updateMask(nonWater.Not())
waterFills = waterFills.map(lambda wf: ee.Image(wf).updateMask(nonWater.Not()))

# create an image with water fills as bands
waterFillsBands = (
ee.ImageCollection(waterFills)
.toBands()
.rename(list(map(lambda s: f"water_fill_{s}", fillPercentilesStr)))
)

fills = waterFills.map(
lambda wf: (
ee.Image.pixelArea()
.mask(wf)
.reduceRegion(
**{"reducer": ee.Reducer.sum(), "geometry": geom, "scale": scale}
)
.get("area")
)
)

areasFilled = fills.map(lambda f: ee.Number(area).add(f))
areasFilledDict = ee.Dictionary(
ee.List(fillPercentilesStr)
.map(lambda s: ee.String("area_filled_p").cat(s))
.zip(areasFilled)
.flatten()
)
i = i.set(areasFilledDict)

filledFractions = ee.Array(fills).divide(ee.Array(areasFilled)).toList()
filledFractionsDict = ee.Dictionary(
ee.List(fillPercentilesStr)
.map(lambda s: ee.String("filled_fraction_p").cat(s))
.zip(ee.List(filledFractions))
.flatten()
)
i = i.set(filledFractionsDict)

return (
i.addBands(waterEdge.rename("water_edge"))
.addBands(ndwi.rename("ndwi"))
.addBands(water.rename("water"))
.addBands(waterFillsBands)
.set(
{
"area": area,
"system:time_start": t,
"ndwi_threshold": th,
"waterOccurrenceExpected": waterOccurrence.mask(waterEdge),
}
)
)


def computeSurfaceWaterAreaJRC(waterbody, start, stop, scale):
geom = ee.Feature(waterbody).geometry()

jrcMonthly = ee.ImageCollection("JRC/GSW1_0/MonthlyHistory")
geom = ee.Feature(waterbody).geometry()

def compute_area(i):
area = ee.Image.pixelArea().mask(i.eq(2)).reduceRegion({
'reducer': ee.Reducer.sum(),
'geometry': geom,
'scale': scale
})
jrcMonthly = ee.ImageCollection("JRC/GSW1_0/MonthlyHistory")

return i.set({'area': area.get('area')})
def compute_area(i):
area = (
ee.Image.pixelArea()
.mask(i.eq(2))
.reduceRegion(
{"reducer": ee.Reducer.sum(), "geometry": geom, "scale": scale}
)
)

water = jrcMonthly.filterDate(start, stop).map(compute_area)
return i.set({"area": area.get("area")})

return water
water = jrcMonthly.filterDate(start, stop).map(compute_area)

return water