import ee
import geemap
from bqplot import pyplot as plt
from bqplot import Bars
Map = geemap.Map()
Map
roi = ee.Geometry.Point([-99.09805297851562, 47.10016154728056])
# set visualization parameters
vis = {
'nir': {'bands': ['N', 'R', 'G']},
'rgb': {'bands': ['R', 'G', 'B']},
'water': {'palette': ['0000FF']},
'ndwi': {
'palette': [
'#ece7f2',
'#d0d1e6',
'#a6bddb',
'#74a9cf',
'#3690c0',
'#0570b0',
'#045a8d',
'#023858',
]
},
'ndvi': {
'palette': [
'#d73027',
'#f46d43',
'#fdae61',
'#fee08b',
'#d9ef8b',
'#a6d96a',
'#66bd63',
'#1a9850',
]
},
}
# search NAIP imagery that has RGBN bands
collection = (
ee.ImageCollection('USDA/NAIP/DOQQ')
.filterBounds(roi)
.filterDate('2009-01-01', '2019-12-31')
.filter(ee.Filter.listContains("system:band_names", "N"))
)
# print(collection.getInfo())
image = collection.first()
Map.centerObject(image)
Map.addLayer(image, vis['nir'], 'NAIP')
polygon = ee.Geometry(image.geometry())
# Compute the histogram of the NIR band. The mean and variance are only FYI.
histogram = image.select('N').reduceRegion(
**{
'reducer': ee.Reducer.histogram(255, 2),
'geometry': polygon,
'scale': 10,
'bestEffort': True,
}
)
hist_dict = histogram.getInfo()
# print(hist_dict)
x = hist_dict['N']['bucketMeans']
y = hist_dict['N']['histogram']
plt.bar(x, y)
plt.show()
# Return the DN that maximizes interclass variance in B5 (in the region).
def otsu(histogram):
counts = ee.Array(ee.Dictionary(histogram).get('histogram'))
means = ee.Array(ee.Dictionary(histogram).get('bucketMeans'))
size = means.length().get([0])
total = counts.reduce(ee.Reducer.sum(), [0]).get([0])
sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0])
mean = sum.divide(total)
indices = ee.List.sequence(1, size)
# Compute between sum of squares, where each mean partitions the data.
def func_xxx(i):
aCounts = counts.slice(0, 0, i)
aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0])
aMeans = means.slice(0, 0, i)
aMean = (
aMeans.multiply(aCounts)
.reduce(ee.Reducer.sum(), [0])
.get([0])
.divide(aCount)
)
bCount = total.subtract(aCount)
bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount)
return aCount.multiply(aMean.subtract(mean).pow(2)).add(
bCount.multiply(bMean.subtract(mean).pow(2))
)
bss = indices.map(func_xxx)
# Return the mean value corresponding to the maximum BSS.
return means.sort(bss).get([-1])
threshold = otsu(histogram.get('N'))
print('threshold', threshold.getInfo())
classA = image.select('N').lt(threshold)
Map.addLayer(classA.mask(classA), {'palette': 'blue'}, 'class A')
Map
def extract_water(image):
histogram = image.select('N').reduceRegion(
**{
'reducer': ee.Reducer.histogram(255, 2),
'geometry': polygon,
'scale': 10,
'bestEffort': True,
}
)
threshold = otsu(histogram.get('N'))
water = image.select('N').lt(threshold).selfMask()
return water.set({"threshold": threshold})
water_images = collection.map(extract_water)
Map.addLayer(water_images.first(), {"palette": "blue"}, "first water")
water_images.aggregate_array("threshold").getInfo()
Map