利用GEE批量下载基于L8影像SR数据集的某一月份的时序RSEI数据

泗水榴莲

已于 2023-09-13 13:40:51 修改

阅读量1k

点赞数 1

分类专栏： GEE学习文章标签： javascript

于 2023-09-01 10:46:02 首次发布

版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。

本文链接：https://blog.csdn.net/qq_33894731/article/details/132612357

版权

GEE学习专栏收录该内容

5 篇文章 3 订阅

订阅专栏

SR数据集属于Landsat8影像的L2级数据，优点是经过了大气处理。

本文基于L8影像的SR数据集，计算了北京市2013年到2022年8月份的RSEI数据。

利用GEE批量下载某一月份时序RSEI数据的具体流程如下：

目录

1. 定义研究区和去云函数

2. 定义SR月份数据和去云处理

3. 计算各个指数（NDVI、WET、NDBSI、LST）并去水体

4. 将各个指数分别归一化

5. 定义PCA函数并计算归一化后的RSEI

6. 批量导出RSEI影像

8. RSEI运行结果

9. 完整代码

1. 定义研究区和去云函数


var roi = table;
Map.addLayer(roi, {}, 'yanjiuqu');

// define cloudMasking function
function maskL8sr(image) {
  var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);
  var saturationMask = image.select('QA_RADSAT').eq(0);

  // Apply the scaling factors to the appropriate bands.
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);

  // Replace the original bands with the scaled ones and apply the masks.
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBands, null, true)
              .updateMask(qaMask)
              .updateMask(saturationMask);
}

2. 定义SR月份数据和去云处理

// define SR dateset
var Imagecol = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
                 .filterBounds(roi)
                 .filterDate('2013-04-08', '2022-12-31')
                 .filterMetadata('CLOUD_COVER', 'less_than', 50)
                 .sort('CLOUD_COVER', true)
                 .map(function(image){
                   return image.set('year', ee.Image(image).date().get('year')).set('month', ee.Image(image).date().get('month'));
                 })
                 .map(maskL8sr);

print('CloudRemoved_Images', Imagecol);

// to prepare mosaiced and clipped dataset
var imgList = ee.List([]);
for (var year = 2013; year < 2023; year++) {
  for (var month = 8; month < 9; month++) {
    var img = Imagecol.filterMetadata('year', 'equals', year)
                      .filterMetadata('month', 'equals', month)
                      .mosaic()
                      .clip(roi);
    var yearmonth = ee.String(ee.Number(year)).cat('-').cat(ee.String(ee.Number(month)));
    var image = img.set('yearmonth', yearmonth);
    var imgList = imgList.add(image);
  }
}

print('Filtered_Images', imgList);

// to filter valid images
var Valid_Imgs = ee.ImageCollection(imgList)
                   .map(function(image) {
                      var bands = image.bandNames();
                      var hasBands = bands.size().gt(0);
                      return ee.Algorithms.If(hasBands, image, null);
                    },true);

print('Valid_Images', Valid_Imgs);

3. 计算各个指数（NDVI、WET、NDBSI、LST）并去水体

// data process
// to define watermask function
function waterMask(image) {
  var ndwi = image.normalizedDifference(['SR_B3', 'SR_B5']);
  var mask = ndwi.lt(0.1);
  return image.updateMask(mask);
}

// to define ndvi function
function getNDVI(image) {
  var ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']);
  return image.addBands(ndvi.rename('NDVI'));
}

// to define wet function
function getWET(image) {
  var WET = image.expression('B*(0.1509) + G*(0.1973) + R*(0.3279) + NIR*(0.3406) + SWIR1*(-0.7112) + SWIR2*(-0.4572)', {
    'B': image.select('SR_B2'),
    'G': image.select('SR_B3'),
    'R': image.select('SR_B4'),
    'NIR': image.select('SR_B5'),
    'SWIR1': image.select('SR_B6'),
    'SWIR2': image.select('SR_B7')
  });
  return image.addBands(WET.rename('WET'));
}

