基于GEE提取站点NDVI序列

// ===============================================
// Step 1: 定义光伏站点的矢量数据
// ===============================================
print('Step 1: 加载光伏站点矢量数据');
var region = ee.FeatureCollection("projects/ee-jiag/assets/InnerMongolia_station_120_180");
print('光伏站点数量:', region.size());

// 检查属性名称
print('属性名称列表:', region.first().propertyNames());
print('第一个要素的属性:', region.first());

// ===============================================
// Step 2: 定义年份和月份范围
// ===============================================
print('Step 2: 定义年份和月份范围');
var startYear = 2014;
var endYear = 2023;
var startMonth = 1;  // 2月
var endMonth = 12;   // 11月
print('年份范围:', startYear, '到', endYear);
print('月份范围:', startMonth, '到', endMonth);

// ===============================================
// Step 3: 创建年份和月份列表
// ===============================================
print('Step 3: 创建年份和月份列表');
var years = ee.List.sequence(startYear, endYear);
var months = ee.List.sequence(startMonth, endMonth);
print('年份列表:', years);
print('月份列表:', months);

// ===============================================
// Step 4: 定义时间范围的起止日期
// ===============================================
print('Step 4: 定义时间范围的起止日期');
var startDate = ee.Date.fromYMD(startYear, startMonth, 1);
var endDate = ee.Date.fromYMD(endYear, endMonth, 30);
print('起始日期:', startDate.format('YYYY-MM-dd'));
print('结束日期:', endDate.format('YYYY-MM-dd'));

// ===============================================
// Step 5: 加载并预处理 Landsat 8 影像集合
// ===============================================
print('Step 5: 加载并预处理 Landsat 8 影像集合');
var landsatCollection = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
  .filterBounds(region)
  .filterDate(startDate, endDate)
  // 过滤云覆盖率（可选）
  .filter(ee.Filter.lt('CLOUD_COVER_LAND', 20))
  // 选择质量较高的像素并计算 NDVI
  .map(function(image) {
    // 应用缩放因子和偏移量
    var opticalBands = image.select(['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'])
                            .multiply(0.0000275).add(-0.2);
    image = image.addBands(opticalBands, null, true);
    
    // 计算 NDVI
    var ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');
    
    // 添加日期属性
    ndvi = ndvi.set('system:time_start', image.get('system:time_start'));
    
    // 应用 NDVI 有效性掩膜
    return ndvi.updateMask(ndvi.gt(-1).and(ndvi.lt(1)));
  });

print('初始 Landsat 影像数量:', landsatCollection.size());

// ===============================================
// Step 6: 验证影像集合中的日期
// ===============================================
print('Step 6: 验证影像集合中的日期');
var imageDates = landsatCollection.aggregate_array('system:time_start')
  .map(function(time) {
    return ee.Date(time).format('YYYY-MM-dd');
  });
print('过滤后的影像日期:', imageDates);

// ===============================================
// Step 7: 定义函数来提取每个站点的 NDVI 时间序列
// ===============================================
print('Step 7: 定义函数来提取每个站点的 NDVI 时间序列');
var extractNDVIForSites = function(image) {
  // 获取影像日期
  var date = ee.Date(image.get('system:time_start'));
  var year = date.get('year');
  var dateString = date.format('YYYY-MM-dd');

  // 计算该影像中每个站点的 NDVI 统计值（mean, median, max, stdDev）
  var ndviBySite = image.reduceRegions({
    collection: region,
    reducer: ee.Reducer.mean().setOutputs(['NDVI_mean'])
               .combine(
                 ee.Reducer.median().setOutputs(['NDVI_median']),
                 '', 
                 true
               )
               .combine(
                 ee.Reducer.max().setOutputs(['NDVI_max']),
                 '', 
                 true
               )
               .combine(
                 ee.Reducer.stdDev().setOutputs(['NDVI_stdDev']),
                 '', 
                 true
               ),
    scale: 30,        // Landsat 8 的空间分辨率为30米
    tileScale: 4      // 性能优化参数
  });

  // 为每个站点添加日期和年份属性，并过滤掉没有统计值的要素
  ndviBySite = ndviBySite.map(function(feature) {
    return feature.set({
      'date': dateString,
      'year': year,
    });
  }).filter(ee.Filter.notNull(['NDVI_mean', 'NDVI_median', 'NDVI_max', 'NDVI_stdDev']));

  return ndviBySite;
};

