GEE初学实操：将Landsat8 LST导入到Sentinel-2RSEI中

1.下载Sentinel-2数据，这里我们只用到B2,B3,B4,B8,B11,B12六个波段，由于空间分辨率不同，我们要将B11，B12两个波段重采样到10米。

// 定义云掩码函数
function maskS2clouds(image) {
  var qa = image.select('QA60');
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));
  return image.updateMask(mask).divide(10000);
}
// 创建影像集合，过滤时间和云量，应用云掩码
var sentinel2_B2348 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
                      .filterBounds(table)
                      .filterDate('2022-01-01', '2022-12-31')
                      .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
                      .map(maskS2clouds)
                      .median().select('B2','B3','B4','B8')
                      .clip(table);
var sentinel2_B1112 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
                      .filterBounds(table)
                      .filterDate('2022-01-01', '2022-12-31')
                      .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30))
                      .map(maskS2clouds)
                      .median().select('B11','B12')
                      .clip(table);
var sentinel2_B2348_50N = sentinel2_B2348.reproject({
  crs: 'EPSG:32650', // UTM Zone 50N
  scale: 10 // Sentinel-2 的分辨率
});
var sentinel2_B1112_50N = sentinel2_B1112.reproject({
  crs: 'EPSG:32650', // UTM Zone 50N
  scale: 10 // Sentinel-2 的分辨率
});
var sentinel2_50N = ee.Image.cat([sentinel2_B2348_50N, sentinel2_B1112_50N]);

2.下载Landsat8波段，这里我们只需要Landsat8的B10波段，重采样到10米。

// 加载Landsat影像数据集
var landsat8_B10 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
              .filterBounds(table)
              .filterDate('2022-06-01', '2022-08-31')
              .filter(ee.Filter.lte('CLOUD_COVER',10))
              .median().select('ST_B10')
              .clip(table);
var landsat8_B10_50N = landsat8_B10.reproject({
  crs: 'EPSG:32650', // UTM Zone 50N
  scale: 10 // Sentinel-2 的分辨率
});

3.计算LST指数。

// 加载Landsat 8影像集合，过滤时间和云量
var landsat8Collection = ee.ImageCollection('LANDSAT/LC08/C02/T1_RT')
  .filterBounds(table)
  .filterDate('2022-06-01', '2022-08-31')
  .filter(ee.Filter.lt('CLOUD_COVER', 10));
// 获取集合中的第一张影像
var landsat8Image = landsat8Collection.first();
// 定义常数
var K1 = ee.Number(774.8853);  // K1 常数
var K2 = ee.Number(1321.0789); // K2 常数
var L = ee.Number(10.895);     // 波长
var rho = ee.Number(1.438e-2); // 普朗克常数除以波尔兹曼常数
// 从影像元数据中获取乘数和加数
var radMultBand10 = ee.Number(landsat8Image.get('RADIANCE_MULT_BAND_10'));
var radAddBand10 = ee.Number(landsat8Image.get('RADIANCE_ADD_BAND_10'));
// 计算辐射亮度
var radiance = landsat8_B10_50N.expression(
  'ML * DN + AL', {
    'ML': radMultBand10,
    'AL': radAddBand10,
    'DN': landsat8_B10_50N.select('ST_B10')
  }
);
// 计算亮度温度
var brightnessTemp = radiance.expression(
  'K2 / log((K1 / L) + 1)', {
    'K1': K1,
    'K2': K2,
    'L': radiance
  }
);
// 假设地表比辐射率
var emissivity = 0.95;
// 计算地表温度
var lst = brightnessTemp.expression(
  'BT / (1 + (L*BT / rho) * log(e))', {
    'BT': brightnessTemp,
    'L': L,
    'rho': rho,
    'e': emissivity
  }).rename('LST');

4.计算RSEI指数。

// 计算NDVI
var ndvi = sentinel2_50N.normalizedDifference(['B8', 'B4']).rename('NDVI');

// 计算缨帽变换WET
var wet = sentinel2_50N.expression('0.1509 * B2 + 0.1973 * B3 + 0.3279 * B4 + 0.3406 * B8 - 0.7112 * B11 - 0.4572 * B12', 
  {
    'B2': sentinel2_50N.select('B2'),
    'B3': sentinel2_50N.select('B3'),
    'B4': sentinel2_50N.select('B4'),
    'B8': sentinel2_50N.select('B8'),
    'B11': sentinel2_50N.select('B11'),
    'B12': sentinel2_50N.select('B12')
  }).rename('WET');

// 计算SI
var si = sentinel2_50N.expression('((B11 + B4)-(B8 + B2)) / ((B11 + B4)+(B8 + B2))',
  {
    'B2': sentinel2_50N.select('B2'),
    'B4': sentinel2_50N.select('B4'),
    'B8': sentinel2_50N.select('B8'),
    'B11': sentinel2_50N.select('B11')
  }).rename('SI');