// to define ndbsi function
function getNDBSI(image) {
  // NDBSI = ( IBI + SI ) / 2
  var IBI = image.expression('(2 * SWIR1 / (SWIR1 + NIR) - (NIR / (NIR + RED) + GREEN / (GREEN + SWIR1))) / (2 * SWIR1 / (SWIR1 + NIR) + (NIR / (NIR + RED) + GREEN / (GREEN + SWIR1)))', {
    'SWIR1': image.select('SR_B6'),
    'NIR': image.select('SR_B5'),
    'RED': image.select('SR_B4'),
    'GREEN': image.select('SR_B3')
  });
  var SI = image.expression('((SWIR1 + RED) - (NIR + BLUE)) / ((SWIR1 + RED) + (NIR + BLUE))', {
    'SWIR1': image.select('SR_B6'),
    'NIR': image.select('SR_B5'),
    'RED': image.select('SR_B4'),
    'BLUE': image.select('SR_B2')
  });
  var NDBSI = (IBI.add(SI)).divide(2);
  return image.addBands(NDBSI.rename('NDBSI'));
}

// to dedfine lst function
function getLST(image) {
  var LST = image.expression('(LST - 273.15)', {
    LST: image.select('ST_B10')
  });
  return image.addBands(LST.rename('LST'));
}

// to get NDVI, WET, NDBSI, and LST
var Indexed_Imgs = Valid_Imgs.map(waterMask)
                             .map(getNDVI)
                             .map(getWET)
                             .map(getNDBSI)
                             .map(getLST);

4. 将各个指数分别归一化

var bandNames = ['NDVI', 'WET', 'NDBSI', 'LST'];
var Prepared_Imgs = Indexed_Imgs.select(bandNames);

// to define norm function
function getNorm(image) {
  var minMax = image.reduceRegion({
    reducer: ee.Reducer.minMax(),
    geometry: roi,
    scale: 30,
    crs: 'EPSG:4326',
    bestEffort: true,
    maxPixels: 1e13
  });
  var yearmonth = image.get('yearmonth');
  var normalize  = ee.ImageCollection.fromImages(
      image.bandNames().map(function(name){
            name = ee.String(name);
            var band = image.select(name);
            return band.unitScale(ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max'))));
        })
    ).toBands().set('yearmonth', yearmonth);
    return normalize;
}

var Normed_Imgs = Prepared_Imgs.map(getNorm);

print('Normed_Images', Normed_Imgs);

5. 定义PCA函数并计算归一化后的RSEI

// to define PCA function
var getPCA = function(image){
  
  var bandNames = image.bandNames();
  var region = roi;
  var yearmonth = image.get('yearmonth');
  // Mean center the data to enable a faster covariance reducer
  // and an SD stretch of the principal components.
  var meanDict = image.reduceRegion({
        reducer:  ee.Reducer.mean(),
        geometry: region,
        scale: 30,
        maxPixels: 1e13
  });
  var means = ee.Image.constant(meanDict.values(bandNames));
  var centered = image.subtract(means).set('yearmonth', yearmonth);
  
  
  // This helper function returns a list of new band names.
  var getNewBandNames = function(prefix, bandNames){
    var seq = ee.List.sequence(1, 4);
    //var seq = ee.List.sequence(1, bandNames.length());
    return seq.map(function(n){
      return ee.String(prefix).cat(ee.Number(n).int());
    });      
  };
  
  // This function accepts mean centered imagery, a scale and
  // a region in which to perform the analysis.  It returns the
  // Principal Components (PC) in the region as a new image.
  var getPrincipalComponents = function(centered, scale, region){
    var year = centered.get('yearmonth');
    var arrays = centered.toArray();
    
    // Compute the covariance of the bands within the region.
    var covar = arrays.reduceRegion({
      reducer: ee.Reducer.centeredCovariance(),
      geometry: region,
      scale: scale,
      bestEffort:true,
      maxPixels: 1e13
    });
        
    // Get the 'array' covariance result and cast to an array.
    // This represents the band-to-band covariance within the region.
    var covarArray = ee.Array(covar.get('array'));
        
    // Perform an eigen analysis and slice apart the values and vectors.
    var eigens = covarArray.eigen();

    // This is a P-length vector of Eigenvalues.
    var eigenValues = eigens.slice(1, 0, 1);
    // This is a PxP matrix with eigenvectors in rows.
    var eigenVectors = eigens.slice(1, 1);

    // Convert the array image to 2D arrays for matrix computations.
    var arrayImage = arrays.toArray(1);
    // Left multiply the image array by the matrix of eigenvectors.
    var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);

    // Turn the square roots of the Eigenvalues into a P-band image.
    var sdImage = ee.Image(eigenValues.sqrt())
    .arrayProject([0]).arrayFlatten([getNewBandNames('SD',bandNames)]);

    // Turn the PCs into a P-band image, normalized by SD.
    return principalComponents
    // Throw out an an unneeded dimension, [[]] -> [].
    .arrayProject([0])
    // Make the one band array image a multi-band image, [] -> image.
    .arrayFlatten([getNewBandNames('PC', bandNames)])
    // Normalize the PCs by their SDs.
    .divide(sdImage)
    .set('yearmonth', yearmonth);
  };
    
  // Get the PCs at the specified scale and in the specified region
  image = getPrincipalComponents(centered, 1000, region);
  return image;
};
  
