Google earth engine 中的投影、重采样、尺度

本文主要翻译自下述GEE官方帮助

https://developers.google.com/earth-engine/guides/scale
https://developers.google.com/earth-engine/guides/projections
https://developers.google.com/earth-engine/guides/resample#resampling

1 尺度(scale)

尺度(scale)指像元分辨率。GEE中进行分析计算时以影像金字塔pyramid (瓦片)形式组织，尺度由输出(output)决定，而不是输入(input)。

pyramid 中给定level上的每个像素都是由下一level上2x2像素块聚合计算出来的。对于连续值图像(continuous valued images)，pyramid 上层的像素值是下一层像素的平均值。对于离散值图像(discrete valued images)，金字塔上层的像素值是下一层像素的采样sample(通常是左上角(top left pixel)的像素)。

pyramid 的最底层代表了图像数据的原始分辨率，之间被不断聚合，直到达到最顶层，一个256x256像素的tile。
当在代码中使用图像数据时，Earth Engine将根据指定的尺度选择最接近的(小于或等于)pyramid level作为输入，必要的话进行重采样(默认使用nearest neighbor)

Earth Engine的分析尺度取决于 a "pull" basis(地图缩放水平)。请求用于计算的输入尺度由输出决定。

var image = ee.Image('LANDSAT/LC08/C01/T1/LC08_044034_20140318').select('B4');

var printAtScale = function(scale) {
  print('Pixel value at '+scale+' meters scale',
    image.reduceRegion({
      reducer: ee.Reducer.first(),
      geometry: image.geometry().centroid(),
      // The scale determines the pyramid level from which to pull the input
      scale: scale
  }).get('B4'));
};

printAtScale(10); // 8883
printAtScale(30); // 8883
printAtScale(50); // 8337
printAtScale(70); // 9215
printAtScale(200); // 8775
printAtScale(500); // 8300

2 投影

GEE中的投影(projection)

Earth Engine的设计使我们在进行计算(computations)时很少需要担心地图投影(map projections)。和尺度(scale)一样，投影(projection)发生于哪一种计算取决于 a "pull" basis(地图缩放水平)。确切地讲，输入(input)投影由输出(output)投影决定。输出的投影又取决于：

1. 函数参数设置(如 crs)

2. Code Editor中的 Map (投影为 maps mercator (EPSG:3857) )，当在Code Editor中展示影像时，输入被请求为Mercator 投影

3. 函数 reproject() 调用

// The input image has a SR-ORG:6974 (sinusoidal) projection.
var image = ee.Image('MODIS/006/MOD13A1/2014_05_09').select(0);

// Normalize the image and add it to the map.
var rescaled = image.unitScale(-2000, 10000);
var visParams = {min: 0.15, max: 0.7};
Map.addLayer(rescaled, visParams, 'Rescaled');

在GEE中运行上述代码时。如下图所示，输入MOD13A1为sinusoidal投影，对数据进行标准化unitScale()后，利用Map在Code Editor中展示。因为Map展示投影为Mercator，所以unitScale()运算在Mercator投影下进行

Earth Engine中，投影由Coordinate Reference System指定 (CRS or the crs 很多方法中的参数)。调用 projection() 查看影像投影。

var image = ee.Image('LANDSAT/LC08/C01/T1/LC08_044034_20140318').select(0);//针对波段，有些影像不同波段具有不同投影
print('Projection, crs, and crs_transform:', image.projection());//影像投影信息
print('Scale in meters:', image.projection().nominalScale());//影像原始投影

GEE中的默认投影(default projection)

只有当计算需要指定投影时才重投影，其他情况下没必要。只有当输出投影不明确时，Earth Engine才会要求指定投影和/或尺度，不明确通常由reduce包含不同投影影像的影像集合(ImageCollection) 导致。

包含不同投影影像的影像集合composite或mosaic生成的影像默认投影是：WGS84(EPSG:4326), 1-degree scale 。因此使用这类影像进行计算时将会报错The default WGS84 projection is invalid for aggregations. Specify a scale or crs & crs_transform. 因为通常并不希望以 1-degree scale 进行aggregation ，故Earth Engine给出上述提醒为输出指定完整的投影定义。

GEE中的重投影(reprojection)

reproject()方法使调用reproject()方法之前的代码中的计算在该方法指定的投影中进行，重采样方法默认为最邻近法。一些需要固定投影的情况如下：

1. Computing gradients (e.g. ee.Terrain.gradient or ee.Terrain.slope)
2. reduceResolution ，当将高分辨率像元聚合为(aggregate )低空间分辨率时

然而需要尽量避免调用reproject()方法。

例如，当重新投影影像并将其添加到map中。如果 reproject()调用中指定的尺度比地图的缩放级别小得多，Earth Engine将请求非常小的像元尺度、非常大的空间范围内的所有输入。这可能会导致一次请求太多数据并导致错误。

// The input image has a SR-ORG:6974 (sinusoidal) projection.
var image = ee.Image('MODIS/006/MOD13A1/2014_05_09').select(0);