// ===============================================
// Step 8: 对影像集合中的每个影像应用 NDVI 提取
// ===============================================
print('Step 8: 对影像集合中的每个影像应用 NDVI 提取');
var ndviTimeSeries = landsatCollection.map(extractNDVIForSites).flatten();
print('NDVI 时间序列总要素数量:', ndviTimeSeries.size());

// ===============================================
// Step 9: 选择需要的属性
// ===============================================
print('Step 9: 选择需要的属性');
// 使用正确的属性名称“编号”
var ndviResults = ndviTimeSeries.select(['编号', 'year', 'date', 'NDVI_mean', 'NDVI_median', 'NDVI_max', 'NDVI_stdDev']);
print('选择后的属性示例:', ndviResults.first());

// ===============================================
// Step 10: 打印部分结果以验证
// ===============================================
print('Step 10: 打印部分 NDVI 时间序列结果');
print('NDVI Time Series for all sites (前 10 个):', ndviResults.limit(10));

// ===============================================
// Step 11: 导出 NDVI 时间序列到 Google Drive
// ===============================================
print('Step 11: 导出 NDVI 时间序列到 Google Drive');
Export.table.toDrive({
  collection: ndviResults,
  description: 'Landsat8_NDVI_2014_2023_1_12_Station120_180',
  fileFormat: 'CSV',
  selectors: ['编号', 'year', 'date', 'NDVI_mean', 'NDVI_median', 'NDVI_max', 'NDVI_stdDev']  // 导出编号、年份、日期和 NDVI 统计值
});
print('导出任务已提交。请在“任务”选项卡中监控进度。');

第二种

输出了每个像元

// ===============================================
// Step 1: 定义光伏站点的矢量数据
// ===============================================
print('Step 1: 加载光伏站点矢量数据');
var region = ee.FeatureCollection("projects/ee-jiag/assets/InnerMongolia_station");
print('光伏站点数量:', region.size());

// 检查属性名称
print('属性名称列表:', region.first().propertyNames());
print('第一个要素的属性:', region.first());

// ===============================================
// Step 2: 定义年份和月份范围
// ===============================================
print('Step 2: 定义年份和月份范围');
var startYear = 2014;
var endYear = 2023;
var startMonth = 1;  // 1月
var endMonth = 12;   // 12月
print('年份范围:', startYear, '到', endYear);
print('月份范围:', startMonth, '到', endMonth);

// ===============================================
// Step 3: 定义时间范围的起止日期
// ===============================================
print('Step 3: 定义时间范围的起止日期');
var startDate = ee.Date.fromYMD(startYear, startMonth, 1);
var endDate = ee.Date.fromYMD(endYear, endMonth, 30);
print('起始日期:', startDate.format('YYYY-MM-dd'));
print('结束日期:', endDate.format('YYYY-MM-dd'));

// ===============================================
// Step 4: 定义云和云影掩膜函数
// ===============================================
print('Step 4: 定义云和云影掩膜函数');

// 云掩膜函数，使用 QA_PIXEL 波段
function maskL8sr(image) {
  var qa = image.select('QA_PIXEL');
  // 位掩码参数，根据官方文档定义
  var cloudBitMask = 1 << 3;
  var cloudShadowBitMask = 1 << 4;

  // 将质量波段的位掩码应用到数据中
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
              .and(qa.bitwiseAnd(cloudShadowBitMask).eq(0));

  // 应用掩膜并返回影像
  return image.updateMask(mask);
}

// ===============================================
// Step 5: 加载并预处理 Landsat 8 影像集合
// ===============================================
print('Step 5: 加载并预处理 Landsat 8 影像集合');

