PCA using raster datasets in R

Question

I have several large rasters that I want to process in a PCA (to produce summary rasters). I have seen several examples whereby people seem to be simply calling prcomp or princomp. However, when I do this, I get the following error message:

Error in as.vector(data): no method for coercing this S4 class to a vector

Example code:

files<-list.files() # a set of rasters
layers<-stack(files) # using the raster package
pca<-prcomp(layers)

I have tried using a raster brick instead of stack but that doesn't seem to the issue. What method do I need to provide the command so that it can convert the raster data to vector format? I understand that there are ways to sample the raster and run the PCA from that, but I would really like to understand why the above method is not working.

Thanks!

Answer 1

The above method is not working simply because prcomp does not know how to deal with a raster object. It only knows how to deal with vectors, and coercing to vector does not work, hence the error.

What you need to do is read each of your files into a vector, and put each of the rasters in a column of a matrix. Each row will then be a time series of values at a single spatial location, and each column will be all the pixels at a certain time step. Note that the exact spatial coordinates are not needed in this approach. This matrix serves as the input of prcomp.

Reading the files can be done using readGDAL, and using as.data.frame to cast the spatial data to data.frame.

Answer 2

Answer to my own question: I ended up doing something slightly different: rather than using every raster cell as input (very large dataset), I took a sample of points, ran the PCA and then saved the output model so that I could make predictions for each grid cell…maybe not the best solution but it works:

rasters <- stack(myRasters)

sr <- sampleRandom(rasters, 5000) # sample 5000 random grid cells

# run PCA on random sample with correlation matrix
# retx=FALSE means don't save PCA scores 
pca <- prcomp(sr, scale=TRUE, retx=FALSE) 

# write PCA model to file 
dput(pca, file=paste("./climate/", name, "/", name, "_pca.csv", sep=""))

x <- predict(rasters, pca, index=1:6) # create new rasters based on PCA predictions