GEE处理哨兵数据并导出

燕南路GISer
已于 2024-04-26 00:20:06 修改
阅读量461
点赞数
文章标签：数据库前端 javascript
于 2024-04-26 00:03:40 首次发布
本文链接：https://blog.csdn.net/xza13155/article/details/138203116
版权
https://code.earthengine.google.com/dc35e6894b8d75711ca0ec449edc0f73

// 定义研究区域和样本数据
var roi = ee.FeatureCollection('users/messi513112/xiamen2');
var samples = ee.FeatureCollection('users/messi513112/Expor_xiament_point'); // 替换为您的shapefile文件的实际路径

//导入sentinel-1和sentinel-2的图像集合
var s1 = ee.ImageCollection('COPERNICUS/S1_GRD');
var s2 = ee.ImageCollection('COPERNICUS/S2_SR');
// 获取NASA数据
var nasaDEM = ee.Image('NASA/NASADEM_HGT/001').clip(roi);


//重采样
var scale = 10;
// 重采样Sentinel-1影像集
var resampledS1 = s1.map(function(image) {
  return image.resample('bilinear').reproject({
    crs: image.select(0).projection().crs(),
    scale: scale
  });
});

// 重采样Sentinel-2影像集
var resampledS2 = s2.map(function(image) {
  return image.resample('bilinear').reproject({
    crs: image.select(0).projection().crs(),
    scale: scale
  });
});

// 重采样NASA DEM数据
var resampledDEM = nasaDEM.resample('bilinear').reproject({
  crs: nasaDEM.projection().crs(),
  scale: scale
});





// 定义时间范围
var start = ee.Date('2018-01-01');
var end = ee.Date('2018-12-31');

// 计算属性 'YSSZ_S' 中每个值的数量。
var valueCounts = samples.aggregate_histogram('YSSZ_S');
// 打印结果。
print(valueCounts);

// 使用简单的云和云影掩膜
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10和11是云和阴影的标识符
  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 vh = resampledS1
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
  .filter(ee.Filter.eq('instrumentMode', 'IW'))
 .map(function(image) {
          var edge = image.lt(-30.0);
          var maskedImage = image.mask().and(edge.not());
          return image.updateMask(maskedImage).clip(roi);
        });
// 筛选不同观测角度的图像
var vhAscending = vh.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'));
var vhDescending = vh.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));
var spring = ee.Filter.date('2018-01-01', '2018-12-31');
// 从不同极化和观测角度的平均值创建合成图像
var composite = ee.Image.cat([
  ee.ImageCollection(vhAscending.select('VH').merge(vhDescending.select('VH'))).filter(spring).mean(),
  ee.ImageCollection(vhAscending.select('VV').merge(vhDescending.select('VV'))).filter(spring).mean(),
]).focal_median();
print('composite',composite);
// 作为极化和后向散射特性的合成显示
Map.addLayer(composite, {min: -20, max: 0}, 'composite');
// 计算纹理特征
var glcm = composite.multiply(100).toInt().glcmTexture({
    'size': 1,
    'kernel': null
});

// 打印纹理特征
print('GLCM:', glcm);

glcm =glcm.select(['VH_asm',
                   'VH_contrast',
                   'VH_corr',
                   'VH_var',
                   'VH_idm',
                   'VH_savg',
                   'VH_sent',
                   'VH_ent',
                   'VH_dvar',
                   'VH_dent',
                   'VH_imcorr1',
                   'VH_imcorr2',
                   'VH_diss',
                   'VH_inertia',
                   'VH_shade',
                   'VH_prom',
                   //'VV_asm',
                   //'VV_contrast',
                  // 'VV_corr',
                   //'VV_var',
                   //'VV_idm',
                   //'VV_savg',
                   //'VV_sent',
                   //'VV_ent',
                   //'VV_dvar',
                   //'VV_dent',
                   //'VV_imcorr1',
                   //'VV_imcorr2',
                   //'VV_diss',
                   //'VV_inertia',
                   //'VV_shade',
                   //'VV_prom',
                    ]);



