Skip to content

Latest commit

 

History

History
380 lines (306 loc) · 16.5 KB

03.cloudmaskTOA.md

File metadata and controls

380 lines (306 loc) · 16.5 KB
Raw

geeguide

03. Cloud Masking of Landsat and Sentinel2 TOA Products

Objective

  • Landsat TOA Cloud and Shadow Mask
  • Sentinel2 TOA Cloud and Shadow Mask

General Information

  • Cloud contamination will lower NDVI values, measures like the timing of ‘green up’ or peak maturity would appear later than they actually occurred (USGS, 2019). Therefore, cloud should be masked effectively and 'aggressively' when the images are used to analysis crop’s phenology.
  • For Landsat TOA product, BQA band can be used to mask out cloud with decent performance. BQA band was generated using CFMask Algorithm. Or ee.Algorithms.Landsat.simpleCloudScore() is a built-in GEE function can provide a shortcut to a rudimentary cloud scoring algorithm GEE Documentation. Behide the scene of the simpleCloudScore algorithm is described in this gee script.
  • However, without a thermal band, current cloud detection methods for Sentinel 2 (eg.Fmask, Sen2Cor, MAJA, etc) could not provide satisfy results. Bright land surface such as white roof, sand beach, bare land, high turbidity surface water, etc. could be miss identified as cloud. Sentinel-2 optimal cloud mask is the main issue in hls.gsfc.nasa.gov initiative (Claverie et.al, 2018)
  • I have been working on a S2 cloud mask based on Supervised Classification and QA60 was used as training data. The results look promissing but it may be updated later.

Core script

Cloud Mask Landsat8

// Landsat 8 Cloud Masking Example

var RADIX = 2;  // Radix for binary (base 2) data.

// Reference a sample Landsat 8 TOA image.
var image = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_043034_20180202');
// Extract the QA band.
var image_qa = image.select('BQA');

var extractQABits = function (qaBand, bitStart, bitEnd) {
  var numBits = bitEnd - bitStart + 1;
  var qaBits = qaBand.rightShift(bitStart).mod(Math.pow(RADIX, numBits));
  //Map.addLayer(qaBits, {min:0, max:(Math.pow(RADIX, numBits)-1)}, 'qaBits');
  return qaBits;
};

// Create a mask for the dual QA bit "Cloud Confidence".
var bitStartCloudConfidence = 5;
var bitEndCloudConfidence = 6;
var qaBitsCloudConfidence = extractQABits(image_qa, bitStartCloudConfidence, bitEndCloudConfidence);
// Test for clouds, based on the Cloud Confidence value.
var testCloudConfidence = qaBitsCloudConfidence.gte(2);

// Create a mask for the dual QA bit "Cloud Shadow Confidence".
var bitStartShadowConfidence = 7;
var bitEndShadowConfidence = 8;
var qaBitsShadowConfidence = extractQABits(image_qa, bitStartShadowConfidence, bitEndShadowConfidence);
// Test for shadows, based on the Cloud Shadow Confidence value.
var testShadowConfidence = qaBitsShadowConfidence.gte(2);

// Calculate a composite mask and apply it to the image.   
var maskComposite = (testCloudConfidence.or(testShadowConfidence)).not();
var imageMasked = image.updateMask(maskComposite);

Map.addLayer(image, {bands:"B4,B3,B2", min:0, max:0.3}, 'original image', false);
Map.addLayer(image.select('BQA'), {min:0, max:10000}, 'BQA', false);
Map.addLayer(testCloudConfidence.mask(testCloudConfidence), {min:0, max:1, palette:'grey,white'}, 'testCloudConfidence');
Map.addLayer(testShadowConfidence.mask(testShadowConfidence), {min:0, max:1, palette:'grey,black'}, 'testShadowConfidence');
Map.addLayer(maskComposite.mask(maskComposite), {min:0, max:1, palette:'green'}, 'maskComposite', false);
Map.addLayer(imageMasked, {bands:"B4,B3,B2", min:0, max:0.3}, 'imageMasked');

//This function even calculates percentage of cloud  coverage over a particular region. Useful for screening too cloudy images

function scoreL8_TOA (image) {
  var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud').rename('cloudScore');
  var cloudiness = cloud.reduceRegion({
    reducer: 'mean', 
    geometry: region, 
    scale: 30,
  });
  return image.set(cloudiness).addBands(cloud);
}

function maskL8_TOA (img){
  var mask = img.select('cloudScore').lt(20);
  return img.updateMask(mask)
}

