I am working with Google Earth Engine and the wonderful geemap package to export channel centerlines using Modified Normalized Difference Water Index (MNDWI; Xu, 2006). Computing the MNDWI is simple enough:

``````import os
import ee
import geemap
import matplotlib.pyplot as plt
import geemap.colormaps as cm
import sys
os.environ["PLANET_API_KEY"] = XXXXXX
ee.Initialize()
geojson = ... #omitted from question due to length
Map = geemap.Map()

collection = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2').filterDate('2018-01-01', '2020-12-31').filterBounds(geemap.geojson_to_ee(geojson)).filterMetadata('CLOUD_COVER', 'less_than', 10).sort("CLOUD_COVER")

im = collection.first()

B3 = im.select('SR_B3')
B6 = im.select('SR_B6')
im = im.clip(geemap.geojson_to_ee(geojson))
mndwi = im.normalizedDifference(['SR_B3', 'SR_B6'])
palette = list(cm.palettes.hsv.n08)
Map

``````

I then Otsu threshold the image to binarize land vs. water. This works well also in the context of Google Earth Engine and geemap. Code is from GEE documentation and Dr. Qiusheng Wu's geemap tutorials.

``````histogram = mndwi.reduceRegion(**{
'reducer': ee.Reducer.histogram(255),
'geometry': geemap.geojson_to_ee(geojson),
'scale': 10,
'bestEffort': True
})
hist_dict = histogram.getInfo()

# 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)
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('nd'))
print('threshold', threshold.getInfo())

classA = mndwi.select('nd').lt(threshold)
``````

This looks great, the mask covers everything except the channel water. Exporting this becomes a sticky issue, however...

When I load the exported masked raster into QGIS, the entire image is black with values of 1.

Here is the exporting code

``````channel = classA.mask(classA.Not())
geemap.ee_export_image(channel, 'channel.tif')
``````

I understand that the mask and MNDWI images are separate layers in Google Earth Engine, but how do I binarize the masked image so that a value of 1 = channel water and export it to a local file?

• Can you please add the code with which you exported the image?
– Matt
Commented Nov 1, 2021 at 12:51
• sure thing, just did! Commented Nov 1, 2021 at 16:02

I was able to download a water mask by selecting pixels greater than the threshold value, using `updateMask` rather than `mask`, and adding the `scale` and `region` parameters in the `geemap.export_image()` function.

Furthermore, I recommend using a variable for your region to save having to call the conversion function each time.

``````roi = geemap.geojson_to_ee(geojson)

classA = mndwi.select('nd').gt(threshold)  ## invert the selection
channel = classA.updateMask(classA)  ## omit `.Not()` here
geemap.ee_export_image(channel, 'channel.tif', scale=30, region=roi.geometry())
``````

Result (using an arbitrary GeoJSON polygon):

• Thanks, Matt. It seems weird that I can't get the inverse of the mask to export. Any suggestions? otherwise I can mark this as correct. Commented Nov 2, 2021 at 18:21
• My apologies, I thought the issue was just the uniform value image. To get only the channel as the mask, you can swap you `lt` (less than) for `gt` (greater than). I have updated my answer
– Matt
Commented Nov 2, 2021 at 18:46
• perfect! this is the correct answer. Commented Nov 2, 2021 at 18:47