// 筛选Sentinel-2数据
var s2filtered = resampledS2
  .filterBounds(roi)
  .filterDate(start, end)
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 40));
// 应用去云函数到筛选后的影像集
var s2CloudRemoved = s2filtered.map(maskS2clouds);
// 打印筛选后的Sentinel-2数据集的元数据，以检查数据集
print('Sentinel-2 filtered metadata:', s2CloudRemoved);
//多时相合成
var s2composite = s2CloudRemoved
  .median() // 多时相合成，这里使用中位数
  .clip(roi); // 裁剪到研究区域
// 显示Sentinel-2多时相合成的结果
Map.addLayer(s2composite, {min: 0, max: 3000, bands: ['B4', 'B3', 'B2']}, 'Sentinel-2 Composite', false);
//提取Sentinel-2的光谱特征
var spectral = s2composite.select(['B4', 'B5', 'B6', 'B7', 'B8']);



// 定义一个函数来计算 NDVI
function addNDVI(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI');
  return ndvi;
}

// 定义一个函数来计算 EVI
function addEVI(image) {
  var evi = image.expression(
    '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
      'NIR': image.select('B8'),
      'RED': image.select('B4'),
      'BLUE': image.select('B2')
    }).rename('EVI');
  return evi;
}

// 定义一个函数来计算 NDWI
function addNDWI(image) {
  var ndwi = image.normalizedDifference(['B3', 'B8']).rename('NDWI');
  return ndwi;
}

// 定义一个函数来计算 NDRE
function addNDRE(image) {
  var ndre = image.normalizedDifference(['B8', 'B5']).rename('NDRE');
  return ndre;
}

// 定义一个函数来计算 SAVI
function addSAVI(image) {
  var savi = image.expression(
    '((NIR - RED) / (NIR + RED + L)) * (1 + L)', {
      'NIR': image.select('B8'),
      'RED': image.select('B4'),
      'L': 0.5
    }).rename('SAVI');
  return savi;
}

// 定义一个函数来计算 MSAVI
function addMSAVI(image) {
  var msavi = image.expression(
    '(2 * NIR + 1 - sqrt((2 * NIR + 1)**2 - 8 * (NIR - RED))) / 2', {
      'NIR': image.select('B8'),
      'RED': image.select('B4')
    }).rename('MSAVI');
  return msavi;
}

// 定义一个函数来计算 ARVI
function addARVI(image) {
  var arvi = image.expression(
    '(NIR - (2 * RED - BLUE)) / (NIR + (2 * RED - BLUE))', {
      'NIR': image.select('B8'),
      'RED': image.select('B4'),
      'BLUE': image.select('B2')
    }).rename('ARVI');
  return arvi;
}




// 应用函数并显示结果
var NDVI = addNDVI(s2composite);
var EVI = addEVI(s2composite);
var NDWI = addNDWI(s2composite);
var NDRE = addNDRE(s2composite);
var SAVI = addSAVI(s2composite);
var MSAVI = addMSAVI(s2composite);
var ARVI = addARVI(s2composite);

// 生长季开始（SOS）
var sos = NDVI.reduce(ee.Reducer.firstNonNull());

// 生长季结束（EOS）
var eos = NDVI.reduce(ee.Reducer.lastNonNull());

// 初花期
var startOfFlowering = NDVI.reduce(ee.Reducer.first().setOutputs(['StartOfFlowering']));

// 叶变色期
var endOfSenescence = NDVI.reduce(ee.Reducer.last().setOutputs(['EndOfSenescence']));

// 准备Sentinel-2数据
var ss2 = ee.ImageCollection('COPERNICUS/S2')
          .filterDate('2018-01-01', '2018-12-31') // 根据设置的时间范围过滤影像
          .filterBounds(roi) // 根据定义的区域过滤影像
          .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 40)) // 过滤云量小于40%的影像
          .map(function(image) {
            // 计算NDVI并将其作为一个波段添加
            var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI');
            // 将system:time_start属性转换为图像并重命名为'time'
            var time = ee.Image.constant(image.get('system:time_start')).toFloat().rename('time');
            // 返回包含NDVI和时间波段的图像
            return ndvi.addBands(time);
          });