var PCA_Imgcol = Normed_Imgs.map(getPCA);
print('PCA_Imgcol', PCA_Imgcol);

// to compute RSEI
var Normed_RSEI = PCA_Imgcol.map(function(image){
  // image = image.addBands(ee.Image(1).rename('constant'));
  // var rsei = image.expression('constant - pc1' , {
  //     constant: image.select('constant'),
  //     pc1: image.select('PC1')
  // });
  var rsei = image.select('PC1');
  var RSEI = getNorm(rsei);
  return image.addBands(RSEI.rename('RSEI'));
});

print('Normed_RSEI', Normed_RSEI);

6. 批量导出RSEI影像

// data download
// 批量下载函数
function exportImage(image, fileName) {  
    Export.image.toDrive({  
      image: image.select('RSEI'),
      description: 'RSEI-' + fileName,
      folder: 'GEE_export',
      region: roi,
      scale: 30,
      crs: 'EPSG:4326',
      maxPixels: 1e13  //最大像元素
  });  
}

// 生成列表，迭代下载
var indexList = Normed_RSEI.reduceColumns(ee.Reducer.toList(), ['system:index']).get('list'); 
indexList.evaluate(function(indexs) {
  for (var i = 0; i < indexs.length; i++) { 
    var image = ee.Image(Normed_RSEI.filter(ee.Filter.eq('system:index', indexs[i])).first());
    var name = 2013 + i + '-8';
    exportImage(image, name);
  }
});

7. 批量RUN

// 运行结果出来后，将代码粘贴到网页控制台

runTasks = function() {
    const evt = new MouseEvent('click', {bubbles: true, cancelable: true, ctrlKey: true})
    $$('.run-button' ,$$('ee-task-pane')[0].shadowRoot).forEach(function(e) {
        e.dispatchEvent(evt)
    })
}
runTasks()

8. RSEI运行结果

9. 完整代码

// 定义成自己的研究区
var roi = table;

// 在图层上显示
Map.addLayer(roi, {}, 'yanjiuqu');

// 定义SR数据的去云函数
function maskL8sr(image) {
  var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);
  var saturationMask = image.select('QA_RADSAT').eq(0);

  // Apply the scaling factors to the appropriate bands.
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);

  // Replace the original bands with the scaled ones and apply the masks.
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBands, null, true)
              .updateMask(qaMask)
              .updateMask(saturationMask);
}

