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
geojson = ... #omitted from question due to length
Map = geemap.Map()
# Add Earth Engine dataset

Map.add_planet_by_month(year=2020, month=8)
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.addLayer(mndwi, {'palette': palette}, "MNDWI")

MNDWI of Bolivian River from LANDSAT-OLI

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]) \
        bCount = total.subtract(aCount)
        bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount)
        return aCount.multiply(aMean.subtract(mean).pow(2)).add(

    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)
Map.addLayer(classA.mask(classA), {'palette': 'blue'}, 'class A')

Binarized Mask based on Otsu 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')

Exported Mask

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
    Nov 1, 2021 at 12:51
  • sure thing, just did!
    – jameshgrn
    Nov 1, 2021 at 16:02

1 Answer 1


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): enter image description here

  • 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.
    – jameshgrn
    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
    Nov 2, 2021 at 18:46
  • 1
    perfect! this is the correct answer.
    – jameshgrn
    Nov 2, 2021 at 18:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.