GEE随记（三）：提取影像的纹理特征

1. 选取感兴趣区域

下面的代码是随机画取一个矩形geometry，并且空心显示

var roi = ee.FeatureCollection(geometry);
Map.centerObject(roi,7);
var styling ={color:'red',fillColor:'00000000',width:2};// display hollow roi
Map.addLayer(roi.style(styling), {}, "outline");

2. 获取矢量范围内的遥感影像，以Landsat 8为例

下面的代码是获取一定时间空间范围内的Landsat 8 影像集ImageCollection，去云后中值合成为一张Image，并将波段名字从 “ SR_B* ” 改为 “ B* ”

//设置去云影像函数(SR)
//包括云周边、卷云、云、云阴影和雪
function masksr(image){
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
  var thermalBand = image.select('ST_B.*').multiply(0.00341802).add(149.0)
  var qa = image.select('QA_PIXEL')
  var CloudBit = 1<<3
  var CloudShadowBit = 1<<4;
  var DilatedCloudBit = 1<<1;
  var CirrusBit =1<<2;
  // 明确条件，设置两个值都为0
  var mask = qa.bitwiseAnd(CloudBit).eq(0)
               .and(qa.bitwiseAnd(CloudShadowBit).eq(0))
               .and(qa.bitwiseAnd(DilatedCloudBit).eq(0))
               .and(qa.bitwiseAnd(CirrusBit).eq(0));
  return image.addBands(opticalBands,null,true)
              .addBands(thermalBand,null,true)
              .updateMask(mask)}
//加载SR影像whole year
var l8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
           .filterBounds(roi)
           .filterDate('2019-01-01','2020-01-01')
           .map(masksr)
var l8_1 = l8.median().clip(roi).select('SR_B[1-7]').rename(['B1','B2','B3','B4','B5','B6','B7'])
Map.addLayer(l8_1.select(['B4','B3','B2']),{max:0.5,gamma:0.8},'rgb')

 3. 计算灰度波段并计算纹理特征

（1）参考文献：Object-Oriented LULC Classification in Google Earth Engine Combining SNIC, GLCM, and Machine Learning Algorithms；在利用GEE云平台调用.glcmTexture()进行纹理特征计算前需要先计算灰度图层：

（2）计算完灰度图层后可以直接调用.glcmTexture()函数：

// compute glcm texture feature
var gray = l8_1.expression(
      '(0.3 * NIR) + (0.59 * R) + (0.11 * G)', {
      'NIR': l8_1.select('B5'),
      'R': l8_1.select('B4'),
      'G': l8_1.select('B3'),
}).rename('gray');

var glcm = gray.unitScale(0,0.30).multiply(100).toInt().glcmTexture({size: 1,kernel:null});
print('glcm',glcm)

 

 （3）得到的是包含18个波段的影像，根据个人需求或参考常用灰度共生矩阵纹理特征选择这个Image中的波段，再进行归一化处理（不过这一步感觉不做也行，具体可以看分类后精度有无提升）：

// Normalized
var band_glcm = glcm.select(['gray_asm','gray_contrast','gray_corr','gray_var','gray_idm',
                             'gray_savg','gray_svar','gray_sent','gray_ent']);
// calculate the min and max value of an image
var minMax = band_glcm.reduceRegion({
  reducer: ee.Reducer.minMax(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e16,
});

var glcm = ee.ImageCollection.fromImages(
  band_glcm.bandNames().map(function(name){
    name = ee.String(name);
    var band = band_glcm.select(name);
    return band.unitScale(ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max'))))
})).toBands().rename(band_glcm.bandNames());

Map.addLayer(glcm.select('gray_asm'),{max:1})

 

4. 将纹理特征加入波段

一般地物分类而言需要多个波段特征，这时在要分类的Image后直接加.addBands(glcm)，就等于把影像的纹理特征也合并到Image当中了，如：

//add texture feature into image
var data = l8_1.addBands(glcm)
print('data',data)