Search code examples
javascriptpcagoogle-earth-engine

Performing PCA per image over imageCollection in Google Earth Engine


I need to perform a PCA per image over a image collection. Then, I want to only keep Principle component axis 1, and add this as a band to every image within my image collection. Ultimately, I want to export a .csv file with GPS sampling locations at row headers and image ID as column headers with mean Principle component axis 1 as values. The idea behind doing this, is that I want a proxy (spectral heterogeneity) to use in further statistical analysis in R.

Here is the code I have thus far:

//Create an test image to extract information to be used during PCA
var testImage =ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_168080_20130407')
.select(['B2', 'B3', 'B4', 'B5', 'B6', 'B7'],
        ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2']);

// Define variables for PCA
var region = Extent;
var scale = testImage.projection().nominalScale();
var bandNames = testImage.bandNames();
Map.centerObject(region);

// Function for performing PCA
function doPCA(image){
  // This code is from https://code.earthengine.google.com/7249153a8a0f5c79eaf562ed45a7adad
var meanDict = image.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: region,
    scale: scale,
    maxPixels: 1e9
});
var means = ee.Image.constant(meanDict.values(bandNames));
var centered = image.subtract(means);

// This helper function returns a list of new band names.
var getNewBandNames = function(prefix) {
  var seq = ee.List.sequence(1, bandNames.length());
  return seq.map(function(b) {
    return ee.String(prefix).cat(ee.Number(b).int());
  });
};

// [START principal_components]
var getPrincipalComponents = function(centered, scale, region) {
  var arrays = centered.toArray();
  var covar = arrays.reduceRegion({
    reducer: ee.Reducer.centeredCovariance(),
    geometry: region,
    scale: scale,
    maxPixels: 1e9
  });
  var covarArray = ee.Array(covar.get('array'));
  var eigens = covarArray.eigen();
  var eigenValues = eigens.slice(1, 0, 1);
  var eigenVectors = eigens.slice(1, 1);
  var arrayImage = arrays.toArray(1);
  var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);
  var sdImage = ee.Image(eigenValues.sqrt())
    .arrayProject([0]).arrayFlatten([getNewBandNames('sd')]);
  return principalComponents
    .arrayProject([0])
    .arrayFlatten([getNewBandNames('pc')])
    .divide(sdImage);
};
var pcImage = getPrincipalComponents(centered, scale, region);
return (pcImage);
}

// map PCA function over collection
var PCA = LandsatCol.map(function(image){return doPCA(image)});
print('pca', PCA);

Extent is my ROI, whereas LandsatCol is a preproccessed image collection. The code here produces an Error when trying to map the PCA over the image collection (second last line of code). The error reads: "Array: Parameter 'values' is required".

Any suggestions on how to deal with this? And how to add Principle component axis 1 as a band per image over the image collection?


Solution

  • I figured it out. The error "Array: Parameter 'values' is required" had to do with sparse matrices, which was a product of filtering, clipping and spesifying regions within to perform PCA. Earth Engine can not work with sparse matrices.

    Here is the working code. LandsatCol is my preproccessed image collection.

    // Display AOI
    var point = ee.Geometry.Point([30.2261, -29.458])
    Map.centerObject(point,10);
    
    // Prepairing imagery for PCA
    var Preped = LandsatCol.map(function(image){
      var orig = image;
      var region = image.geometry();
      var scale = 30;
      var bandNames = ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'];
      var meanDict = image.reduceRegion({
        reducer: ee.Reducer.mean(),
        geometry: region,
        scale: scale,
        maxPixels: 1e9
      });
      var means = ee.Image.constant(meanDict.values(bandNames));
      var centered = image.subtract(means);
      var getNewBandNames = function(prefix) {
      var seq = ee.List.sequence(1, 6);
      return seq.map(function(b) {
        return ee.String(prefix).cat(ee.Number(b).int());
        });
      };
    
      // PCA function
      var getPrincipalComponents = function(centered, scale, region) {
        var arrays = centered.toArray();
        var covar = arrays.reduceRegion({
          reducer: ee.Reducer.centeredCovariance(),
          geometry: region,
          scale: scale,
          maxPixels: 1e9
        });
        var covarArray = ee.Array(covar.get('array'));
        var eigens = covarArray.eigen();
        var eigenValues = eigens.slice(1, 0, 1);
        var eigenVectors = eigens.slice(1, 1);
        var arrayImage = arrays.toArray(1);
        var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);
        var sdImage = ee.Image(eigenValues.sqrt())
        .arrayProject([0]).arrayFlatten([getNewBandNames('sd')]);
        return principalComponents.arrayProject([0])
        .arrayFlatten([getNewBandNames('pc')])
        .divide(sdImage);
        };
    
      var pcImage = getPrincipalComponents(centered, scale, region);
      return ee.Image(image.addBands(pcImage));
    });
    print("PCA imagery: ",Preped);