Find Water and Remove Clouds with Fmask on Google Earth Engine

Visualizing Fmask pixel classifications and the near infrared band with cloud affected pixels removed.

Fmask is a fantastic algorithm for identifying clouds, water, snow and clear ground pixels in your satellite data. It’s possible to run Fmask yourself, but if you’re working with Landsat data on Google Earth Engine you’re in luck — Google provides calculated Fmask output for the Landsat 4, 5, 7 and 8 satellites. Here I’ll give a quick run through on how to use Fmask layers on Earth Engine to mask clouds (and cloud shadows) and calculate surface water area.

First, grab a mosaic of data that includes the fmask classification band to work with. Here I’m using Landsat 8.

var polygon = ee.Geometry.Polygon([[
 [-119.33364769821117, 46.05123532178373],
 [-119.3233620672313, 45.869732769408905],
 [-119.04111088542663, 45.873079023065166],
 [-119.0396574679861, 46.045448840018565]
]]);
var mosaic = ee.ImageCollection(‘LANDSAT/LC8_L1T_TOA_FMASK’)
 .filterBounds(polygon)
 .filterDate(‘2016–08–01’, ‘2016–08–30’)
 .mosaic();

Next make a binary image, cloudMask, that contains zeros (denoting invalid pixels) wherever the fmask band has identified pixels as clouds or cloud shadows.

// Fmask classification values
var FMASK_CLOUD = 4;
var FMASK_CLOUD_SHADOW = 2;
var fmask = mosaic.select('fmask');
var cloudMask = fmask.neq(FMASK_CLOUD)
  .and(fmask.neq(FMASK_CLOUD_SHADOW));

Finally, mask out the cloudy pixels in the mosaic. The updateMask function will add our mask of cloud pixels to the mosaic’s existing mask, leaving us with a clean image to work with.

var maskedMosaic = mosaic.updateMask(cloudMask);

Here’s a link to the complete code. If you don’t have Earth Engine, that link will let you know how to sign up and you can view this gist with the same code in the meantime.

Visualizing

Now that we have cloud free data, lets have some fun with it. For this section you may want to tinker with the code yourself.

I like to visualize the Fmask band to get a sense of how it’s classifying. The ordering of the color palette here matches the numerical values Fmask uses: 0 - Ground (green), 1 - Water (blue), 2 - Cloud Shadow (black), 3 - Snow (cyan), 4 - Clouds (white).

Map.setCenter(-119.34, 45.97, 8);
Map.addLayer(fmask, {min:0, max:4, palette:'green, blue, black, cyan, white'}, 'Fmask');

You can also visualize your cleaned up data, lets look at only the cloud free near infrared pixels from band 5.

Map.addLayer(maskedMosaic.select('B5'), {min:0, max:0.7, palette:'yellow, green'}, 'Masked NIR');

Calculating Statistics

We can also use Fmask to grab some basic stats about land/water cover within a region. The frequencyHistogram will retrieve counts of all pixel classes within our area of interest. We can use that to calculate percent cover in a region or the total area of pixels classified as water.

var PIXEL_SCALE = 30; // Meters. Resolution of most Landsat 8 bands
var PIXEL_AREA = PIXEL_SCALE * PIXEL_SCALE;
var FMASK_WATER = 1;
var regionCoverHistogram = fmask.reduceRegion(
  ee.Reducer.frequencyHistogram(), polygon, PIXEL_SCALE);
print('Fmask class pixel counts', regionCoverHistogram);
// The return value of .get() is untyped, 
// cast it with ee.Number and ee.Dictionary.
var waterPixelCount =
  ee.Dictionary(regionCoverHistogram.get('fmask'))
    .get(FMASK_WATER.toString());
var waterArea = ee.Number(waterPixelCount).multiply(PIXEL_AREA);
print('Water Area (sq meters) in region', waterArea);