数组转换

Earth Engine 支持数组转换,例如转置、转置和伪转置。 例如,假设对一组图片的时间序列进行普通最小二乘 (OLS) 回归。在以下示例中,将包含预测器和响应波段的图片转换为数组图片,然后通过三种方式“求解”以获取最小二乘系数估计值。首先,组装图片数据并转换为数组:

Code Editor (JavaScript)

// Scales and masks Landsat 8 surface reflectance images.
function prepSrL8(image) {
  // Develop masks for unwanted pixels (fill, cloud, 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);
}

// Load a Landsat 8 surface reflectance image collection.
var collection = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
  // Filter to get only two years of data.
  .filterDate('2019-04-01', '2021-04-01')
  // Filter to get only imagery at a point of interest.
  .filterBounds(ee.Geometry.Point(-122.08709, 36.9732))
  // Prepare images by mapping the prepSrL8 function over the collection.
  .map(prepSrL8)
  // Select NIR and red bands only.
  .select(['SR_B5', 'SR_B4'])
  // Sort the collection in chronological order.
  .sort('system:time_start', true);

// This function computes the predictors and the response from the input.
var makeVariables = function(image) {
  // Compute time of the image in fractional years relative to the Epoch.
  var year = ee.Image(image.date().difference(ee.Date('1970-01-01'), 'year'));
  // Compute the season in radians, one cycle per year.
  var season = year.multiply(2 * Math.PI);
  // Return an image of the predictors followed by the response.
  return image.select()
    .addBands(ee.Image(1))                                  // 0. constant
    .addBands(year.rename('t'))                             // 1. linear trend
    .addBands(season.sin().rename('sin'))                   // 2. seasonal
    .addBands(season.cos().rename('cos'))                   // 3. seasonal
    .addBands(image.normalizedDifference().rename('NDVI'))  // 4. response
    .toFloat();
};

// Define the axes of variation in the collection array.
var imageAxis = 0;
var bandAxis = 1;

// Convert the collection to an array.
var array = collection.map(makeVariables).toArray();

// Check the length of the image axis (number of images).
var arrayLength = array.arrayLength(imageAxis);
// Update the mask to ensure that the number of images is greater than or
// equal to the number of predictors (the linear model is solvable).
array = array.updateMask(arrayLength.gt(4));

// Get slices of the array according to positions along the band axis.
var predictors = array.arraySlice(bandAxis, 0, 4);
var response = array.arraySlice(bandAxis, 4);

Python 设置

如需了解 Python API 以及如何使用 geemap 进行交互式开发,请参阅 Python 环境页面。

import ee
import geemap.core as geemap

Colab (Python)

import math


# Scales and masks Landsat 8 surface reflectance images.
def prep_sr_l8(image):
  # Develop masks for unwanted pixels (fill, cloud, cloud shadow).
  qa_mask = image.select('QA_PIXEL').bitwiseAnd(int('11111', 2)).eq(0)
  saturation_mask = image.select('QA_RADSAT').eq(0)

  # Apply the scaling factors to the appropriate bands.
  optical_bands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
  thermal_bands = 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(optical_bands, None, True)
      .addBands(thermal_bands, None, True)
      .updateMask(qa_mask)
      .updateMask(saturation_mask)
  )


# Load a Landsat 8 surface reflectance image collection.
collection = (
    ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    # Filter to get only two years of data.
    .filterDate('2019-04-01', '2021-04-01')
    # Filter to get only imagery at a point of interest.
    .filterBounds(ee.Geometry.Point(-122.08709, 36.9732))
    # Prepare images by mapping the prep_sr_l8 function over the collection.
    .map(prep_sr_l8)
    # Select NIR and red bands only.
    .select(['SR_B5', 'SR_B4'])
    # Sort the collection in chronological order.
    .sort('system:time_start', True)
)