Cloud Mask Sentinel 2

//CloudScore originally written by Matt Hancher and adapted for S2 data by Ian Housman
//////////////////////////////////////////////////////////
///https://code.earthengine.google.com/c0b316ba2b56121b82318ffd1e5c42de
//User Params
var startYear = 2016;
var endYear = 2017;
var startJulian = 190;
var endJulian = 250;
var cloudThresh =20;//Ranges from 1-100.Lower value will mask more pixels out. Generally 10-30 works well with 20 being used most commonly 
var cloudHeights = ee.List.sequence(200,10000,250);//Height of clouds to use to project cloud shadows
var irSumThresh =0.35;//Sum of IR bands to include as shadows within TDOM and the shadow shift method (lower number masks out less)
var dilatePixels = 2; //Pixels to dilate around clouds
var contractPixels = 1;//Pixels to reduce cloud mask and dark shadows by to reduce inclusion of single-pixel comission errors

//////////////////////////////////////////////////////////
var vizParams = {'min': 0.05,'max': [0.3,0.6,0.35],   'bands':'swir1,nir,red'};
// var vizParams = {bands: ['red', 'green', 'blue'], min: 0, max: 0.3};
//////////////////////////////////////////////////////////////////////////
var rescale = function(img, exp, thresholds) {
    return img.expression(exp, {img: img})
        .subtract(thresholds[0]).divide(thresholds[1] - thresholds[0]);
  };
  
////////////////////////////////////////
////////////////////////////////////////
// Cloud masking algorithm for Sentinel2
//Built on ideas from Landsat cloudScore algorithm
//Currently in beta and may need tweaking for individual study areas
function sentinelCloudScore(img) {
  

  // Compute several indicators of cloudyness and take the minimum of them.
  var score = ee.Image(1);
  
  // Clouds are reasonably bright in the blue and cirrus bands.
  score = score.min(rescale(img, 'img.blue', [0.1, 0.5]));
  score = score.min(rescale(img, 'img.cb', [0.1, 0.3]));
  score = score.min(rescale(img, 'img.cb + img.cirrus', [0.15, 0.2]));
  
  // Clouds are reasonably bright in all visible bands.
  score = score.min(rescale(img, 'img.red + img.green + img.blue', [0.2, 0.8]));

  
  //Clouds are moist
  var ndmi = img.normalizedDifference(['nir','swir1']);
  score=score.min(rescale(ndmi, 'img', [-0.1, 0.1]));
  
  // However, clouds are not snow.
  var ndsi = img.normalizedDifference(['green', 'swir1']);
  score=score.min(rescale(ndsi, 'img', [0.8, 0.6]));
  
  score = score.multiply(100).byte();
 
  return img.addBands(score.rename('cloudScore'));
}
//////////////////////////////////////////////////////////////////////////
// Function to mask clouds using the Sentinel-2 QA band.
function maskS2clouds(image) {
  var qa = image.select('QA60').int16();
  
  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = Math.pow(2, 10);
  var cirrusBitMask = Math.pow(2, 11);
  
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
             qa.bitwiseAnd(cirrusBitMask).eq(0));

  // Return the masked and scaled data.
  return image.updateMask(mask);
}
//////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//Function for finding dark outliers in time series
//Masks pixels that are dark, and dark outliers
function simpleTDOM2(c){
  var shadowSumBands = ['nir','swir1'];
  var irSumThresh = 0.4;
  var zShadowThresh = -1.2;
  //Get some pixel-wise stats for the time series
  var irStdDev = c.select(shadowSumBands).reduce(ee.Reducer.stdDev());
  var irMean = c.select(shadowSumBands).mean();
  var bandNames = ee.Image(c.first()).bandNames();
  print('bandNames',bandNames);
  //Mask out dark dark outliers
  c = c.map(function(img){
    var z = img.select(shadowSumBands).subtract(irMean).divide(irStdDev);
    var irSum = img.select(shadowSumBands).reduce(ee.Reducer.sum());
    var m = z.lt(zShadowThresh).reduce(ee.Reducer.sum()).eq(2).and(irSum.lt(irSumThresh)).not();
    
    return img.updateMask(img.mask().and(m));
  });
  
  return c.select(bandNames);
}
////////////////////////////////////////////////////////
/////////////////////////////////////////////
/***
 * Implementation of Basic cloud shadow shift
 * 
 * Author: Gennadii Donchyts
 * License: Apache 2.0
 */
