How to export masked MNDWI raster from Google Earth Engine

Question

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

classA = mndwi.select('nd').lt(threshold)
Map.addLayer(classA.mask(classA), {'palette': 'blue'}, 'class A')

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?

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):