1

I am working on a project to analyze land area for a coastal area with several nearshore islands. To approach this I am using Landsat 8, creating a land/water mask based on an NDWI threshold, and reclassifying to either land or water (only the first two steps are in the code pasted below as this is the focus of my question). Google has provided helpful documentation for this, although using a single image, not an image collection. Similarly, this post describes a similar question applied to MODIS, but the NDWI.map() function does not work in my example.

Is there a method or best practice to either (1) map the NDWI function across an image collection (done already) and classify this NDWI variable to land/water, or (2) convert the image collection subset to a single mosaicked image, then calculate the NDWI and classify using the mosaicked image as an input? I assume option 1 would be more efficient, but the last code block below notes that “ndwi.map is not a function”. Is there another classification approach that could be applied to the NDWI output to convert it into land/water classes?

//////////////////////////////////////////////////////////////
// Asset List
//////////////////////////////////////////////////////////////


// Outline for AOI (subset); loaded as asset on user account;
// A general polygon will work for testing purposes
var AOI_subset = AOI_subset;

// Image collection for Landsat 8 TOA
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA');
// Filter by AOI
var spatialFiltered = l8.filterBounds(AOI_subset);
// Filter by date
var temporalFiltered = spatialFiltered.filterDate('2017-01-01', '2018-12-31');
// Filter by cloud cover
var sorted = temporalFiltered.sort('CLOUD_COVER');
var scene = ee.Image(sorted.first());

// Get the median over time, in each band, in each pixel.
var median = temporalFiltered.median();
// Make a variable of visualization parameters.
var visParams = {bands: ['B4', 'B3', 'B2'], max: 0.1};

// Create an NDWI image, define visualization parameters and display.
var ndwi = median.normalizedDifference(['B3', 'B5']);
var ndwiViz = {min: 0.01, max: 0.75, palette: ['2ca25f', '0000FF']};

// Reclassify each image in NDWI collection and calculate composite
var NWDIreclassified = ndwi.map(function(img){
  return img.gte(0.4);
});

//////////////////////////////////////////////////////////////
// Otsu Segmentation
//////////////////////////////////////////////////////////////

// Use Otsu (1979) method for image segmentation to create water-land mask
// Based on code from Nicholas Clinton (Google): https://code.earthengine.google.com/6f3abb73b6642198485e36000024827a
// Not all code used from above link because was giving an error

// Compute the histogram of the NIR band.
var histogram = scene.select('B5').reduceRegion({
  reducer: ee.Reducer.histogram(255, 2),
  geometry: AOI_subset, 
  scale: 30,
  bestEffort: true
});
  • It is very well doable to map the otsu thresholding over an image collection. Please provide the part of the code until you have the image collection with NDWI to get help. – Kuik Feb 14 at 10:56
  • @Kuik, apologies, I forgot to paste in my code after writing the question. Edited to include now. – Dan Feb 14 at 17:17
1

So, you are asking two question at once. You would like to use the Otsu threshold to calculate a water mask on a per-image basis. Besides, you would like a watermask in a similar way for the median image. In your code, you have mixed them up. For that reason, your code is not working.

You tried to filter by cloud cover, however, as you do it now you will remain the full collection and just rearrange the order of the images. To filter by cloudcover, use:

var sorted = temporalFiltered.filter(ee.Filter.lt('CLOUD_COVER', 10));

First, make a function to calculate NDWI for a random image:

// Create a function to calculate NDWI per image
function calcNDWI(image){
  return image.normalizedDifference(['B3', 'B5']).rename('NDWI');
}

Then, define the Otsu thresholding function of Nicolas Clinton inside your code (which I have done, see the link). To get the Otsu threshold and water mask for the full image collection, map over you imageCollection, which you named 'sorted':

var mappedCollection = sorted.map(function(image){
  var NDWI = calcNDWI(image);
  var histogram = ee.Dictionary(NDWI.reduceRegion({
    reducer: ee.Reducer.histogram(100),
    geometry: AOI_subset, 
    scale: 30,
    bestEffort: true
  }).get('NDWI'));
  var threshold = otsu(histogram);
  var waterMask = NDWI.lt(ee.Number(threshold)).rename('waterMask');
  return image.addBands([NDWI, waterMask]).set('OtsuThresh', threshold);
});

For the median, you can do that in a similar way. But, the median is a SINGLE image, so you cannot map over it (which gave the error):

var NDWImedian = calcNDWI(median);
var histogram = ee.Dictionary(NDWImedian.reduceRegion({
  reducer: ee.Reducer.histogram(100),
  geometry: AOI_subset, 
  scale: 30,
  bestEffort: true
}).get('NDWI'));
var threshold = otsu(histogram);
var waterMask = NDWImedian.lt(ee.Number(threshold)).rename('waterMask');
var median = median.addBands([NDWImedian, waterMask]).set('OtsuThresh', threshold);

Then visualize some stuff in the console and map. Note that inside a mapped function, occuring on the server-side, you cannot print something (e.g. the histogram), so I did that only for the median image:

// Print
print('size full collection', sorted.size());
print('The full collection', mappedCollection);
print('Histogram of the median', ui.Chart.image.histogram(median.select('NDWI'), AOI_subset, 30));
print('Otsu thresholds', mappedCollection.aggregate_array('OtsuThresh'));

// Make a variable of visualization parameters.
var visParams = {bands: ['B4', 'B3', 'B2'], min: 0, max: 0.35};
var ndwiViz = {min: 0.01, max: 0.75, palette: ['2ca25f', '0000FF'], bands: 'NDWI'};

Map.addLayer(median.updateMask(median.select('waterMask')), visParams, 'True color median water masked');
Map.addLayer(mappedCollection.first().updateMask(mappedCollection.first().select('waterMask')), visParams, 'True color first image water masked');

Link script

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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