GEE C15 Clouds and Image Compositing 去云和影像合成

学习要点：

1.理解和应用特定卫星的云掩膜功能

2.整合不同传感器的影像

3.使用焦点函数填补数据空白

一、Cloud Filter and Cloud Mask 云过滤和云掩膜

1.1 大尺度国际边界（LSIB）数据集

 

1.2 影像缩放因子

• image.select('SR_B.')：选择影像的表面反射率波段(SR_B)。点号(.)表示模式匹配，意味着它选择所有以SR_B开头的波段。
• .multiply(0.0000275).add(-0.2)：对选定的光学波段应用缩放因子和偏移量。这个操作是将数字量化值转换为反射率值所必需的。
• image.select('ST_B.*')：使用正则表达式选择表面温度波段(ST_B)，匹配任何以ST_B开头的波段名称。
• .multiply(0.00341802).add(149.0)：对热波段应用缩放因子和偏移量，将数字量化值转换为开尔文温度值。
• 函数返回带有缩放后的光学和热波段的原始影像（addBands中的true参数意味着替换原有波段）。

 

 1.3 Landsat5-8 QA波段

质量评估波段（the Quality Assessment (QA) Band）。根据USGS Landsat Missions webpage，“QA通过标示哪个像素可能受仪器或云层影响，从而提高了科学研究的完整性，这已经被有效地利用起来。”简而言之QA波段让最终用户更容易地识别 “坏”像素、挑选出“好”数据、提高更准确更精确的结果。这个新波段可以有很多用途，如辨别积雪地区的城市区域。

// Define the cloud mask function. 
function maskSrClouds(image) { 
    // Bit 0 - Fill 
    // Bit 1 - Dilated Cloud 
    // Bit 2 - Cirrus 
    // Bit 3 - Cloud 
    // Bit 4 - Cloud Shadow 
    var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0); 

    var saturationMask = image.select('QA_RADSAT').eq(0); 
    return image.updateMask(qaMask).updateMask(saturationMask); 
} 

代码解释：

1. 定义云和阴影掩膜：函数通过分析QA_PIXEL波段（质量评估像素波段）来创建一个掩膜，以识别和移除云、云影、雪、以及其他可能干扰图像分析的因素。

   • image.select('QA_PIXEL')：选择图像中的QA_PIXEL波段。这个波段包含了关于像素质量的编码信息，例如是否被云、云影或其他因素覆盖。
   • .bitwiseAnd(parseInt('11111', 2))：通过位运算与操作（bitwiseAnd）和一个二进制掩膜（11111，转换为十进制后与QA_PIXEL波段进行比较）来识别特定的质量标志。这个二进制数代表了云、云影等的不同质量标志位。
     • Bit 0 - Fill：填充位，通常用于指示无数据区域。
     • Bit 1 - Dilated Cloud：扩展云，指示周围可能被云覆盖的区域。
     • Bit 2 - Cirrus：卷云，高层大气中的薄云。
     • Bit 3 - Cloud：云。
     • Bit 4 - Cloud Shadow：云影。
   • .eq(0)：此操作用于选择那些所有相关位都为0的像素，即不包含上述任何质量问题的像素。

2. 定义饱和掩膜：另外，通过检查QA_RADSAT波段（辐射饱和度和地形遮蔽的质量评估）来创建一个饱和掩膜。

   • image.select('QA_RADSAT').eq(0)：选择QA_RADSAT波段，并找出其中值为0的像素，即没有辐射饱和或地形遮蔽问题的像素。

3. 应用掩膜：最后，函数将这些掩膜应用于原始图像，以排除被云、云影、饱和等问题覆盖的像素。

• .updateMask(qaMask)：首先应用由QA_PIXEL波段生成的掩膜，移除那些包含云、云影等问题的像素。
• .updateMask(saturationMask)：然后应用由QA_RADSAT波段生成的掩膜，移除那些存在辐射饱和或地形遮蔽问题的像素。

1.bitwiseAnd 函数

位运算，二进制表示中的位为1 (1 × 1 = 1);否则，结果为0 (1 × 0 = 0和0 × 0 = 0)。

2.parseInt 函数

parseInt函数解析字符串参数(在本例中为5个字符的字符串' 11111 ')，并返回指定的以2为基数的整数。

3.完整代码

