# Array Transformations

Earth Engine supports array transformations such as transpose, inverse and pseudo-inverse.
As an example, consider an ordinary least squares (OLS) regression of a time series of
images. In the following example, an image with bands for predictors and a response is
converted to an array image, then "solved" to obtain least squares coefficients estimates
three ways. First, assemble the image data and convert to arrays:

### Code Editor (JavaScript)

```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 setup

See the [Python Environment](https://developers.google.com/earth-engine/guides/python_install) page for information on the Python API and using
`geemap` for interactive development.

```python
import ee
import geemap.core as geemap
```

### Colab (Python)

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

Note that `arraySlice()` returns all the images in the time series for the
range of indices specified along the `bandAxis` (the 1-axis). At this point,
matrix algebra can be used to solve for the OLS coefficients:

### Code Editor (JavaScript)

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

See the [Python Environment](https://developers.google.com/earth-engine/guides/python_install) page for information on the Python API and using
`geemap` for interactive development.

```python
import ee
import geemap.core as geemap
```

### Colab (Python)

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

Although this method works, it is inefficient and makes for difficult to read code. A
better way is to use the `pseudoInverse()` method
(`matrixPseudoInverse()` for an array image):

### Code Editor (JavaScript)

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

See the [Python Environment](https://developers.google.com/earth-engine/guides/python_install) page for information on the Python API and using
`geemap` for interactive development.

```python
import ee
import geemap.core as geemap
```

### Colab (Python)

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

From a readability and computational efficiency perspective, the best way to get the OLS
coefficients is `solve()` (`matrixSolve()` for an array image). The
`solve()` function determines how to best solve the system from characteristics
of the inputs, using the pseudo-inverse for overdetermined systems, the inverse for square
matrices and special techniques for nearly singular matrices:

### Code Editor (JavaScript)

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

See the [Python Environment](https://developers.google.com/earth-engine/guides/python_install) page for information on the Python API and using
`geemap` for interactive development.

```python
import ee
import geemap.core as geemap
```

### Colab (Python)

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

To get a multi-band image, project the array image into a lower dimensional space, then
flatten it:

### Code Editor (JavaScript)

```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 setup

See the [Python Environment](https://developers.google.com/earth-engine/guides/python_install) page for information on the Python API and using
`geemap` for interactive development.

```python
import ee
import geemap.core as geemap
```

### Colab (Python)

```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 the outputs of the three methods and observe that the resultant matrix of
coefficients is the same regardless of the solver. That `solve()` is flexible
and efficient makes it a good choice for general purpose linear modeling.