GEE 9：Earth Engine Reducers 的基本操作

目錄

• • 1.Image 、ImageCollection and Regions Reducers（圖層和區域的相關操作）
  • • 1.1 Image Reductions（處理單個圖層）
    • 1.2 ImageCollection Reductions（處理圖層集）
    • 1.3 Greenest pixel (maximum NDVI) composite（了解去雲操作）
    • 1.4 ImageCollection reductions（注意事項）
    • 1.5 Image Region reductions
    • 1.6 Image Region reductions, cont.
    • 1.7 Image RegionS reductions
  • 2. Image Neighborhood reductions (convolution)（卷積Kernel）
  • 3.FeatureCollection column reductions（屬性運算）

1.Image 、ImageCollection and Regions Reducers（圖層和區域的相關操作）

Reducers are the way to aggregate data over time, space, bands, arrays and other data structures in Earth Engine.

Reducers 是在 Earth Engine 中聚合時間、空間、波段、數組和其他資料結構的資料的方式。 ee.Reducer 類指定資料的聚合方式。此類中的Reducer可以指定用于聚合的簡單統計資料（最小值、最大值、平均值、中值、标準差等），或輸入資料的更複雜的組合（例如直方圖、線性回歸、清單）。時間處理器：imageCollection.reduce()；空間 處理器：image.reduceRegion()、image.reduceNeighborhood()；波段處理器：image.reduce()； FeatureCollection 的屬性空間 ：featureCollection.reduceColumns()、以 aggregate_ 開頭的 FeatureCollection 方法。

顧名思義，Reducer翻譯為減少器，即将多個波段圖像減少到一個波段圖像的方法。

Reducers take an input dataset and produce a single output.

ee.Reducer()指定資料的聚合方式。

1.1 Image Reductions（處理單個圖層）

使用 image.reduce() 對一張圖像進行簡化，并且會對該圖像的所有波段進行處理，輸出圖像隻有一個波段，比如ee.Reducer.max()，但是ee.Reducer.minMax()有兩個波段輸出。

// Reducers: Image.reduce()
var image = l8sr.first()
var max = image.reduce(ee.Reducer.max());

Map.centerObject(image, 11);
Map.addLayer(max, {min: 1000, max: 50000}, 'max');
           

結果顯示：

1.2 ImageCollection Reductions（處理圖層集）

由ImageCollection表示的圖像時間序列取中值，輸出中的每個像素由每個波段在該位置的集合中所有像素的中值組成。按照逐像素和逐波段進行計算（pixel-wise、band-wise）。

var median = imageCollection.reduce(ee.Reducer.median(), 2); #注意這裡的2，parallelScale參數
           

注意：2 這裡設定的是parallelScale（用來限制記憶體使用的比例因子），使用更大的 parallelScale（例如 2 或 4）可能會啟用預設情況下耗盡記憶體的計算。

// Reducers: ImageCollection.reduce()
var median = l8sr.reduce(ee.Reducer.median(), 2);

var visParams = {
  bands: ['SR_B4_median', 'SR_B3_median', 'SR_B2_median'], 
  min: 7000, 
  max: 20000
}

Map.setCenter(116, 30, 8);
Map.addLayer(median, visParams, 'median');
           

結果顯示：

1.3 Greenest pixel (maximum NDVI) composite（了解去雲操作）

使用Landsat8圖層，計算NDVI，計算方法為(B5-B4)/(B5+B4)，并結合B2, B3, B4。max(4) reducer 有 4 個輸入，4 個輸出：保持元組具有第一個 (NDVI) 輸入的最大值，将 4 個波段圖像的堆疊轉換為單個 4 波段圖像。

計算流程：

• 加載Landsat8圖像，并使用filterDate()選擇時間，時間為2020年
• 使用 maskL8sr() 函數對圖層進行處理：去雲、去除飽和像素；适當縮放波段；計算并插入一個NDVI波段
• 選擇 ImageCollection 的四個波段（NDVI、B2、B3、B4），并使用 ee.Reducer.max(4) 将它們縮小為單個圖像。

// Computes a greenest-pixel composite of Landsat 8 SR images.
function maskL8sr(image) {
  var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);
  var saturationMask = image.select('QA_RADSAT').eq(0);

  // Apply the scaling factors to the appropriate bands.
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);

  // Replace the original bands with the scaled ones and apply the masks.
  return image.addBands(opticalBands, null, true)
      .addBands(thermalBands, null, true)
      .addBands(image.normalizedDifference(['SR_B5', 'SR_B4']))
      .updateMask(qaMask)
      .updateMask(saturationMask);
}

// The bands we'll have in our composite, with the NDVI band first.
var bands = ['nd', 'SR_B2', 'SR_B3', 'SR_B4'];

