[GRASS-user] i.pca vs. r.covar/m.eigensystem/r.mapcalc
Dylan Beaudette
dylan.beaudette at gmail.com
Wed Nov 5 10:06:12 EST 2008
On Wed, Nov 5, 2008 at 5:39 AM, Nikos Alexandris
<nikos.alexandris at felis.uni-freiburg.de> wrote:
> 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
>
Sounds like NAs are getting in the way. This is a limitation of
several R functions, that becomes particularly annoying when working
with GIS data. Here are some tips on dealing with NA in GIS data:
http://casoilresource.lawr.ucdavis.edu/drupal/node/664
Cheers,
Dylan
More information about the grass-user
mailing list