// ---------- Section 1 -----------------


// Define the AOI. 
// var country = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
//     .filter(ee.Filter.equals('country_na', 'China')); 
var country = table;
print(table);
// Center the Map. The second parameter is zoom level. 
Map.centerObject(country, 13);
Map.addLayer(country,{},'GZL baohuqu',1,0.3);

// Define time variables. 
var startDate = '2019-01-01'; 
var endDate = '2019-12-31'; 
// Load and filter the Landsat 8 collection. 
var landsat8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') 
    .filterBounds(country) 
    .filterDate(startDate, endDate); 
print('landsat8原始影像:', landsat8);
// Apply scaling factors. 
function applyScaleFactors(image) { 
    var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2); 
    var thermalBands = image.select('ST_B.*').multiply(0.00341802) .add(149.0); 
    return image.addBands(opticalBands, null, true) 
                .addBands(thermalBands, null, true); 
} 
landsat8 = landsat8.map(applyScaleFactors); 
print('landsat8对反射率SR_B和温度ST_B缩放因子后的影像：', landsat8);
print('Size landsat8 collection',landsat8.size());


// Create composite. 
var composite = landsat8.median().clip(country); 
//the function of landsat8.median() and landsat8.reduce(ee.Reducer.median()) are same
//var composite = landsat8.reduce(ee.Reducer.median()).clip(country);
var visParams = { 
    bands: ['SR_B4', 'SR_B3', 'SR_B2'], 
    min: 0, 
    max: 0.2 
}; 
Map.addLayer(composite, visParams, 'L8 Composite');


// Filter by the CLOUD_COVER property. 
var landsat8FiltClouds = landsat8 
    .filterBounds(country) 
    .filterDate(startDate, endDate) 
    .filter(ee.Filter.lessThan('CLOUD_COVER', 50)); 
print('Size landsat8FiltClouds collection', landsat8FiltClouds.size()); 
// Create a composite from the filtered imagery. 
var compositeFiltClouds = landsat8FiltClouds.median().clip(country); 
Map.addLayer(compositeFiltClouds, visParams, 'L8 Composite cloud filter'); 

//QAmask
// Define the cloud mask function. 
function maskSrClouds(image) { 
    // Bit 0 - Fill 
    // Bit 1 - Dilated Cloud 
    // Bit 2 - Cirrus 
    // Bit 3 - Cloud 
    // Bit 4 - Cloud Shadow 
    var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0); 

    var saturationMask = image.select('QA_RADSAT').eq(0); 
    return image.updateMask(qaMask).updateMask(saturationMask); 
} 
// Apply the cloud mask to the collection. 
var landsat8FiltMasked = landsat8FiltClouds.map(maskSrClouds); 
print('Size landsat8FiltClouds collection maskQA', landsat8FiltClouds.size()); 
// Create a composite. 
var landsat8compositeMasked = landsat8FiltMasked.median().clip(country); 
Map.addLayer(landsat8compositeMasked, visParams, 'L8 composite masked'); 

 

二、Incorporating Data from Other Satellites 整合其他传感器的影像

2.1 代码

// ---------- Section 1 -----------------


// Define the AOI. 
// var country = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
//     .filter(ee.Filter.equals('country_na', 'China')); 
var country = table;
print(table);
// Center the Map. The second parameter is zoom level. 
Map.centerObject(country, 13);
Map.addLayer(country,{},'GZL baohuqu',1,0.3);

// Define time variables. 
var startDate = '2019-01-01'; 
var endDate = '2019-12-31'; 
// Load and filter the Landsat 8 collection. 
var landsat8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') 
    .filterBounds(country) 
    .filterDate(startDate, endDate); 
print('landsat8原始影像:', landsat8);
// Apply scaling factors. 
function applyScaleFactors(image) { 
    var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2); 
    var thermalBands = image.select('ST_B.*').multiply(0.00341802) .add(149.0); 
    return image.addBands(opticalBands, null, true) 
                .addBands(thermalBands, null, true); 
} 
landsat8 = landsat8.map(applyScaleFactors); 
print('landsat8对反射率SR_B和温度ST_B缩放因子后的影像：', landsat8);
print('Size landsat8 collection',landsat8.size());


