[GRASS-stats] Raster maps from GRASS to R and back to GRASS?

[Nikos Alexandris]

Hi GRASS&R experts,

After playing around with R's functions for PCA, I need to

1. load GRASS raster maps in R
2. Do something with the data
   -involved workaround for NA's
3. Get back the NA's
4. Write back the transformed data in GRASS as new raster maps.


So, this is rather a (general?) GRASS - R question: how is it done,
after having calculated the principal components for a given set of
GRASS raster maps in R, to create new GRASS raster maps with them? I
mean how is it done to write the transformed values back to the right
coordinates?

Is there a clear example on the web using the writeRAST6()function?


Below an example.
Regards, Nikos
-------------------------------------------------------------------------

 # launch R from within GRASS
R


 # load spgrass6, and projection parameters
library(spgrass6)
G <- gmeta6()


 # load data
modis_2006_250m <- readRAST6 (c('MOD2006_239_250_sur_refl_b01',
'MOD2006_239_250_sur_refl_b02'))


 # workaround NA's as explained by Dylan [1]
modis_2006_250m.values <- which(!is.na(modis_2006_250m at data
$MOD2006_239_250_sur_refl_b01) & !is.na(modis_2006_250m at data
$MOD2006_239_250_sur_refl_b02))


 # create object without NA's
mod06_250m <- modis_2006_250m at data[modis_2006_250m.values, ]


 # perform pca based on SVD (default with center=TRUE, scale=FALSE)
prcomp.mod06_250m <- prcomp(mod06_250m)


 # check eigenvalues
summary(prcomp.mod06_250m)

Importance of components:
                           PC1      PC2
Standard deviation     3183.53 330.0584
Proportion of Variance    0.99   0.0106
Cumulative Proportion     0.99   1.0000


 # check eigenvectors (=rotation)
prcomp.mod06_250m$rotation

                                    PC1        PC2
MOD2006_239_250_sur_refl_b01 -0.3640966 -0.9313612
MOD2006_239_250_sur_refl_b02 -0.9313612  0.3640966


 # write back to GRASS the transformed data as new RASTER maps?
--> Any help here on how to get the values to the right coordinates? Is
it enough to just use some <--




=======================================================================
Extra info: comparing the resulted PC's from both R and GRASS
=======================================================================
 # compare with i.pca
i.pca input="MOD06_250_b01,MOD06_250_b02" output="TEST06_250m_b12"
rescale=0,0

Eigen values, (vectors), and [percent importance]:
   PC1 3481240.24 (  0.36  0.93 ) [ 98.92% ]
   PC2  37875.46 (  0.93 -0.36 ) [  1.08% ]


   ### Note ###
The signs of the eigenvectors exactly the opposite. Therefore we expect
that the Principal Components deriving from R's "prcomp()" will be the
opposite of the Components deriving from "i.pca"
   ### #### ###


 # get summary for PC's produced by "prcomp()"
summary(prcomp.mod06_250m$x)
      PC1                 PC2          
 Min.   :-10167.54   Min.   :-3919.99  
 1st Qu.: -3422.88   1st Qu.: -167.49  
 Median : -3119.27   Median :   45.32  
 Mean   : -3135.72   Mean   :   21.60  
 3rd Qu.: -2851.33   3rd Qu.:  239.13  
 Max.   :   -64.22   Max.   : 1611.15


 # get variances for PC's produce by "prcomp()"
var(prcomp.mod06_250m$x[, 1])
[1] 302107.8

var(prcomp.mod06_250m$x[, 2])
[1] 108471.9


 # compare with PC's from GRASS

### kept are only commonly reported parameters by GRASS' "r.univar" and
R's "summary()" ###


 # for GRASS' PC1
r.univar TEST2006_250m_b12.1

[...]
n: 350879
minimum: 64.2261
maximum: 10166.7
...
mean: 3139.65
mean of absolute values: 3139.65
standard deviation: 543.404
variance: 295287


 # for GRASS' PC2
r.univar TEST2006_250m_b12.2

 100%
total null and non-null cells: 1023328
total null cells: 672449

Of the non-null cells:
----------------------
n: 350879
minimum: -1609.93
maximum: 3922.17

mean: -22.0414
...
variance: 108444

######## Result is #########
### almost identical :-) ###
############################