Transformações de matriz

O Earth Engine oferece suporte a transformações de matrizes, como transposição, inversa e pseudo-inversa. Como exemplo, considere uma regressão de mínimos quadrados ordinários (OLS) de uma série temporal de imagens. No exemplo a seguir, uma imagem com bandas para preditores e uma resposta é convertida em uma imagem de matriz e, em seguida, "resolvida" para estimar os coeficientes de mínimos quadrados de três maneiras. Primeiro, monte os dados da imagem e converta em matrizes:

Editor de código (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);

Configuração do Python

Consulte a página Ambiente Python para informações sobre a API Python e o uso de geemap para desenvolvimento interativo.

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() retorna todas as imagens na série temporal para o intervalo de índices especificado ao longo de bandAxis (eixo 1). Nesse ponto, a álgebra matricial pode ser usada para resolver os coeficientes de OLS:

Editor de código (JavaScript)

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

Configuração do Python

Consulte a página Ambiente Python para informações sobre a API Python e o uso de geemap para desenvolvimento interativo.

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

Embora esse método funcione, ele é ineficiente e dificulta a leitura do código. Uma maneira melhor é usar o método pseudoInverse() (matrixPseudoInverse() para uma imagem de matriz):

Editor de código (JavaScript)

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

Configuração do Python

Consulte a página Ambiente Python para informações sobre a API Python e o uso de geemap para desenvolvimento interativo.

import ee
import geemap.core as geemap

Colab (Python)

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

Do ponto de vista da legibilidade e da eficiência computacional, a melhor maneira de obter os coeficientes de OLS é solve() (matrixSolve() para uma imagem de matriz). A função solve() determina a melhor forma de resolver o sistema com base nas características das entradas, usando o pseudoinverso para sistemas superdeterminados, o inverso para matrizes quadradas e técnicas especiais para matrizes quase singulares:

Editor de código (JavaScript)

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

Configuração do Python

Consulte a página Ambiente Python para informações sobre a API Python e o uso de geemap para desenvolvimento interativo.

import ee
import geemap.core as geemap

Colab (Python)

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

Para conseguir uma imagem de várias bandas, projete a imagem do array em um espaço dimensional menor e, em seguida, aplaine-a:

Editor de código (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']
]);

Configuração do Python

Consulte a página Ambiente Python para informações sobre a API Python e o uso de geemap para desenvolvimento interativo.

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

Examine as saídas dos três métodos e observe que a matriz resultante de coeficientes é a mesma, independente do solucionador. O fato de solve() ser flexível e eficiente faz com que ele seja uma boa escolha para modelagem linear de uso geral.