利用Google Earth Engine索引优化海量点SMAP L4时间序列采样

针对超30万个聚集点位的SMAP L4时间序列提取需求，核心思路是先给每个点位标记所属SMAP像素的唯一标识，再对去重后的像素提取时间序列，最后关联回原始点位，避免重复计算同一像素。以下是具体实现步骤：

1. 获取SMAP投影与分辨率 #

首先获取SMAP L4数据的投影参数，确保后续点位坐标转换匹配像素网格：

// 加载SMAP L4影像集并提取投影信息
var smapCollection = ee.ImageCollection('NASA_USDA/HSL/SMAP_soil_moisture');
var smapFirst = smapCollection.first();
var smapProjection = smapFirst.projection();
var smapScale = smapProjection.nominalScale(); // SMAP像素分辨率为3600米

2. 给采样点添加像素唯一标识 #

通过坐标转换，为每个采样点计算所属SMAP像素的行列号或像素中心经纬度，作为唯一标识：

方法1：用行列号生成像素ID #

var pointsWithPixelId = yourPointCollection.map(function(point) {
  // 将WGS84点位转换到SMAP投影坐标系
  var projectedPoint = point.transform(smapProjection, 1);
  var projectedCoords = projectedPoint.coordinates();
  
  // 计算该点所在像素的行列号（投影坐标除以分辨率后取整）
  var pixelX = ee.Number(projectedCoords.get(0)).divide(smapScale).floor();
  var pixelY = ee.Number(projectedCoords.get(1)).divide(smapScale).floor();
  
  // 生成唯一像素ID（行列号拼接）
  var pixelId = pixelX.format('%d').cat('_').cat(pixelY.format('%d'));
  
  return point.set({
    pixel_x: pixelX,
    pixel_y: pixelY,
    pixel_id: pixelId
  });
});

方法2：用像素中心经纬度作为标识 #

如果需要更直观的地理坐标作为关联依据，可计算像素中心的WGS84坐标：

var pointsWithPixelCenter = yourPointCollection.map(function(point) {
  var projectedPoint = point.transform(smapProjection, 1);
  var projectedCoords = projectedPoint.coordinates();
  
  // 计算像素左上角坐标，再偏移半分辨率得到中心坐标
  var pixelLeftX = ee.Number(projectedCoords.get(0)).divide(smapScale).floor().multiply(smapScale);
  var pixelTopY = ee.Number(projectedCoords.get(1)).divide(smapScale).floor().multiply(smapScale);
  var centerX = pixelLeftX.add(smapScale.divide(2));
  var centerY = pixelTopY.add(smapScale.divide(2));
  
  // 转换回WGS84坐标系
  var pixelCenter = ee.Geometry.Point([centerX, centerY], smapProjection).transform('EPSG:4326');
  
  return point.set({
    pixel_center_lon: pixelCenter.coordinates().get(0),
    pixel_center_lat: pixelCenter.coordinates().get(1)
  });
});

3. 提取唯一像素集合 #

基于像素ID去重，得到需要提取时间序列的唯一像素集合：

// 基于pixel_id去重，得到唯一像素点位
var uniquePixels = pointsWithPixelId.distinct('pixel_id');

4. 提取唯一像素的SMAP时间序列 #

对去重后的像素批量提取时间序列，避免重复计算：

var smapTimeSeries = uniquePixels.map(function(pixel) {
  var pixelId = pixel.get('pixel_id');
  // 对每个像素遍历SMAP影像集，提取对应日期的土壤湿度值
  var pixelTs = smapCollection.map(function(image) {
    var ssmValue = image.reduceRegion({
      reducer: ee.Reducer.first(),
      geometry: pixel.geometry(),
      scale: smapScale,
      crs: smapProjection
    }).get('ssm'); // 按需选择SMAP的波段，如ssm（表层土壤湿度）、susm（深层土壤湿度）
    
    return ee.Feature(null, {
      pixel_id: pixelId,
      date: image.date().format('YYYY-MM-dd'),
      ssm: ssmValue
    });
  });
  return pixelTs.flatten();
}).flatten();

5. 将时间序列关联回原始点位 #

通过Join操作，将像素时间序列与原始点位按像素ID关联：

// 构建匹配过滤器：按pixel_id关联
var joinFilter = ee.Filter.equals({
  leftField: 'pixel_id',
  rightField: 'pixel_id'
});

// 执行内连接，关联点位与时间序列
var joinedFeatures = ee.Join.inner('point_attrs', 'ts_attrs').apply(
  pointsWithPixelId, 
  smapTimeSeries, 
  joinFilter
);

// 整理最终结果：合并点位属性与时间序列属性
var finalResults = joinedFeatures.map(function(joinedFeature) {
  var pointProps = joinedFeature.get('point_attrs');
  var tsProps = joinedFeature.get('ts_attrs');
  return ee.Feature(null, ee.Dictionary(pointProps).combine(tsProps));
});

关于ee.Image.PixelCoordinates的说明 #

ee.Image.PixelCoordinates是生成包含像素行列号的影像，若直接用reduceRegions提取每个点位的行列号，本质还是对每个点位单独计算，无法避免重复处理同一像素的多个点位。而上述方法先批量计算点位的像素标识再去重，能大幅减少后续时间序列的计算量，更适合大数量聚集点位的场景。

内容的提问来源于stack exchange，提问作者paul_hutton_maps