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

[Nikos Alexandris]

On Tue, 2008-10-28 at 09:39 -0700, Dylan Beaudette wrote:
> 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

[...]

Hi Dylan!

Sorry for replying so late and thank you for your detailed post. I now
see what I was doing wrong when trying to produce pc's with R. Here is
how I try to replicate your example using my data:

#load library(spgrass6) and data
library(spgrass6)
x <-
readRAST6(c('MOD2007_242_500_sur_refl_b02','MOD2007_242_500_sur_refl_b06','MOD2007_242_500_sur_refl_b07'))

# prcomp (gives me an error) ## maybe because of NULL's ??
x.pca <- prcomp(x at data)
Error in svd(x, nu = 0) : infinite or missing values in 'x'

----

# Trying with a subset with no NULL's ## resolution is set to 250m
i.pca in=b2,b6,b7 out=test_pca_b267

Eigen values:
( 0.79 0.61 0.06 )
( 0.44 -0.50 -0.74 )
( 0.42 -0.61 0.67 )

# set resolution to 500m and repeat i.pca
g.region res=500 -pa
i.pca in=b2,b6,b7 out=test_pca_b267_500m

Eigen values:
( 0.79 0.61 0.06 )
( 0.44 -0.50 -0.74 )
( 0.42 -0.61 0.67 )

## Well, this subset is rather small, thus no differences between 250m
and 500m (?)
g.region -p

[...]
cells:      7455
[...]

# trying with R
x <- readRAST6(c('b2','b6','b7'))
x.pca <- prcomp(x at data)

Error in svd(x, nu = 0) : infinite or missing values in 'x'


## again the same error!

I'll try another time again.

Kind regards, Nikos