// 对ImageCollection应用linearFit reducer
var linearFit = ss2.select(['NDVI', 'time']).reduce(ee.Reducer.linearFit());



// 计算地形特征
var terrain = ee.Algorithms.Terrain(resampledDEM);

// 提取地形特征波段
var elevation = terrain.select('elevation'); // 海拔
var slope = terrain.select('slope'); // 坡度
var aspect = terrain.select('aspect'); // 坡向
var hillshade = terrain.select('hillshade'); // 山体阴影



var img1 = ee.ImageCollection('COPERNICUS/S2')
    .filterDate('2018-2-13', '2018-2-14')
    .filterBounds(roi)
    .median();
  
var img6 = ee.ImageCollection('COPERNICUS/S2')
    .filterDate('2018-9-26', '2018-9-27')
    .filterBounds(roi)
    .median();

// 计算SVI
var svi = img1.addBands(img6)
            .rename(['B4_1', 'B8_1', 'B4_6', 'B8_6'])
            .expression(
              '((1/6) * timeDiff * (LB8 - LB4) * (B4_1 + 2 * B8_1 + 2 * B4_6 + B8_6))/100000', {
                'timeDiff': 225,
                //'centralWavelengths': centralWavelengths,
                'B4_1': img1.select('B4'),
                'B8_1': img1.select('B8'),
                'B4_6': img6.select('B4'),
                'B8_6': img6.select('B8'),
                'LB8':835.1,
                'LB4':664.5
              }
            ).rename('SVI');
            
// 计算B3与B8之间的SVI
var sviB3B8 = img1.addBands(img6)
  .rename(['B3_1', 'B8_1', 'B3_6', 'B8_6'])
  .expression(
    '(((1/6) * timeDiff * (835.1 - 560.0) * (B3_1 + 2 * B8_1 + 2 * B3_6 + B8_6))/10000)-6800', {
      'timeDiff': 225,
      'B3_1': img1.select('B3'),
      'B8_1': img1.select('B8'),
      'B3_6': img6.select('B3'),
      'B8_6': img6.select('B8')
    }
  ).rename('SVI_B3_B8');

// 计算B5与B8之间的SVI
var sviB5B8 = img1.addBands(img6)
  .rename(['B5_1', 'B8_1', 'B5_6', 'B8_6'])
  .expression(
    '((1/6) * timeDiff * (835.1 - 703.9) * (B5_1 + 2 * B8_1 + 2 * B5_6 + B8_6))/10000-2700', {
      'timeDiff': 225,
      'B5_1': img1.select('B5'),
      'B8_1': img1.select('B8'),
      'B5_6': img6.select('B5'),
      'B8_6': img6.select('B8')
    }
  ).rename('SVI_B5_B8');

// 计算B6与B8之间的SVI
var sviB6B8 = img1.addBands(img6)
  .rename(['B6_1', 'B8_1', 'B6_6', 'B8_6'])
  .expression(
    '((1/6) * timeDiff * (835.1 - 740.2) * (B6_1 + 2 * B8_1 + 2 * B6_6 + B8_6))/10000-2700', {
      'timeDiff': 225,
      'B6_1': img1.select('B6'),
      'B8_1': img1.select('B8'),
      'B6_6': img6.select('B6'),
      'B8_6': img6.select('B8')
    }
  ).rename('SVI_B6_B8');

// 计算B7与B8之间的SVI
var sviB7B8 = img1.addBands(img6)
  .rename(['B7_1', 'B8_1', 'B7_6', 'B8_6'])
  .expression(
    '((1/6) * timeDiff * (835.1 - 782.5) * (B7_1 + 2 * B8_1 + 2 * B7_6 + B8_6))/10000-1000', {
      'timeDiff': 225,
      'B7_1': img1.select('B7'),
      'B8_1': img1.select('B8'),
      'B7_6': img6.select('B7'),
      'B8_6': img6.select('B8')
    }
  ).rename('SVI_B7_B8');