function projectShadows(cloudMask,image,cloudHeights){
  var meanAzimuth = image.get('MEAN_SOLAR_AZIMUTH_ANGLE');
  var meanZenith = image.get('MEAN_SOLAR_ZENITH_ANGLE');
  ///////////////////////////////////////////////////////
  // print('a',meanAzimuth);
  // print('z',meanZenith)
  
  //Find dark pixels
  var darkPixels = image.select(['nir','swir1','swir2']).reduce(ee.Reducer.sum()).lt(irSumThresh)
    .focal_min(contractPixels).focal_max(dilatePixels)
  ;//.gte(1);
  
  
  //Get scale of image
  var nominalScale = cloudMask.projection().nominalScale();
  //Find where cloud shadows should be based on solar geometry
  //Convert to radians
  var azR =ee.Number(meanAzimuth).add(180).multiply(Math.PI).divide(180.0);
  var zenR  =ee.Number(meanZenith).multiply(Math.PI).divide(180.0);
  
  
 
  //Find the shadows
  var shadows = cloudHeights.map(function(cloudHeight){
    cloudHeight = ee.Number(cloudHeight);
    
    var shadowCastedDistance = zenR.tan().multiply(cloudHeight);//Distance shadow is cast
    var x = azR.sin().multiply(shadowCastedDistance).divide(nominalScale);//X distance of shadow
    var y = azR.cos().multiply(shadowCastedDistance).divide(nominalScale);//Y distance of shadow
    // print(x,y)
   
    return cloudMask.changeProj(cloudMask.projection(), cloudMask.projection().translate(x, y));
    
    
  });
  
  
  var shadowMask = ee.ImageCollection.fromImages(shadows).max();
  // Map.addLayer(cloudMask.updateMask(cloudMask),{'min':1,'max':1,'palette':'88F'},'Cloud mask');
  // Map.addLayer(shadowMask.updateMask(shadowMask),{'min':1,'max':1,'palette':'880'},'Shadow mask');
  
  //Create shadow mask
  shadowMask = shadowMask.and(cloudMask.not());
  shadowMask = shadowMask.and(darkPixels).focal_min(contractPixels).focal_max(dilatePixels);
  
  var cloudShadowMask = shadowMask.or(cloudMask);
  
  image = image.updateMask(cloudShadowMask.not()).addBands(shadowMask.rename(['cloudShadowMask']));
  return image;
}
//////////////////////////////////////////////////////
//Function to bust clouds from S2 image
function bustClouds(img){
  img = sentinelCloudScore(img);
  img = img.updateMask(img.select(['cloudScore']).gt(cloudThresh).focal_min(contractPixels).focal_max(dilatePixels).not());
  return img;
}
//////////////////////////////////////////////////////
//Function for wrapping the entire process to be applied across collection
function wrapIt(img){
  img = sentinelCloudScore(img);
  var cloudMask = img.select(['cloudScore']).gt(cloudThresh)
    .focal_min(contractPixels).focal_max(dilatePixels)

  img = projectShadows(cloudMask,img,cloudHeights);

  return img;
}
//////////////////////////////////////////////////////
//Function to find unique values of a field in a collection
function uniqueValues(collection,field){
    var values  =ee.Dictionary(collection.reduceColumns(ee.Reducer.frequencyHistogram(),[field]).get('histogram')).keys();
    
    return values;
  }
//////////////////////////////////////////////////////
//Function to simplify data into daily mosaics
function dailyMosaics(imgs){
  //Simplify date to exclude time of day
  imgs = imgs.map(function(img){
  var d = ee.Date(img.get('system:time_start'));
  var day = d.get('day');
  var m = d.get('month');
  var y = d.get('year');
  var simpleDate = ee.Date.fromYMD(y,m,day);
  return img.set('simpleTime',simpleDate.millis());
  });
  
  //Find the unique days
  var days = uniqueValues(imgs,'simpleTime');
  
  imgs = days.map(function(d){
    d = ee.Number.parse(d);
    d = ee.Date(d);
    var t = imgs.filterDate(d,d.advance(1,'day'));
    var f = ee.Image(t.first());
    t = t.mosaic();
    t = t.set('system:time_start',d.millis());
    t = t.copyProperties(f);
    return t;
    });
    imgs = ee.ImageCollection.fromImages(imgs);
    
    return imgs;
}
//////////////////////////////////////////////////////
//Get some s2 data
var s2s = ee.ImageCollection('COPERNICUS/S2')
                  .filter(ee.Filter.calendarRange(startYear,endYear,'year'))
                  .filter(ee.Filter.calendarRange(startJulian,endJulian))
                  .filterBounds(geometry)
                  .map(function(img){
                    
                    var t = img.select([ 'B1','B2','B3','B4','B5','B6','B7','B8','B8A', 'B9','B10', 'B11','B12']).divide(10000);//Rescale to 0-1
                    t = t.addBands(img.select(['QA60']));
                    var out = t.copyProperties(img).copyProperties(img,['system:time_start']);
                  return out;
                    })
                    .select(['QA60', 'B1','B2','B3','B4','B5','B6','B7','B8','B8A', 'B9','B10', 'B11','B12'],['QA60','cb', 'blue', 'green', 'red', 're1','re2','re3','nir', 'nir2', 'waterVapor', 'cirrus','swir1', 'swir2']);

