2

I am currently trying to obtain two rasters that are the Principal component axes 1 and 2 for a list of bioclimatic variables that I am using. I am using the rasterPCA() function that exists in the RStoolbox in R, but I am unable to write it to a raster. Below is the code and the error that I am getting:

library(raster)
library(RStoolbox)

##Loading the variables and then creating a raster stack, followed by a PCA
bio3 <- raster("C:\\Users\\rameshv\\bio3")
bio4 <- raster("C:\\Users\\rameshv\\bio4")
bio5 <- raster("C:\\Users\\rameshv\\bio5")
bio6 <- raster("C:\\Users\\rameshv\\bio6")
bio14 <- raster("C:\\Users\\rameshv\\bio_14")

pres_stack <- stack(bio3, bio4,bio5,bio6,bio14)

pre_pca <-  rasterPCA(pres_stack, nComp = 2)

writeRaster(pre_pca, "C:\\Users\\rameshv\\PC.tif")

Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘writeRaster’ for signature ‘"rasterPCA", "character"’

I even tried a modification of the above code

writeRaster(pre_pca,"C:\\Users\\rameshv\\4_PCAforR\\PC.asc", format="ascii", bylayer=T)

Getting the same error as above: Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘writeRaster’ for signature ‘"rasterPCA", "character"’

Any suggestions?

Vijay Ramesh
  • 1,332
  • 15
  • 46
  • The out put of the function is a list object that contains the rasters. You can have the function result in a raster brick my using pre_pca <- rasterPCA(pres_stack, nComp = 2)$map – Jeffrey Evans Mar 07 '17 at 20:03
  • Thanks @JeffreyEvans. That is definitely a quicker workaround. Failed to notice that it is a list object. – Vijay Ramesh Mar 07 '17 at 20:05
  • @JeffreyEvans Could I nudge you towards this? This query was based off of one of your suggestions: http://gis.stackexchange.com/questions/230970/how-many-polygons-does-a-single-polygon-intersect-with – Vijay Ramesh Mar 07 '17 at 20:26

1 Answers1

3

After digging in to the RSToolbox package, I realized that the rasterPCA() command writes in a "large rasterPCA" that consists of 3 elements and I failed to call the right element.

Shown below is the output from using the rasterPCA() command:

> pre_pca
  $call
  rasterPCA(img = pres_stack, nComp = 2)

  $model
  Call:
  princomp(cor = spca, covmat = covMat[[1]])

  Standard deviations:
  Comp.1       Comp.2       Comp.3       Comp.4       Comp.5       Comp.6  
  113.53111875  22.39603549   8.08710815   4.11482299   2.21497438   0.79337709    

  6  variables and  83300 observations.

  $map
  class       : RasterBrick 
  dimensions  :   (nrow, ncol, ncell, nlayers)
  resolution  : 0.04166667, 0.04166667  (x, y)
  extent      :  (xmin, xmax, ymin, ymax)
  coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 
  data source : in memory
  names       :        PC1,        PC2 
  min values  : -870.80335,  -49.97994 
  max values  :  141.77369,   75.40234 


 attr(,"class")
 [1] "rasterPCA" "RStoolbox"

By calling $map when using the writeRaster command does the trick.

writeRaster(pre_pca$map,"C:\\Users\\rameshv\\4_PCAforR\\PC.asc", format="ascii", bylayer=T)

Now I have an output of two rasters with a .asc extension that can be visualized in ArcMap now.

Vijay Ramesh
  • 1,332
  • 15
  • 46