var imgs = l8sr.filterDate('2020-01-01', '2020-12-31');

// Use setOutputs() to keep the band names (otherwise they're
// placed with "max", "max2", "max3", etc.).
var greenestPixel =
      imgs.map(maskL8sr)
          .select(bands)
          .reduce(ee.Reducer.max(4).setOutputs(bands));

var visParams = {
  bands: ['SR_B4', 'SR_B3', 'SR_B2'],
  min: 0,
  max: 0.5,
  gamma: 1.4,
};

Map.addLayer(greenestPixel, visParams);

// Optionallly display a simple median composite for comparison.
Map.addLayer(imgs.map(maskL8sr).median(), visParams, 'median');

           

文法注釋：

• parseInt(‘11111’, 2)：強制轉換，将’11111’這個2進制數字轉化成了十進制數字31。
• bitwiseAnd() 即“按位與”：對應的兩個二進位都為1時，結果位才為1，如十進制位的6（0110）和10（1010），結果為0010。
• .eq(0) 是将影像中等于0的像素變為1，其它像素變為0。

去雲操作：

• 第一步，選取影像的 ‘QA_PIXEL’ 波段；
• 第二步，對其中的每一個像素和 ‘0000 0000 0001 1111’ 進行按位與運算（不到16位就往左補0）；
• 第三步，将值等于0的像素變為1，其它像素變為0。我們不難看出，第二步中，隻有’XXXX XXXX XXX0 0000’（X代表0或1）這樣的像素，在進行位運算後的結果為0，最終的結果為1。結合第一步中展示的表格，我們最終得到的影像中，值為0的像素就是函數中注釋提到的需要篩去的像素。

結果顯示：

median：

maxNDVI

1.4 ImageCollection reductions（注意事項）

這些 reducer 在 ImageCollection.reduce() 中不起作用，因為它們沒有數字輸出：

ee.Reducer.frequencyHistogram()

ee.Reducer.histogram()

ee.Reducer.toList()

如果不首先将每個輸入轉換為一維數組，這些縮減器将無法在 ImageCollection.reduce() 中工作：

ee.Reducer.covariance()

ee.Reducer.centeredCovariance()

1.5 Image Region reductions

var mean = image.reduceRegion({
  reducer: ee.Reducer.mean(), 
  geometry: santaCruz.geometry(),
  scale: 30,
  maxPixels: 1e9,
});
           

解釋：

reducer: 需要使用的 reducer

geometry: 資料範圍

scale: 每個像素的分辨率

crs: 投影，如果沒有則設定為圖像第一個波段的投影

maxPixels: 圖像最大像素數量。預設為10,000,000個. 這是防止意外運作大量計算的保護機制。

// Reducers: Image.reduceRegion()
var santaCruz = counties.filter(ee.Filter.eq('NAME', 'Santa Cruz')).first();

var image = l8sr.filterBounds(santaCruz.geometry()).first()

var mean = image.reduceRegion({
  reducer: ee.Reducer.mean(), 
  geometry: santaCruz.geometry(),
  scale: 30,
  maxPixels: 1e9,
});

print(mean)
           

The result is a set of <band, band_mean> tuples.

1.6 Image Region reductions, cont.

var meanDictionary = composite.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: region.geometry(),
  scale: 30,
  maxPixels: 1e10,
  tileScale: 4
});
           

此示例使用更大的 maxPixels 值以及使用 tileScale 參數。tileScale 是一個介于 0.1 和 16（預設為 1）之間的比例因子，用于調整聚合切片大小。設定較大的 tileScale（例如 2 或 4）會使用較小的圖塊，并且可能會啟用預設情況下耗盡記憶體的計算。（類似于并行計算，将需要計算的圖層分為2^n個）

提示：

• 計算被分成每個像素的子計算，每個子計算都可以由不同的伺服器完成。
• 增加 tileScale 會使每台伺服器的計算量更小，進而減少每台伺服器記憶體不足的可能性。
• 注意： 将計算拆分為更多像素的設定時間可能會使整體計算花費更長的時間。

// See:
// https://code.earthengine.google.com/?scriptPath=Examples%3ACloud%20Masking%2FLandsat8%20Surface%20Reflectance
function maskL8sr(image) {
  // Bit 0 - Fill
  // Bit 1 - Dilated Cloud
  // Bit 2 - Cirrus
  // Bit 3 - Cloud
  // Bit 4 - Cloud Shadow
  var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);
  var saturationMask = image.select('QA_RADSAT').eq(0);

  // Apply the scaling factors to the appropriate bands.
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);

  // Replace the original bands with the scaled ones and apply the masks.
  return image.addBands(opticalBands, null, true)
      .addBands(thermalBands, null, true)
      .updateMask(qaMask)
      .updateMask(saturationMask);
}