// Operations *before* the reproject call will be done in the projection
// specified by reproject().  The output results in another reprojection.
var reprojected = image
    .unitScale(-2000, 10000)
    .reproject('EPSG:4326', null, 500);
Map.addLayer(reprojected, visParams, 'Reprojected');

上述代码，输入MOD13A1为sinusoidal投影，对数据进行标准化 unitScale() 后调用reproject()方法，所以unitScale()在WGS84投影下进行。最后利用Map在Code Editor中展示时因为Map展示投影为Mercator，又需要将结果重投影为Mercator。

3 重采样(resample)

Earth Engine 重投影默认执行最近邻重采样。可以使用resample() 或 reduceResolution() 更改。具体地说，当这些方法应用于输入图像时，输入图像所需的任何重投影都将使用指定的重采样或聚合方法完成。

resample()

resample()使用指定的重采样方法('bilinear' or 'bicubic') 用于下一次重投影。由于Earth Engine根据输出投影请求输入的投影，因此隐式的( implicit)重投影可能在对输入图像进行任何其它操作之前发生。因此，直接在输入图像上调用resample()。

// Load a Landsat image over San Francisco, California, UAS.
var landsat = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_044034_20160323');

// Set display and visualization parameters.
Map.setCenter(-122.37383, 37.6193, 15);
var visParams = {bands: ['B4', 'B3', 'B2'], max: 0.3};

// Display the Landsat image using the default nearest neighbor resampling.
// when reprojecting to Mercator for the Code Editor map.
Map.addLayer(landsat, visParams, 'original image');

// Force the next reprojection on this image to use bicubic resampling.
var resampled = landsat.resample('bicubic');

// Display the Landsat image using bicubic resampling.
Map.addLayer(resampled, visParams, 'resampled');

reduceResolution()

当我们的目标不是在重投影时进行重采样，而是将小尺度像元采用不同的投影聚合为一个更大尺度的像元，需要使用此功能。这一功能适用于对比具有不同尺度的数据集，例如对比30米的Landsat数据与500米的MODIS数据。具体使用方法看下面代码。

// Load a MODIS EVI image.
var modis = ee.Image(ee.ImageCollection('MODIS/006/MOD13A1').first())
    .select('EVI');

// Display the EVI image near La Honda, California.
Map.setCenter(-122.3616, 37.5331, 12);
// 这里在zoom level为12的pyramids水平下，默认最邻近采样，显示MODIS数据
// MODIS数据原始投影为sinusoidal投影，加载到地图上时转化为Mercator
Map.addLayer(modis, {min: 2000, max: 5000}, 'MODIS EVI');//projection: SR-ORG:6974

// Get information about the MODIS projection.
var modisProjection = modis.projection();
print('MODIS projection:', modisProjection);//projection: SR-ORG:6974

// Load and display forest cover data at 30 meters resolution.
var forest = ee.Image('UMD/hansen/global_forest_change_2015')
    .select('treecover2000');
Map.addLayer(forest, {max: 80}, 'forest cover 30 m');//projection: EPSG:4326

// Get the forest cover data at MODIS scale and projection.
var forestMean = forest
    // Force the next reprojection to aggregate instead of resampling.
    .reduceResolution({
      reducer: ee.Reducer.mean(),
      maxPixels: 1024
    })
    // Request the data at the scale and projection of the MODIS image.
    .reproject({
      crs: modisProjection
    });

// Display the aggregated, reprojected forest cover data.
Map.addLayer(forestMean, {max: 80}, 'forest cover at MODIS scale');//projection: SR-ORG:6974

聚合像元值，依据面积加权思想合成。如下图所示，输出像元的面积为a，与b相交的输入像元的权重计算为b/a，与c相交的输入像元的权重计算为c/a。

在进行一些应用时，使用除mean reducer之外的reducer时需要注意上述合成规则避免出现错误。例如

例如，要计算每个像元的森林面积，使用mean reducer来计算所覆盖像素的百分比，然后乘以面积(而不是计算较小像素的面积，然后用sum reducer来将它们相加，因为这里使用sum的话会利用上述合成规则有一个加权过程导致不准)。

// Compute forest area per MODIS pixel.
var forestArea = forest.gt(0)
    // Force the next reprojection to aggregate instead of resampling.
    .reduceResolution({
      reducer: ee.Reducer.mean(),
      maxPixels: 1024
    })
    // The reduce resolution returns the fraction of the MODIS pixel
    // that's covered by 30 meter forest pixels.  Convert to area
    // after the reduceResolution() call.
    .multiply(ee.Image.pixelArea())
    // Request the data at the scale and projection of the MODIS image.
    .reproject({
      crs: modisProjection
    });
Map.addLayer(forestArea, {max: 500 * 500}, 'forested area at MODIS scale');

欢迎关注本人公众号，获取更多有深度、有内容、有趣味的信息。
公众号：“囚室”