var landsatCollection = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
  .filterBounds(region)
  .filterDate(startDate, endDate)
  // 过滤云覆盖率（可选）
  .filter(ee.Filter.lt('CLOUD_COVER_LAND', 20))
  // 应用云和云影掩膜函数
  .map(maskL8sr)
  // 选择质量较高的像素并计算 NDVI
  .map(function(image) {
    // 应用缩放因子和偏移量
    var opticalBands = image.select(['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'])
                            .multiply(0.0000275).add(-0.2);
    image = image.addBands(opticalBands, null, true);
    
    // 计算 NDVI
    var ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');
    
    // 添加日期属性
    ndvi = ndvi.set('system:time_start', image.get('system:time_start'));
    
    // 应用 NDVI 有效性掩膜
    return ndvi.updateMask(ndvi.gt(-1).and(ndvi.lt(1)));
  });

print('初始 Landsat 影像数量:', landsatCollection.size());

// ===============================================
// Step 6: 定义函数来提取每个站点的每个像元的 NDVI 值
// ===============================================
print('Step 6: 定义函数来提取每个站点的每个像元的 NDVI 值');

var extractNDVIForPixels = function(image) {
  // 获取影像日期
  var date = ee.Date(image.get('system:time_start')).format('YYYY-MM-dd');
  var timestamp = image.get('system:time_start');

  // 对于每个站点，提取其内的所有像元的 NDVI 值
  var ndviByPixel = region.map(function(feature) {
    // 获取站点编号
    var stationId = feature.get('编号');

    // 在站点范围内对 NDVI 影像进行采样，获取每个像元的 NDVI 值
    var pixels = image.sample({
      region: feature.geometry(),
      scale: 30,
      projection: 'EPSG:4326',
      geometries: true  // 保留像元的几何信息
    }).map(function(pixel) {
      // 为每个像元添加站点编号和日期属性
      return pixel.set({
        'station_id': stationId,
        'date': date,
        'timestamp': timestamp
      });
    });

    return pixels;
  }).flatten();  // 将结果展平成一个 FeatureCollection

  return ndviByPixel;
};

// ===============================================
// Step 7: 对影像集合中的每个影像应用 NDVI 像元提取
// ===============================================
print('Step 7: 对影像集合中的每个影像应用 NDVI 像元提取');

var ndviPixelTimeSeries = landsatCollection.map(extractNDVIForPixels).flatten();
print('NDVI 像元时间序列总要素数量:', ndviPixelTimeSeries.size());

// ===============================================
// Step 8: 选择需要的属性
// ===============================================
print('Step 8: 选择需要的属性');

// 选择需要的属性，包括站点编号、日期、像元坐标和 NDVI 值
var ndviResults = ndviPixelTimeSeries.select(['station_id', 'date', 'NDVI', '.geo']);
print('选择后的属性示例:', ndviResults.first());

// ===============================================
// Step 9: 打印部分结果以验证
// ===============================================
print('Step 9: 打印部分 NDVI 像元时间序列结果');
print('NDVI Pixel Time Series for all sites (前 10 个):', ndviResults.limit(10));

// ===============================================
// Step 10: 导出 NDVI 像元时间序列到 Google Drive
// ===============================================
print('Step 10: 导出 NDVI 像元时间序列到 Google Drive');

// 导出时，需要注意数据量可能很大，可以考虑导出为 GeoJSON 格式
Export.table.toDrive({
  collection: ndviResults,
  description: 'Landsat8_NDVI_Pixel_2014_2023_Station',
  fileFormat: 'CSV',
  selectors: ['station_id', 'date', 'NDVI', '.geo']  // 导出站点编号、日期、NDVI 值和几何信息
});
print('导出任务已提交。请在“任务”选项卡中监控进度。');

IQR异常值过滤结果

