Wednesday, 24 June 2015

google earth engine - Cloud mask in Surface Reflectance Landsat 8 test


I've noticed that my cloud mask wasn't working, so I've tried this simple test: https://code.earthengine.google.com/50699c2eaa1a873ccd28f26c583c5a45


But my data uses the Surface Reflectance Landsat 8 imagery so, I just changed to that, and since there's no fmask band in this, I've changed to the 'pixel_qa' that also make those distinctions. I thought it was the same, but isn't working.


Code:


//Choose country using GEE Feature Collection

var region = ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw').filterMetadata('Country', 'equals', 'Portugal');


//Add region outline to layer ‐ for selected countries

Map.addLayer(region,{}, 'Portugal');

var landsat8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')

// Fmask classification values
var FMASK_CLEAR_GROUND = 0;
var FMASK_WATER = 2;
var FMASK_CLOUD_SHADOW = 3;

var FMASK_SNOW = 4;
var FMASK_CLOUD = 5;

var mosaic = landsat8
.filterBounds(region)
.filterDate('2017-08-01', '2017-10-11')
.mosaic();

// Update the mask on our mosaic to mask cloud and cloud shadow pixels
var fmask = mosaic.select('pixel_qa');

var cloudMask = fmask.neq(FMASK_CLOUD).and(fmask.neq(FMASK_CLOUD_SHADOW));
var maskedMosaic = mosaic.updateMask(cloudMask);

Map.addLayer(fmask, {min:0, max:5, palette:'green, blue, black, cyan, pink, white'}, 'Fmask');
Map.addLayer(maskedMosaic.select('B4'), {min:0, max:0.5, palette:'yellow, green'}, 'Masked NIR');

Answer



Surface reflectance doesn't have fmask, neither cfmask (used for cloud masking in old Landsat scenes).


You need to use Quality Assessment layer, this topic I answered before with MODIS in R. This case apply at the same way, you need to decode bit flags from pixel_qa band.


Check the documentation, you need to use Bit 3 and 5 to mask clouds shadows and cloud:


enter image description here



Function used is:


var getQABits = function(image, start, end, newName) {
// Compute the bits we need to extract.
var pattern = 0;
for (var i = start; i <= end; i++) {
pattern += Math.pow(2, i);
}
// Return a single band image of the extracted QA bits, giving the band
// a new name.
return image.select([0], [newName])

.bitwiseAnd(pattern)
.rightShift(start);
};

Applied to your code:


//Choose country using GEE Feature Collection

var region = ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw').filterMetadata('Country', 'equals', 'Portugal');

//Add region outline to layer ‐ for selected countries


Map.addLayer(region,{}, 'Portugal');

var landsat8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')

// Fmask classification values var FMASK_CLEAR_GROUND = 0; var FMASK_WATER = 2; var FMASK_CLOUD_SHADOW = 3; var FMASK_SNOW = 4; var FMASK_CLOUD = 5;

var mosaic = landsat8 .filterBounds(region) .filterDate('2017-08-01', '2017-10-11') .mosaic();

var getQABits = function(image, start, end, newName) {

// Compute the bits we need to extract.
var pattern = 0;
for (var i = start; i <= end; i++) {
pattern += Math.pow(2, i);
}
// Return a single band image of the extracted QA bits, giving the band
// a new name.
return image.select([0], [newName])
.bitwiseAnd(pattern)
.rightShift(start);

};

// A function to mask out cloudy pixels.
var cloud_shadows = function(image) {
// Select the QA band.
var QA = image.select(['pixel_qa']);
// Get the internal_cloud_algorithm_flag bit.
return getQABits(QA, 3,3, 'Cloud_shadows').eq(0);
// Return an image masking out cloudy areas.
};


// A function to mask out cloudy pixels.
var clouds = function(image) {
// Select the QA band.
var QA = image.select(['pixel_qa']);
// Get the internal_cloud_algorithm_flag bit.
return getQABits(QA, 5,5, 'Cloud').eq(0);
// Return an image masking out cloudy areas.
};


var maskClouds = function(image) {
var cs = cloud_shadows(image);
var c = clouds(image);
image = image.updateMask(cs);
return image.updateMask(c);
};

var mosaic_free = maskClouds(mosaic);

var visParams = {bands: ['B4', 'B3', 'B2'],min: [0,0,0],max: [2000, 2000, 2000]};


Map.addLayer(mosaic, visParams, 'With clouds');
Map.addLayer(mosaic_free, visParams, 'Cloud free');

Link: https://code.earthengine.google.com/d653edd684a02416d3910182cc465684


With clouds:


enter image description here


Without clouds:


enter image description here


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...