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