// ===============================================
// Step 1: Load photovoltaic station vector data
// ===============================================
print('Step 1: Load photovoltaic station vector data');
var region = ee.FeatureCollection("projects/ee-jiag/assets/InnerMongolia_station");
print('Number of photovoltaic stations:', region.size());

// Check property names
print('List of property names:', region.first().propertyNames());
print('First feature properties:', region.first());

// ===============================================
// Step 2: Define time range
// ===============================================
print('Step 2: Define time range');
var startDate = ee.Date('2014-04-01');
var endDate = ee.Date('2023-12-31');
print('Start date:', startDate.format('YYYY-MM-dd'));
print('End date:', endDate.format('YYYY-MM-dd'));

// ===============================================
// Step 3: Define cloud and cloud shadow mask function
// ===============================================
print('Step 3: Define cloud and cloud shadow mask function');

// Cloud mask function using QA_PIXEL band
function maskL8sr(image) {
  var qa = image.select('QA_PIXEL');
  // Bitmask parameters as per official documentation
  var cloudBitMask = 1 << 3;
  var cloudShadowBitMask = 1 << 4;

  // Apply the quality band bitmask to data
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
              .and(qa.bitwiseAnd(cloudShadowBitMask).eq(0));

  // Apply mask and return the image
  return image.updateMask(mask);
}

// ===============================================
// Step 4: Load and preprocess Landsat 8 image collection
// ===============================================
print('Step 4: Load and preprocess Landsat 8 image collection');

var landsatCollection = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
  .filterBounds(region)
  .filterDate(startDate, endDate)
  // Filter cloud cover (optional)
  .filter(ee.Filter.lt('CLOUD_COVER_LAND', 20))
  // Apply cloud and cloud shadow mask function
  .map(maskL8sr)
  // Apply scaling factors and offsets, and compute NDVI
  .map(function(image) {
    // Apply scaling factors and offsets
    var opticalBands = image.select(['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'])
                            .multiply(0.0000275).add(-0.2);
    image = image.addBands(opticalBands, null, true);

    // Compute NDVI
    var ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');

    // Add date property
    var date = ee.Date(image.get('system:time_start')).format('YYYY-MM-dd');
    ndvi = ndvi.set('date', date);

    // Apply NDVI validity mask
    return ndvi.updateMask(ndvi.gt(-1).and(ndvi.lt(1)));
  });

print('Initial number of Landsat images:', landsatCollection.size());

// ===============================================
// Step 5: Define function to extract NDVI values and perform IQR outlier removal
// ===============================================
print('Step 5: Define function to extract NDVI values and perform IQR outlier removal');

var extractNDVIForImage = function(image) {
  var date = image.get('date');
  var timestamp = image.get('system:time_start');

  // For each station, extract all pixel NDVI values
  var statsByStation = region.map(function(feature) {
    // Get station ID
    var stationId = feature.get('编号');

    // Sample NDVI values within the station geometry
    var pixels = image.sample({
      region: feature.geometry(),
      scale: 30,
      projection: 'EPSG:4326',
      geometries: false  // No need for geometry information
    }).aggregate_array('NDVI');

    var ndviValues = ee.List(pixels);

    // Check if there are NDVI values
    var ndviSize = ndviValues.size();

    // Proceed only if there are NDVI values
    return ee.Algorithms.If(ndviSize.gt(0), (function() {
      // Convert NDVI values to numeric list
      ndviValues = ndviValues.map(function(value) {
        return ee.Number(value);
      });

      // Sort NDVI values
      var ndviSorted = ndviValues.sort();

      // Compute Q1 and Q3 indices
      var q1Index = ndviSize.multiply(0.25).subtract(0.5).int();
      var q3Index = ndviSize.multiply(0.75).subtract(0.5).int();

      // Get Q1 and Q3 values
      var q1 = ee.Number(ndviSorted.get(q1Index));
      var q3 = ee.Number(ndviSorted.get(q3Index));

      // Compute IQR
      var iqr = q3.subtract(q1);

      // Define outlier bounds (1.5 times IQR)
      var lowerBound = q1.subtract(iqr.multiply(1.5));
      var upperBound = q3.add(iqr.multiply(1.5));

      // Filter out outliers
      var filteredNdviValues = ndviValues.filter(ee.Filter.gte('item', lowerBound))
                                         .filter(ee.Filter.lte('item', upperBound));

      // Compute statistics from filtered NDVI values
      var ndviMean = ee.Number(filteredNdviValues.reduce(ee.Reducer.mean()));
      var ndviMedian = ee.Number(filteredNdviValues.reduce(ee.Reducer.median()));
      var ndviMax = ee.Number(filteredNdviValues.reduce(ee.Reducer.max()));
      var ndviStdDev = ee.Number(filteredNdviValues.reduce(ee.Reducer.stdDev()));

      // Create a feature with station ID, date, and statistics
      return ee.Feature(null, {
        'station_id': stationId,
        'date': date,
        'NDVI_mean': ndviMean,
        'NDVI_median': ndviMedian,
        'NDVI_max': ndviMax,
        'NDVI_stdDev': ndviStdDev
      });
    })(), null);  // Return null if no NDVI values
  }, true);  // Set dropNulls to true in the map function

  // Return the FeatureCollection of statistics
  return statsByStation;
};

