On Wed, Nov 5, 2008 at 5:39 AM, Nikos Alexandris <[EMAIL PROTECTED]> 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([EMAIL PROTECTED]) > 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([EMAIL PROTECTED]) > > 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 _______________________________________________ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user