// 计算B8A与B8之间的SVI
var sviB8AB8 = img1.addBands(img6)
  .rename(['B8A_1', 'B8_1', 'B8A_6', 'B8_6'])
  .expression(
    '((1/6) * timeDiff * (864.8 - 835.1) * (B8A_1 + 2 * B8_1 + 2 * B8A_6 + B8_6))/10000', {
      'timeDiff': 225,
      'B8A_1': img1.select('B8A'),
      'B8_1': img1.select('B8'),
      'B8A_6': img6.select('B8A'),
      'B8_6': img6.select('B8')
    }
  ).rename('SVI_B8A_B8');

// 特征融合 想导出的特征
var features = glcm//.addBands(spectral)
                   .addBands(NDVI)
                   .addBands(EVI)
                   //.addBands(NDWI)
                   //.addBands(NDRE)
                   //.addBands(SAVI)
                   //.addBands(MSAVI)
                   //.addBands(ARVI)
                   //.addBands(sos)
                   //.addBands(eos)
                   //.addBands(startOfFlowering)
                   //.addBands(endOfSenescence)
                   //.addBands(linearFit)
                   //.addBands(elevation)
                   //.addBands(slope)
                   //.addBands(aspect)
                   //.addBands(hillshade)
                   .addBands(svi)
                   //.addBands(sviB3B8)
                   //.addBands(sviB5B8)
                   //.addBands(sviB6B8)
                   //.addBands(sviB7B8)
                   //.addBands(sviB8AB8);
                   
                   
print(features);

// 假设'image'是您要导出的影像
var imageToExport = features.toFloat(); // 将所有波段转换为Float32类型

// 定义导出参数
var exportParams = {
  image: imageToExport,
  description: 'export_image',
  scale: 10,
  region: roi, // 确保已定义导出区域的几何形状
  maxPixels: 1e13
};

// 开始导出图像
Export.image.toDrive(exportParams);


// 定义您想要保留的样本数量。
var numSamples =200;
// 从特征集合中随机选择样本。
var randomSamples = samples.randomColumn('random').sort('random').limit(numSamples);
print(randomSamples);

// 训练数据准备
var training = features.sampleRegions({
  collection: randomSamples,
  properties: ['YSSZ_S'],
  scale: 10,
  tileScale:16
})


// 打印训练数据集，以检查样本和特征
print('Training data:', training);
// 随机化样本数据集并添加随机列
var samplesWithRandom = training.randomColumn('random');
// 设置训练和验证样本的比例
var split = 0.5; 
var trainingPartition = samplesWithRandom.filter(ee.Filter.lt('random', split));
var testingPartition = samplesWithRandom.filter(ee.Filter.gte('random', split));
// 训练随机森林分类器
var classifier = ee.Classifier.smileRandomForest(50).train({
  features: trainingPartition,
  classProperty: 'YSSZ_S',
  inputProperties: features.bandNames()
});
print(classifier.explain()); 
// 对整个区域进行分类
var classified = features.classify(classifier);

// 获取分类结果的特征集
var classifiedFeatures = ee.FeatureCollection(classified);

// 使用验证样本集评估模型性能
var test = testingPartition.classify(classifier);
// 计算混淆矩阵
var confusionMatrix = test.errorMatrix('YSSZ_S', 'classification');
// 打印混淆矩阵和精度评估指标
print('Confusion Matrix:', confusionMatrix);
print('Overall Accuracy:', confusionMatrix.accuracy());
print('Kappa Accuracy:', confusionMatrix.kappa());

// 定义一个颜色列表，每个类别一个颜色
var palette = [
  'ffd966', // 浅黄色，代表针叶林
  '8db360', // 橄榄绿，代表落叶阔叶林
  '006400', // 深绿色，代表常绿阔叶林
  'ffbb22', // 金色，代表混交林
  'f13c20', // 天蓝色，代表红树林
  '92cddc'  // 地红色，代表灌木丛
];
// 添加分类结果的图层到地图上
Map.addLayer(classified, {min: 1, max: 8, palette: palette}, 'Classification Result');

// 导出分类结果
Export.image.toDrive({
  image: classified,
  description: 'xiamen_classification',
  folder: 'xiamen_classification',
  scale: 10,
  region: roi,
  maxPixels: 1e13
});