//_______________________________________________________________________________________
// 定义SR数据集并进行去云处理
var Imagecol = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
                 .filterBounds(roi)
                 .filterDate('2013-04-08', '2022-12-31')
                 .filterMetadata('CLOUD_COVER', 'less_than', 50) // 根据需要替换阈值
                 .sort('CLOUD_COVER', true)
                 .map(function(image){
                   return image.set('year', ee.Image(image).date().get('year')).set('month', ee.Image(image).date().get('month'));
                 })
                 .map(maskL8sr);

print('CloudRemoved_Images', Imagecol);

// 筛选想要的月份数据并进行裁剪
var imgList = ee.List([]);
for (var year = 2013; year < 2023; year++) {
  for (var month = 8; month < 9; month++) { // 根据需要替换月份
    var img = Imagecol.filterMetadata('year', 'equals', year)
                      .filterMetadata('month', 'equals', month)
                      .mosaic()
                      .clip(roi);
    var yearmonth = ee.String(ee.Number(year)).cat('-').cat(ee.String(ee.Number(month)));
    var image = img.set('yearmonth', yearmonth);
    var imgList = imgList.add(image);
  }
}

print('Filtered_Images', imgList);

// 筛选出有效的数据集
var Valid_Imgs = ee.ImageCollection(imgList)
                   .map(function(image) {
                    var bands = image.bandNames();
                    var hasBands = bands.size().gt(0);
                    return ee.Algorithms.If(hasBands, image, null);
                   },true);

print('Valid_Images', Valid_Imgs);

//_______________________________________________________________________________________
// 数据处理
// 定义去水体函数
function waterMask(image) {
  var ndwi = image.normalizedDifference(['SR_B3', 'SR_B5']);
  var mask = ndwi.lt(0.1);
  return image.updateMask(mask);
}

// 定义NDVI函数
function getNDVI(image) {
  var ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']);
  return image.addBands(ndvi.rename('NDVI'));
}

// 定义WET函数
function getWET(image) {
  var WET = image.expression('B*(0.1509) + G*(0.1973) + R*(0.3279) + NIR*(0.3406) + SWIR1*(-0.7112) + SWIR2*(-0.4572)', {
    'B': image.select('SR_B2'),
    'G': image.select('SR_B3'),
    'R': image.select('SR_B4'),
    'NIR': image.select('SR_B5'),
    'SWIR1': image.select('SR_B6'),
    'SWIR2': image.select('SR_B7')
  });
  return image.addBands(WET.rename('WET'));
}

// 定义NDBSI函数
function getNDBSI(image) {
  // NDBSI = ( IBI + SI ) / 2
  var IBI = image.expression('(2 * SWIR1 / (SWIR1 + NIR) - (NIR / (NIR + RED) + GREEN / (GREEN + SWIR1))) / (2 * SWIR1 / (SWIR1 + NIR) + (NIR / (NIR + RED) + GREEN / (GREEN + SWIR1)))', {
    'SWIR1': image.select('SR_B6'),
    'NIR': image.select('SR_B5'),
    'RED': image.select('SR_B4'),
    'GREEN': image.select('SR_B3')
  });
  var SI = image.expression('((SWIR1 + RED) - (NIR + BLUE)) / ((SWIR1 + RED) + (NIR + BLUE))', {
    'SWIR1': image.select('SR_B6'),
    'NIR': image.select('SR_B5'),
    'RED': image.select('SR_B4'),
    'BLUE': image.select('SR_B2')
  });
  var NDBSI = (IBI.add(SI)).divide(2);
  return image.addBands(NDBSI.rename('NDBSI'));
}

// 定义LST函数
function getLST(image) {
  var LST = image.expression('(LST - 273.15)', {
    LST: image.select('ST_B10')
  });
  return image.addBands(LST.rename('LST'));
}