# This function computes the predictors and the response from the input.
def make_variables(image):
  # Compute time of the image in fractional years relative to the Epoch.
  year = ee.Image(image.date().difference(ee.Date('1970-01-01'), 'year'))
  # Compute the season in radians, one cycle per year.
  season = year.multiply(2 * math.pi)
  # Return an image of the predictors followed by the response.
  return (
      image.select()
      .addBands(ee.Image(1))  # 0. constant
      .addBands(year.rename('t'))  # 1. linear trend
      .addBands(season.sin().rename('sin'))  # 2. seasonal
      .addBands(season.cos().rename('cos'))  # 3. seasonal
      .addBands(image.normalizedDifference().rename('NDVI'))  # 4. response
      .toFloat()
  )


# Define the axes of variation in the collection array.
image_axis = 0
band_axis = 1

# Convert the collection to an array.
array = collection.map(make_variables).toArray()

# Check the length of the image axis (number of images).
array_length = array.arrayLength(image_axis)
# Update the mask to ensure that the number of images is greater than or
# equal to the number of predictors (the linear model is solvable).
array = array.updateMask(array_length.gt(4))

# Get slices of the array according to positions along the band axis.
predictors = array.arraySlice(band_axis, 0, 4)
response = array.arraySlice(band_axis, 4)

请注意,arraySlice() 会返回时间序列中对于沿 bandAxis(1 轴)指定的索引范围的所有图片。此时,可以使用矩阵代数来求解 OLS 系数:

Code Editor (JavaScript)

// Compute coefficients the hard way.
var coefficients1 = predictors.arrayTranspose().matrixMultiply(predictors)
  .matrixInverse().matrixMultiply(predictors.arrayTranspose())
    .matrixMultiply(response);

Python 设置

如需了解 Python API 以及如何使用 geemap 进行交互式开发,请参阅 Python 环境页面。

import ee
import geemap.core as geemap

Colab (Python)

# Compute coefficients the hard way.
coefficients_1 = (
    predictors.arrayTranspose()
    .matrixMultiply(predictors)
    .matrixInverse()
    .matrixMultiply(predictors.arrayTranspose())
    .matrixMultiply(response)
)

虽然这种方法可行,但效率不高,并且会导致代码难以阅读。更好的方法是使用 pseudoInverse() 方法(对于数组图片,请使用 matrixPseudoInverse()):

Code Editor (JavaScript)

// Compute coefficients the easy way.
var coefficients2 = predictors.matrixPseudoInverse()
  .matrixMultiply(response);

Python 设置

如需了解 Python API 以及如何使用 geemap 进行交互式开发,请参阅 Python 环境页面。

import ee
import geemap.core as geemap

Colab (Python)

# Compute coefficients the easy way.
coefficients_2 = predictors.matrixPseudoInverse().matrixMultiply(response)

从可读性和计算效率的角度来看,获取 OLS 系数的最佳方法是 solve()(对于数组图片,则为 matrixSolve())。solve() 函数会根据输入的特征确定如何最佳地求解系统,对于过秩系统,使用伪逆;对于方阵,使用逆;对于接近奇异的矩阵,使用特殊技术:

Code Editor (JavaScript)

// Compute coefficients the easiest way.
var coefficients3 = predictors.matrixSolve(response);

Python 设置

如需了解 Python API 以及如何使用 geemap 进行交互式开发,请参阅 Python 环境页面。

import ee
import geemap.core as geemap

Colab (Python)

# Compute coefficients the easiest way.
coefficients_3 = predictors.matrixSolve(response)

如需获取多波段图像,请将数组图像投影到更低维的空间,然后将其展平:

Code Editor (JavaScript)

// Turn the results into a multi-band image.
var coefficientsImage = coefficients3
  // Get rid of the extra dimensions.
  .arrayProject([0])
  .arrayFlatten([
    ['constant', 'trend', 'sin', 'cos']
]);

Python 设置

如需了解 Python API 以及如何使用 geemap 进行交互式开发,请参阅 Python 环境页面。

import ee
import geemap.core as geemap

Colab (Python)

# Turn the results into a multi-band image.
coefficients_image = (
    coefficients_3
    # Get rid of the extra dimensions.
    .arrayProject([0]).arrayFlatten([['constant', 'trend', 'sin', 'cos']])
)

检查这三种方法的输出,您会发现无论使用哪种求解器,系数的最终矩阵都是相同的。solve() 灵活且高效,因此非常适合通用线性建模。