// ===============================================
// Step 6: Apply NDVI extraction to each image in the collection
// ===============================================
print('Step 6: Apply NDVI extraction to each image in the collection');

var ndviTimeSeries = landsatCollection.map(extractNDVIForImage).flatten();
print('Total number of NDVI time series features:', ndviTimeSeries.size());

// ===============================================
// Step 7: Print sample results for verification
// ===============================================
print('Step 7: Print sample NDVI time series results');
print('NDVI Time Series for all sites (first 10):', ndviTimeSeries.limit(10));

// ===============================================
// Step 8: Export NDVI time series to Google Drive
// ===============================================
print('Step 8: Export NDVI time series to Google Drive');

Export.table.toDrive({
  collection: ndviTimeSeries,
  description: 'Landsat8_NDVI_IQR_2014_06_2014_08_Station_120_180',
  fileFormat: 'CSV',
  selectors: ['station_id', 'date', 'NDVI_mean', 'NDVI_median', 'NDVI_max', 'NDVI_stdDev']
});

print('Export task has been submitted. Please monitor progress in the "Tasks" tab.');

LUCC提取

// ===============================================
// Step 1: 加载光伏站点矢量数据
// ===============================================
var region = ee.FeatureCollection("projects/ee-jiag/assets/InnerMongolia_station_60");

// ===============================================
// Step 2: 加载 ESA WorldCover v200 数据集
// ===============================================
var worldCover = ee.Image("ESA/WorldCover/v200/2021").select('Map');

// ===============================================
// Step 3: 使用多数决法提取每个站点的主要土地利用类型
// ===============================================
var extractMajorityLandCoverForSites = function(feature) {
  // 使用reduceRegion而不是reduceRegions, 只对每个站点计算一次
  var majorityLandCover = worldCover.reduceRegion({
    reducer: ee.Reducer.mode(), // 使用mode选择出现最多的土地利用类型
    geometry: feature.geometry(), // 计算每个站点的几何区域
    scale: 10,  // ESA WorldCover 数据的分辨率是10米
    maxPixels: 1e6 // 最大像元数设置
  });
  
  // 返回每个站点的主导土地利用类型
  return feature.set('land_cover', majorityLandCover.get('Map'));
};

// ===============================================
// Step 4: 应用提取土地利用类型的函数到每个站点
// ===============================================
var landCoverBySite = region.map(extractMajorityLandCoverForSites);
print('土地利用类型结果:', landCoverBySite.limit(10));

// ===============================================
// Step 5: 导出结果到 Google Drive
// ===============================================
Export.table.toDrive({
  collection: landCoverBySite,
  description: 'ESA_WorldCover_Majority_LandUse_Station60',
  fileFormat: 'CSV',
  selectors: ['编号', 'land_cover']  // 导出站点编号和多数决土地利用类型
});