[GRASS-stats] Raster maps from GRASS to R and back to GRASS?
Dylan Beaudette
debeaudette at ucdavis.edu
Tue Mar 3 11:53:43 EST 2009
On Tuesday 03 March 2009, Nikos Alexandris wrote:
> 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?
>
Does this help:
http://casoilresource.lawr.ucdavis.edu/drupal/node/664
Cheers,
Dylan
> 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 :-) ###
> ############################
>
> _______________________________________________
> grass-stats mailing list
> grass-stats at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-stats
--
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341
More information about the grass-stats
mailing list