// 计算4个指数并去水体
var Indexed_Imgs = Valid_Imgs.map(waterMask)
                             .map(getNDVI)
                             .map(getWET)
                             .map(getNDBSI)
                             .map(getLST);

//_______________________________________________________________________________________
// 筛选出所需波段
var bandNames = ['NDVI', 'WET', 'NDBSI', 'LST'];
var Prepared_Imgs = Indexed_Imgs.select(bandNames);

// 定义归一化函数
function getNorm(image) {
  var minMax = image.reduceRegion({
    reducer: ee.Reducer.minMax(),
    geometry: roi,
    scale: 30,
    crs: 'EPSG:4326',
    bestEffort: true,
    maxPixels: 1e13
  });
  var yearmonth = image.get('yearmonth');
  var normalize  = ee.ImageCollection.fromImages(
      image.bandNames().map(function(name){
            name = ee.String(name);
            var band = image.select(name);
            return band.unitScale(ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max'))));
        })
    ).toBands().set('yearmonth', yearmonth);
    return normalize;
}

// 对指数分别归一化
var Normed_Imgs = Prepared_Imgs.map(getNorm);

print('Normed_Images', Normed_Imgs);

//_______________________________________________________________________________________
// 定义PCA函数
var getPCA = function(image){
  
  var bandNames = image.bandNames();
  var region = roi;
  var yearmonth = image.get('yearmonth');
  // Mean center the data to enable a faster covariance reducer
  // and an SD stretch of the principal components.
  var meanDict = image.reduceRegion({
        reducer:  ee.Reducer.mean(),
        geometry: region,
        scale: 30,
        maxPixels: 1e13
  });
  var means = ee.Image.constant(meanDict.values(bandNames));
  var centered = image.subtract(means).set('yearmonth', yearmonth);
  
  
  // This helper function returns a list of new band names.
  var getNewBandNames = function(prefix, bandNames){
    var seq = ee.List.sequence(1, 4);
    //var seq = ee.List.sequence(1, bandNames.length());
    return seq.map(function(n){
      return ee.String(prefix).cat(ee.Number(n).int());
    });      
  };
  
  // This function accepts mean centered imagery, a scale and
  // a region in which to perform the analysis.  It returns the
  // Principal Components (PC) in the region as a new image.
  var getPrincipalComponents = function(centered, scale, region){
    var year = centered.get('yearmonth');
    var arrays = centered.toArray();
    
    // Compute the covariance of the bands within the region.
    var covar = arrays.reduceRegion({
      reducer: ee.Reducer.centeredCovariance(),
      geometry: region,
      scale: scale,
      bestEffort:true,
      maxPixels: 1e13
    });
        
    // Get the 'array' covariance result and cast to an array.
    // This represents the band-to-band covariance within the region.
    var covarArray = ee.Array(covar.get('array'));
        
    // Perform an eigen analysis and slice apart the values and vectors.
    var eigens = covarArray.eigen();

    // This is a P-length vector of Eigenvalues.
    var eigenValues = eigens.slice(1, 0, 1);
    // This is a PxP matrix with eigenvectors in rows.
    var eigenVectors = eigens.slice(1, 1);

    // Convert the array image to 2D arrays for matrix computations.
    var arrayImage = arrays.toArray(1);
    // Left multiply the image array by the matrix of eigenvectors.
    var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);

    // Turn the square roots of the Eigenvalues into a P-band image.
    var sdImage = ee.Image(eigenValues.sqrt())
    .arrayProject([0]).arrayFlatten([getNewBandNames('SD',bandNames)]);

    // Turn the PCs into a P-band image, normalized by SD.
    return principalComponents
    // Throw out an an unneeded dimension, [[]] -> [].
    .arrayProject([0])
    // Make the one band array image a multi-band image, [] -> image.
    .arrayFlatten([getNewBandNames('PC', bandNames)])
    // Normalize the PCs by their SDs.
    .divide(sdImage)
    .set('yearmonth', yearmonth);
  };
    
  // Get the PCs at the specified scale and in the specified region
  image = getPrincipalComponents(centered, 1000, region);
  return image;
};
  
