#3: create a cloud free composite

Remote sensing Leave a comment

A modular land cover system part 3: cloud masking

Step 1: create a new script and call it cloud_module

Step 2: copy the code below and store it in the script


var cloudScoreThresh = 30;
var cloudScorePctl = 0;
var contractPixels = 1.5;
var dilatePixels = 2.5;

////////////////////////////////////////////////////////////////////////////////
// Compute a cloud score and adds a band that represents the cloud mask.
// Cloud masking algorithm for Sentinel2. Built on ideas from Landsat cloudScore
// algorithm. Currently in beta and may need tweaking for individual study areas
exports.sentinelCloudScore = function(s2s) {

  function getCloudScore(img){
    // Compute several indicators of cloudyness and take the minimum of them.
    var score = ee.Image(1);
    var blueCirrusScore = ee.Image(0);

    // Clouds are reasonably bright in the blue or cirrus bands.
    //Use .max as a pseudo OR conditional
    blueCirrusScore = blueCirrusScore.max(rescale(img, 'img.blue', [0.1, 0.5]));
    blueCirrusScore = blueCirrusScore.max(rescale(img, 'img.cb', [0.1, 0.5]));
    blueCirrusScore = blueCirrusScore.max(rescale(img, 'img.cirrus', [0.1, 0.3]));
    score = score.min(blueCirrusScore);

    // 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 reasonably bright in all infrared bands.
    score = score.min(rescale(img, 'img.nir + img.swir1 + img.swir2', [0.3, 0.8]));

    // 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();
    score = score.clamp(0,100);
    return img.addBands(score.rename(['cloudScore']));
    }

  function maskScore(img){
    //var cloudMask = img.select(['cloudScore']).subtract(minCloudScore).lt(cloudScoreThresh)
    // .focal_max(contractPixels).focal_min(dilatePixels).rename('cloudMask');
    var cloudMask = img.select(['cloudScore']).lt(cloudScoreThresh).focal_max(contractPixels).focal_min(dilatePixels).rename('cloudMask');
    return img.updateMask(cloudMask).addBands(cloudMask);
    }

  s2s = s2s.map(getCloudScore);

  // Find low cloud score pctl for each pixel to avoid comission errors
  var minCloudScore = s2s.select(['cloudScore']).reduce(ee.Reducer.percentile([cloudScorePctl]));

  s2s = s2s.map(maskScore);

  return s2s;
};

////////////////////////////////////////////////////////////////////////////////
// Function to mask clouds using the Sentinel-2 QA band.
exports.QAMaskCloud= function(collection) {

function maskClouds(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);
}

collection = collection.map(maskClouds);

return collection;

};

////////////////////////////////////////////////////////////////////////////////
// Helper function to apply an expression and linearly rescale the output.
// Used in the sentinelCloudScore function below.
function rescale(img, exp, thresholds) {
return img.expression(exp, {img: img})
.subtract(thresholds[0]).divide(thresholds[1] - thresholds[0]);
}

Step 4: Open the main script and add the lines below

print("applying QA cloud mask");
s2 = clouds.QAMaskCloud(s2);
Map.addLayer(ee.Image(s2.first()),{min:0,max:??,bands:"??,nir,red"},"after QA cloudmask");
print(ee.Image(s2.first()));

print("applying sentinel cloud mask");
s2 = clouds.sentinelCloudScore(s2);
Map.addLayer(ee.Image(s2.??()),{min:0,max:0.6,bands:"swir2,nir,??"},"after s2 cloudmask");
print(ee.Image(s2.first()));

Step 5: Run the script

link

Leave a Reply