//Convert to daily mosaics to avoid redundent observations in MGRS overlap areas and edge artifacts for shadow masking
s2s = dailyMosaics(s2s);
print(s2s);
//////////////////////////////////////////////////
//Optional- View S2 daily mosaics
// days.getInfo().map(function(d){
//   var ds = new Date(parseInt(d))
//   d = ee.Date(ee.Number.parse(d))
  
//   print(ds)
//   var s2sT = ee.Image(s2s.filterDate(d,d.advance(1,'minute')).first());
 
//   Map.addLayer(s2sT,vizParams,ds,false);
// })
/////////////////////////////////////////////////////
//Look at individual image
// var s2 = ee.Image(s2s.first());
// Map.addLayer(s2,vizParams,'Before Masking',false);

// s2 = sentinelCloudScore(s2);
// var cloudScore = s2.select('cloudScore');
// Map.addLayer(cloudScore,{'min':0,'max':100},'CloudScore',false);


// var s2CloudMasked  = bustClouds(s2);
// Map.addLayer(s2CloudMasked,vizParams,'After CloudScore Masking',false);

// var s2CloudShadowMasked = wrapIt(s2);
// Map.addLayer(s2CloudShadowMasked,vizParams,'After CloudScore and Shadow Masking',false);

// var s2MaskedQA = maskS2clouds(s2) ;
// print(s2)
// Map.addLayer(s2MaskedQA,vizParams,'After QA Masking',false);


// var s2TDOMMasked = ee.Image(simpleTDOM2(s2s.map(bustClouds)).first());
// Map.addLayer(s2TDOMMasked,vizParams,'After CloudScore and TDOM Masking',false);


/////////////////////////////////////////////
//Look at mosaics
//Get the raw mosaic and median
Map.addLayer(s2s.mosaic(),vizParams,'Raw Mosaic', false);
Map.addLayer(s2s.median(),vizParams,'Raw Median', false);


//Bust clouds using BQA method
var s2MaskedQA = s2s.map(maskS2clouds);
Map.addLayer(s2MaskedQA.mosaic(),vizParams,'QA Cloud Masked Mosaic',false);
Map.addLayer(s2MaskedQA.median(),vizParams,'QA Cloud Masked Median',false);



//Bust clouds using cloudScore method
var s2sMosaic = s2s.map(bustClouds);
Map.addLayer(s2sMosaic.mosaic(),vizParams,'CloudScore Masked Mosaic', false);
Map.addLayer(s2sMosaic.median(),vizParams,'CloudScore Masked Median', false);


//Bust clouds using cloudScore and shadows using TDOM
var s2TDOM = simpleTDOM2(s2sMosaic);
Map.addLayer(s2TDOM.mosaic(),vizParams,'CloudScore and TDOM Masked Mosaic', false);
Map.addLayer(s2TDOM.median(),vizParams,'CloudScore and TDOM Masked Median', false);

//Bust clouds using cloudScore and shadows using shadow shift method
var s2sMosaicCloudsProjectedDark = s2s.map(wrapIt)
Map.addLayer(s2sMosaicCloudsProjectedDark.mosaic(), vizParams, 'Cloud Masked + Projected Shadows + Dark + Mosaic ',false);
Map.addLayer(s2sMosaicCloudsProjectedDark.median(), vizParams, 'Cloud Masked + Projected Shadows + Dark + Median ',false);

Visualization and Checking

References

  1. Cloud masking Sentinel 2, Google Earth Engine Developers.
  2. Landsat 8 BQA - using cloud confidence to create a cloud mask, Stack Exchange