var PCA_Imgcol = Normed_Imgs.map(getPCA);
print('PCA_Imgcol', PCA_Imgcol);

// 计算归一化的RSEI
var Normed_RSEI = PCA_Imgcol.map(function(image){

  // 利用'1-PC1'表示RSEI
  // image = image.addBands(ee.Image(1).rename('constant'));
  // var rsei = image.expression('constant - pc1' , {
  //     constant: image.select('constant'),
  //     pc1: image.select('PC1')
  // });
  
  // 利用'PC1'表示RSEI
  var rsei = image.select('PC1');
  var RSEI = getNorm(rsei);
  return image.addBands(RSEI.rename('RSEI'));
});

print('Normed_RSEI', Normed_RSEI);

//_______________________________________________________________________________________
// 数据下载
// 批量下载函数
function exportImage(image, fileName) {  
    Export.image.toDrive({  
      image: image.select('RSEI'),
      description: 'RSEI-' + fileName,
      folder: 'GEE_export',
      region: roi,
      scale: 30,
      crs: 'EPSG:4326',
      maxPixels: 1e13  //最大像元素
  });  
}

// 生成列表，迭代下载
var indexList = Normed_RSEI.reduceColumns(ee.Reducer.toList(), ['system:index']).get('list'); 
indexList.evaluate(function(indexs) {
  for (var i = 0; i < indexs.length; i++) { 
    var image = ee.Image(Normed_RSEI.filter(ee.Filter.eq('system:index', indexs[i])).first());
    var name = 2013 + i + '-8';
    exportImage(image, name);
  }
});

参考了以下文章：

PCA函数参考：利用GEE计算城市遥感生态指数(RSEI)——Landsat 8为例_gee计算rsei_何处惹尘埃？的博客-CSDN博客GEE，RSEI，城市生态状况快速监测与评价_gee计算rseihttps://blog.csdn.net/qq_35490698/article/details/125751534

批量导出参考：GEE学习笔记001-ImageCollection（影像集合）的批量裁剪和批量导出_gee导出imagecollection_小郭学地信的博客-CSDN博客在GEE（Google Earth Engine）中，遥感影像可以是单张的Image，也可以是多张影像组成的影像集合（ImageCollection）。其中，Image的裁剪很简单，直接调用clip函数即可，但是针对ImageCollection，不能直接使用clip函数。_gee导出imagecollectionhttps://blog.csdn.net/guoshihui123/article/details/129242873批量RUN参考：新版GEE批量导出点run代码_P0级推荐算法工程师的博客-CSDN博客偶然在GEE的论坛上看到杨志强老师写的新版GEE批量导出点run代码，分享给大家，杨志强老师还实现了LandTrendr算法（膜拜大神），论坛讨论链接在这儿了。Exporting images without clicking on RUN button (google.com)当导出run的列表全部显示后，打开浏览器的console命令行运行窗口，复制粘贴杨老师的批量点run代码，再enter，耐心等待即可runTasks = function() { const evt = new Mhttps://blog.csdn.net/SunStrongInChina/article/details/120707460

关注

1
点赞
踩
31

收藏

觉得还不错? 一键收藏
3
评论
利用GEE批量下载基于L8影像SR数据集的某一月份的时序RSEI数据

SR数据集属于Landsat8影像的L2级数据，优点是经过了大气处理。本文基于L8影像的SR数据集，计算了北京市2013年到2022年8月份的RSEI数据。利用GEE批量下载某一月份时序RSEI数据的具体流程如下：
复制链接

扫一扫

专栏目录

评论 3

被折叠的条评论为什么被折叠?

到【灌水乐园】发言

查看更多评论

添加红包

成就一亿技术人!

hope_wisdom

发出的红包

实付元

使用余额支付

点击重新获取

扫码支付

钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。