Google Earth Engine – How to Get Class with Highest Probability from Bands

Question

I'm trying to use Google's Dynamic World classified rasters to get a land/water map for another tool I've built. I'm having trouble doing so.
So far, I have essentially copied this code (in Python syntax) into a notebook, which is working well. For completeness.. It looks like this:

import ee
import geemap

Map = geemap.Map()
Map

COL_FILTER = ee.Filter.And(
    ee.Filter.bounds(ee.Geometry.Point(-76.409826, 36.898152)),
    ee.Filter.date('2021-04-02', '2022-04-03'))
dwCol = ee.ImageCollection("GOOGLE/DYNAMICWORLD/V1").filter(COL_FILTER)
dwImage = ee.Image(dwCol.first())
CLASS_NAMES = [
    'water', 'trees', 'grass', 'flooded_vegetation', 'crops',
    'shrub_and_scrub', 'built', 'bare', 'snow_and_ice'];
VIS_PALETTE = [
    '#419BDF', '#397D49', '#88B053', '#7A87C6',
    '#E49635', '#DFC35A', '#C4281B', '#A59B8F',
    '#B39FE1'];
dwRgb = dwImage.select('label').visualize(min=0, max=8, palette=VIS_PALETTE).divide(255)
top1Prob = dwImage.select(CLASS_NAMES).reduce(ee.Reducer.max());
top1ProbHillshade = ee.Terrain.hillshade(top1Prob.multiply(100)).divide(255);
dwRgbHillshade = dwRgb.multiply(top1ProbHillshade);
Map.setCenter(-76.409826, 36.898152, 12);
Map.addLayer(dwRgbHillshade, {min: 0, max: 0.65}, 'Dynamic World')

which creates a map that looks like this:

I can tell visually that this is correct, where the water is blue, land is something else. The problem is when I try to get land/water maps out of this. I've been struggling to understand GEE's APIs in a way that can do this. I was looking at reducer, but it doesn't have much documentation on how I can reduce the bands into a single classification. Right now the shape is three dimensional, where the third dimension, the band (axis=2), all have values. At least.. it is for dwImage, which is what I was hoping to reduce on. I could not find any examples that created a custom reducing function to show how to do this in GEE.

I got desperate and tried writing this to a raster file which worked. Using the rectangular region I drew on the map above..

feature = Map.draw_last_feature
roi = feature.geometry()
clipped = dwImage.clip(roi).unmask()
out_dir = os.path.join(os.path.expanduser('~'), 'Downloads')
filename = os.path.join(out_dir, 'dwclasses.tif')
clipped = dwImage.clip(roi).unmask()
geemap.ee_export_image(
    clipped, filename=filename, scale=20, region=roi, file_per_band=False
)
# now to process it...
dataset = gdal.Open(filename, gdal.GA_ReadOnly)
bands = [
    dataset.GetRasterBand(k + 1).ReadAsArray() for k in range(dataset.RasterCount)
]
stacked = np.dstack(bands)[:,:,:9] # -1 to remove alpha band
band_with_highest = np.amax(stacked, axis=2)
water_land = band_with_highest == stacked[:, :, 1] # where the highest probability equalled water, set as True, else it is False and is land
mask = np.where(np.amax(stacked, axis=2) == np.zeros(stacked.shape[:2]), np.nan, 1)
fig, ax = plt.subplots(figsize=(15,15))
ax.imshow(water_land* mask)

Which outputs… (note, its inverted vertically)

but this is visually wrong too… I tried outputting the dwRgb and top1Prob rasters out but they weren't too helpful from what I can tell. The closest I got was by thresholding the water probability like..

fig, ax = plt.subplots(figsize=(15,15))
ax.imshow(stacked[:,:,1] > .05 * mask)

While that looks fine, its an unreliable cheat…

Answer 1

When working with a temporal stack of Dynamic World images, it is better to perform your analysis on the probability bands, and only convert to a categorical "label" as a last step (see top_probability in the following code example).

(By contrast, aggregating the categorical "label" band of individual images can produce noisy results if several classes have similar, high probabilities.)

# Copyright 2022 Google LLC.
# SPDX-License-Identifier: Apache-2.0

dw = ee.ImageCollection('GOOGLE/DYNAMICWORLD/V1')

probability_bands = [
  'water', 'trees', 'grass', 'flooded_vegetation', 'crops',
  'shrub_and_scrub', 'built', 'bare', 'snow_and_ice',
]
palette = [
  '#419BDF', '#397D49', '#88B053', '#7A87C6', '#E49635', 
  '#DFC35A', '#C4281B', '#A59B8F', '#B39FE1'
]

start_date = '2019-04-01'
end_date = '2019-07-01'

# Filter the Dynamic World image collection by time
dw_time_interval = dw.filter(ee.Filter.date(start_date, end_date))

# Select probability bands 
dw_time_series = dw_time_interval.select(probability_bands)

# Create a multi-band image summarizing probability 
# for each band across the time-period
mean_probability = dw_time_series.reduce(ee.Reducer.mean())

# Create a single band image containing the class with the top probability
top_probability = mean_probability.toArray().arrayArgmax().arrayGet(0).rename('label')

Link to full code example: https://colab.research.google.com/gist/tylere/40aa0ea2f67653959dbef6b05a11b7c9/dynamic-world-top-probability.ipynb#scrollTo=rmrPAJoYbWZr