[GRASS-user] i.pca vs. r.covar/m.eigensystem/r.mapcalc

[Dylan Beaudette]

On Monday 27 October 2008, Nikos Alexandris wrote:
> My apologies for BUMPing :-)
>
> On Sun, 2008-10-12 at 00:20 +0200, Nikos Alexandris wrote:
> > I wanted to check/verify some i.pca results so I calculated 2 principal
> > components using r.covar, m.eigensystem and r.mapcalc. The results are
> > not the same with i.pca's output.
> >
> > (Maybe this post http://www.nabble.com/Eigen-vectors-tt15686112.html was
> > about the same issue)
> >
> > # i.pca performed on (3) modis surface reflectance bands (2, 6 and 7  -
> > 500m spatial resolution)
>
> *** and region resolution set to 500m of course***
>
> > r.info pca.2007.242.500.b267.1 -h
> > [...]
> > Eigen vectors:
> >    ( -0.64 -0.65 -0.42 )
> >    ( -0.71 0.29 0.64 )
> >    ( 0.29 -0.70 0.65 )
> > [...]
>
> I repeated the same i.pca but setting this time a different region
> resolution (250m) than before (500m) and I get slightly different
> eigenvectors:
>
> *** with region resolution=250m ***
> Eigen values:
> ( -0.64 -0.65 -0.42 )
> ( -0.71 0.28 0.64 )
> ( -0.30 0.71 -0.64 )
>
>
> I understand that the region "resolution" influences the results of
> i.pca. OK, then how would someone re-produce the same PCs doing the math
> in R and not in GRASS?
>
> Could someone explain or give any hints to learn more about i.pca?
>
> Thank you, Nikos
>

Hi,

Remember that the region settings may cause re-sampling of your image due to 
changes in resolution or grid-cell alignment. Here is a re-enactment of your 
example above in GRASS, then in R- using some color imagery. Although I am 
not an expert in PCA, I think that the methods in R/GRASS here are 
comparable. In the output from R, the row and column names of the returned 
matrix are printed, and it looks like the rotation vectors (eigen vectors) 
for the 3rd PC are the least stable. This seems to make sense as (in this 
case) the 3rd PC accounts for the least amount of variance-- and is thus 
probably just picking up noise.

Not sure about applying the rotation manually and how that compares with 
i.pca ...

Cheers,

Dylan


------------------- example ----------------------

# intial test at 1m res
g.region res=1 -a
i.pca in=naip.red,naip.green,naip.blue out=pca --o
d.rgb r=pca.2 g=pca.1 b=pca.3

Eigen values:
( -0.74 -0.55 -0.38 )
( 0.67 -0.67 -0.31 )
( 0.08 0.49 -0.87 )

# now try 10m res
g.region res=10 -a
i.pca in=naip.red,naip.green,naip.blue out=pca --o
d.rgb r=pca.2 g=pca.1 b=pca.3

Eigen values:
( -0.73 -0.56 -0.38 )
( 0.68 -0.65 -0.34 )
( 0.06 0.51 -0.86 )


## now try it in R

# init libraries
library(spgrass6)

# read in same data used in GRASS PCA calc
system('g.region res=1 -a')
x <- readRAST6(c('naip.red','naip.green', 'naip.blue'))

# use prcomp, as it directly returns rotation matrix
# same as Eigen values?
x.pca <- prcomp(x at data)

# transpose result for printing in same order as GRASS
signif(t(x.pca$rotation), 2)

    naip.red naip.green naip.blue
PC1   -0.730      -0.56     -0.38
PC2    0.680      -0.62     -0.39
PC3    0.023       0.54     -0.84



# do again, at reduces res:
system('g.region res=10 -a')

x <- readRAST6(c('naip.red','naip.green', 'naip.blue'))

# use prcomp, as it directly returns rotation matrix
# same as Eigen values- yes, see the manual page for prcomp
x.pca <- prcomp(x at data)

# transpose result for printing in same order as GRASS
signif(t(x.pca$rotation), 2)

PC1   -0.730      -0.56     -0.38
PC2    0.680      -0.64     -0.36
PC3    0.041       0.53     -0.85







-- 
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341