// Create composite. 
var composite = landsat8.median().clip(country); 
//the function of landsat8.median() and landsat8.reduce(ee.Reducer.median()) are same
//var composite = landsat8.reduce(ee.Reducer.median()).clip(country);
var visParams = { 
    bands: ['SR_B4', 'SR_B3', 'SR_B2'], 
    min: 0, 
    max: 0.2 
}; 
Map.addLayer(composite, visParams, 'L8 Composite');


// Filter by the CLOUD_COVER property. 
var landsat8FiltClouds = landsat8 
    .filterBounds(country) 
    .filterDate(startDate, endDate) 
    .filter(ee.Filter.lessThan('CLOUD_COVER', 50)); 
print('Size landsat8FiltClouds collection', landsat8FiltClouds.size()); 
// Create a composite from the filtered imagery. 
var compositeFiltClouds = landsat8FiltClouds.median().clip(country); 
Map.addLayer(compositeFiltClouds, visParams, 'L8 Composite cloud filter'); 

//QAmask
// Define the cloud mask function. 
function maskSrClouds(image) { 
    // Bit 0 - Fill 
    // Bit 1 - Dilated Cloud 
    // Bit 2 - Cirrus 
    // Bit 3 - Cloud 
    // Bit 4 - Cloud Shadow 
    var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0); 

    var saturationMask = image.select('QA_RADSAT').eq(0); 
    return image.updateMask(qaMask).updateMask(saturationMask); 
} 
// Apply the cloud mask to the collection. 
var landsat8FiltMasked = landsat8FiltClouds.map(maskSrClouds); 
print('Size landsat8FiltClouds collection maskQA', landsat8FiltClouds.size()); 
// Create a composite. 
var landsat8compositeMasked = landsat8FiltMasked.median().clip(country); 
Map.addLayer(landsat8compositeMasked, visParams, 'L8 composite masked'); 


// ---------- Section 2 -----------------
// Define Landsat 7 Level 2, Collection 2, Tier 1 collection. 
 
// Scaling factors for L7. 
function applyScaleFactorsL7(image) { 
    var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2); 
    var thermalBand = image.select('ST_B6').multiply(0.00341802).add(149.0); 
    return image.addBands(opticalBands, null, true) .addBands(thermalBand, null, true); 
} 
// Filter collection, apply cloud mask, and scaling factors. 
var landsat7FiltMasked = landsat7 
    .filterBounds(country) 
    .filterDate(startDate, endDate) 
    .filter(ee.Filter.lessThan('CLOUD_COVER', 50)) 
    .map(maskSrClouds) 
    .map(applyScaleFactorsL7); 
print('Size landsat7FiltMasked collection maskQA', landsat7FiltMasked.size());     
// Create composite. 
var landsat7compositeMasked = landsat7FiltMasked .median() .clip(country); 
Map.addLayer(landsat7compositeMasked, 
    { 
        bands: ['SR_B3', 'SR_B2', 'SR_B1'], 
        min: 0, 
        max: 0.2 
    }, 
    'L7 composite masked'); 


//合并landsat7和landsat8去云后的影像：
//第一步：把landsat7的波段名与landsat8 的波段名统一，landsat8的432波段对应landsat7的321波段。
//第二步：合并数据集，要确保两个数据集的波段名以及波段数量一致，合并时选择波段。   
// Since Landsat 7 and 8 have different band designations, 
// let's create a function to rename L7 bands to match to L8. 
function rename(image) { 
    return image.select( 
        ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B7'], 
        ['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7']); 
} 
// Apply the rename function. 
var landsat7FiltMaskedRenamed = landsat7FiltMasked.map(rename);

// Merge Landsat collections. 
var landsat78 = landsat7FiltMaskedRenamed .merge(landsat8FiltMasked.select( 
    ['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'])) 
    .map(function(img) { 
        return img.toFloat(); 
    }); 
print('Merged collections', landsat78); 
print('Size Merged collections', landsat78.size()); 

// Create Landsat 7 and 8 image composite and add to the Map. 
var landsat78composite = landsat78.median().clip(country); 
Map.addLayer(landsat78composite, visParams, 'L7 and L8 composite'); 

 

三、GEE-BAP 

Best-Available-Pixel Compositing Earth Engine Application

3.1 应用程序访问地址：https://saveriofrancini.users.earthengine.app/view/bap-gee