var collection = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
                     .filterDate('2020-01-01', '2021-01-01')
                     .map(maskL8sr);

var composite = collection.median();


// Load an input region: Sierra Nevada.
var region = ee.Feature(ee.FeatureCollection('EPA/Ecoregions/2013/L3')
  .filter(ee.Filter.eq('us_l3name', 'Sierra Nevada'))
  .first());

// Reduce the region. The region parameter is the Feature geometry.
var meanDictionary = composite.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: region.geometry(),
  scale: 30,
  maxPixels: 1e10,
  tileScale: 4
});

// The result is a Dictionary.  Print it.
print(meanDictionary);
           

1.7 Image RegionS reductions

var featuresWithMeans = image.reduceRegions({
  collection: maineCounties,
  reducer: ee.Reducer.mean(),
  scale: 30,
  tileScale: 4
});
           

// Load a FeatureCollection of counties in Maine.
var maineCounties = ee.FeatureCollection('TIGER/2016/Counties')
  .filter(ee.Filter.eq('STATEFP', '23'));

var image = l8sr
    .filterBounds(maineCounties)
    .filterDate('2020-01-01', '2020-12-31')
    .median()

// Add reducer output to the Features in the collection.
var maineMeansFeatures = image.reduceRegions({
  collection: maineCounties,
  reducer: ee.Reducer.mean(),
  scale: 30,
  tileScale: 2
});

// Print the first feature, to illustrate the result.
print(ee.Feature(maineMeansFeatures.first()).select(image.bandNames()));
           

使用 reduceRegions() 執行圖像區域縮減。傳回輸入特征，每個特征都增加了相應的 reducer 輸出。

2. Image Neighborhood reductions (convolution)（卷積Kernel）

// Reducers: Image.reduceNeighborhood()

// Region in the redwood forest
var redwoods = ee.Geometry.Rectangle(-124.0665, 41.0739, -123.934, 41.2029);

// Input NAIP imagery
//National Agriculture Imagery Program (NAIP) ImageCollection
var naipCollection = ee.ImageCollection('USDA/NAIP/DOQQ')
  .filterBounds(redwoods)
  .filterDate('2012-01-01', '2012-12-31');
var naip = naipCollection.mosaic();

// Compute NDVI from the NAIP
var naipNDVI = naip.normalizedDifference(['N', 'R']);

// Standard Deviation (SD) as texture of NDVI
// a circle kernel of radius 7 pixels
var texture = naipNDVI.reduceNeighborhood({
  reducer: ee.Reducer.stdDev(), 
  kernel: ee.Kernel.circle(7), 
});

// Display
Map.centerObject(redwoods, 12);
Map.addLayer(naip, {}, 'NAIP input imagery');
Map.addLayer(naipNDVI, {min: -1, max: 1, palette:['FF0000', '00FF00']}, 'NDVI');
Map.addLayer(texture, {min:0, max: 0.3}, 'SD of NDVI');
           

NAIP:

NDVI:

SD：表示NDVI的紋理

3.FeatureCollection column reductions（屬性運算）

腳本流程：

1.加載美國人口普查資料，并在Benton County選擇需要的屬性

2.選擇人口和房屋單元兩個屬性

3.實作reduceColumn操作：selectors設定為人口和房屋單元屬性；reducer設定為ee.Reducer.sum().repeat(2)。

Notes: Unlike imageCollection.reduce(), in which reducers are automatically repeated for each band, reducers on a FeatureCollection must be explicitly repeated using repeat(). Specifically, repeat the reducer m times for m inputs.

// Consider a FeatureCollection of US census blocks with census data as attributes.
// The variables of interest are total population and total housing units.
// You can get their sum(s) by supplying a summing reducer argument to 
// reduceColumns() and printing the result.

// Load US cenus data as a FeatureCollection.
var census = ee.FeatureCollection('TIGER/2010/Blocks');

// Filter the collection to include only Benton County, OR.
var benton = census.filter(
  ee.Filter.and(
    ee.Filter.eq('statefp10', '41'),
    ee.Filter.eq('countyfp10', '003')
  )
);

// Display Benton County cenus blocks.
Map.setCenter(-123.27, 44.57, 13);
Map.addLayer(benton);

// Compute sums of the specified properties.
var properties = ['pop10', 'housing10'];
var sums = benton
    .filter(ee.Filter.notNull(properties))
    .reduceColumns({
      reducer: ee.Reducer.sum().repeat(2),
      selectors: properties
    });

// Print the resultant Dictionary.
print(sums);
           

結果顯示：