// 计算IBI
var ibi = sentinel2_50N.expression('(2 * B11 / (B11 + B8) - (B8 / (B8 + B4) + B3 / (B3 + B11))) / (2 * B11 / (B11 + B8) + (B8 / (B8 + B4) + B3 / (B3 + B11)))',
  {
    'B3': sentinel2_50N.select('B3'),
    'B4': sentinel2_50N.select('B4'),
    'B8': sentinel2_50N.select('B8'),
    'B11': sentinel2_50N.select('B11')
  }).rename('IBI');
// 将 si 和 ibi 添加到影像中
sentinel2_50N = sentinel2_50N.addBands([si, ibi]);
// 计算NDSI
var ndsi = sentinel2_50N.expression('(si + ibi) / 2',
  {
    'si': sentinel2_50N.select('SI'),
    'ibi': sentinel2_50N.select('IBI'),
  }).rename('NDSI');
  
// 将所有指数堆叠为多波段影像
var indices = ee.Image.cat([ndvi, wet, ndsi, lst]);

// 标准化各个指数
  var meanDict = indices.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: table.geometry(),
    scale: 10,
    maxPixels: 1e13
  });
  var stdDevDict = indices.reduceRegion({
    reducer: ee.Reducer.stdDev(),
    geometry: table.geometry(),
    scale: 10,
    maxPixels: 1e13
  });
  
// 遍历每个波段进行标准化
var standardizedIndices = indices.bandNames().map(function(bandName) {
  bandName = ee.String(bandName);
  var band = indices.select(bandName);
  var mean = ee.Number(meanDict.get(bandName));
  var stdDev = ee.Number(stdDevDict.get(bandName));
  return band.subtract(mean).divide(stdDev).rename(bandName);
});
print('standardizedIndices details:', standardizedIndices);

// 将标准化后的波段组合成一个影像
standardizedIndices = ee.ImageCollection.fromImages(standardizedIndices);
print('standardizedIndices details:', standardizedIndices);

  
// 计算协方差矩阵
var covar = standardizedIndices.map(function(image) {
  return image.reduceRegion({
    reducer: ee.Reducer.covariance(),
    geometry: table.geometry(),
    scale: 10,
    maxPixels: 1e13
  });
});
var covarArray = ee.Array(covar.get('array'));
  
// 计算特征值和特征向量
var eigens = covarArray.eigen();

// 获取特征向量
var eigenVectors = eigens.slice(1, 1);

// 计算主成分
var pca = standardizedIndices.map(function(img) {
  return img.toArray().matrixMultiply(eigenVectors).arrayProject([0]).arrayFlatten([['pc1', 'pc2', 'pc3', 'pc4', 'pc5']]);
});
  
// 将计算出的所有主成分图像合并成一个多波段图像
pca = ee.ImageCollection.fromImages(pca.flatten());
  
// 获取第一主成分
var pc1 = pca.limit(1).mosaic();

// 计算RSEI
var minMax = pc1.reduceRegion({
  reducer: ee.Reducer.minMax(),
  geometry: table.geometry(),
  scale: 10,
  maxPixels: 1e13
});
 var min = ee.Number(minMax.get('pc1_min'));
var max = ee.Number(minMax.get('pc1_max'));
var rsei = pc1.subtract(min).divide(max.subtract(min)).rename('RSEI');

// 导出RSEI影像
Export.rsei.toDrive({
  image: rseiImage.select('RSEI'),
  description: 'RSEI_2022',
  scale: 10,
  maxPixels: 1e13,
  region: table,
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true
  }
});

5.查看影像投影及分辨率代码

// 获取影像的投影信息

var projection = image.select('B4').projection(); // 替换为感兴趣的波段

// 打印投影信息

print('Projection:', projection);

// 获取影像的分辨率

var scale = projection.nominalScale();

// 打印分辨率

print('Resolution (scale):', scale);

6.查看像元最大最小值代码

var minMax = image.reduceRegion({
  reducer: ee.Reducer.minMax(),
  geometry: table.geometry(),
  scale: 10,
  maxPixels: 1e9
});
// 打印最小值和最大值
print('Min and Max:', minMax);

7.各个指数可视化

Map.addLayer(ndvi, {min: 0, max: 1, palette: ['blue', 'white', 'green']}, 'NDVI');
Map.addLayer(wet, {min: 0, max: 1, palette: ['blue', 'white', 'green']}, 'WET');
Map.addLayer(ndbsi, {min: 0, max: 1, palette: ['blue', 'white', 'green']}, 'NDSI');
Map.addLayer(si, {min: 0, max: 1, palette: ['blue', 'white', 'green']}, 'SI');
Map.addLayer(ibi, {min: 0, max: 1, palette: ['blue', 'white', 'green']}, 'IBI');
Map.addLayer(lst, {min: 0, max: 1, palette: ['blue', 'white', 'green']}, 'LST');