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');

歡迎關注本人公眾號，獲取更多有深度、有內容、有趣味的信息。
公眾號：“囚室”