[GRASS-stats] Raster maps from GRASS to R and back to GRASS?

Nikos Alexandris nikos.alexandris at felis.uni-freiburg.de
Mon Apr 13 09:21:07 EDT 2009


GRASS_&_R-experts,
I really need some hand-of-help here.

What I actually do is nothing but a PCA on a combination of different
bands. The fact is that bands a, b, c are of equal "length" in R and
have _more_ non-NULLs than bands d, e, f.

That is, bands d, e and f are of equal length compared to each other but
of different length comparing to bands a, b and c.

I suspect this is the root of the problem. Any suggestion?
Thanks, Nikos

> Nikos:
> > I need a confirmation that I am doing the NA-work-around the _right_
> > way.
> 
> Sorry, I fat-fingered the Send button before. For simplicity I tried to
> use small names instead of the long "original" names. Anyhow, I think
> it's clear what I am doing and what the problem is.
> 
> 
> The problem is:
> > I get several NULL's in the raster maps I write back in GRASS in
> > locations (=pixels) were the input raster map(s) have a NON-NULL value.
> 
> For simplicity I tried to use small names instead of the long "original"
> names. Anyhow, I think it's clear what I am doing and what the problem
> is.
> 
> Cheers,  Nikos
> 
> # launch grass; launch R from within grass
> # load libraries and projection information
> library(spgrass6); G <- gmeta6()
> 
> # load data
> raw <- readRAST6 (c('a', 'b', 'c', 'd', 'e', 'f'))
> 
> # create vector pointing to values
> values <- which(!is.na(raw at data$a) & !is.na(raw at data$b) & !
> is.na(raw at data$c) & !is.na(raw at data$d) & !is.na(raw at data$e) & !
> is.na(raw at data$f))
> 
> # get data
> data <- raw at data[values, ]
> 
> # check structure
> str(data)
> 
> # perform pca with prcomp(): SVD - SCALED
> prcomp.data <- prcomp(data, scale=TRUE)
> prcomp.data$rotation
> summary(prcomp.data)
> 
> # add new columns to the original sp object #
> raw at data$pc1 <- NA
> raw at data$pc2 <- NA
> raw at data$pc3 <- NA
> raw at data$pc4 <- NA
> raw at data$pc5 <- NA
> raw at data$pc6 <- NA
> str(raw)
> 
> # fill-in the PCs
> raw at data$pc1[values] <- prcomp.data$x[,"PC1"]
> raw at data$pc2[values] <- prcomp.data$x[,"PC1"]
> raw at data$pc3[values] <- prcomp.data$x[,"PC1"]
> raw at data$pc4[values] <- prcomp.data$x[,"PC1"]
> raw at data$pc5[values] <- prcomp.data$x[,"PC1"]
> raw at data$pc6[values] <- prcomp.data$x[,"PC1"]
> 
> # write grass raster maps
> writeRAST6(raw, zcol=7, vname="prcomp1_0607", overwrite=TRUE)
> writeRAST6(raw, zcol=8, vname="prcomp1_0607", overwrite=TRUE)
> writeRAST6(raw, zcol=9, vname="prcomp1_0607", overwrite=TRUE)
> writeRAST6(raw, zcol=10, vname="prcomp1_0607", overwrite=TRUE)
> writeRAST6(raw, zcol=11, vname="prcomp1_0607", overwrite=TRUE)
> writeRAST6(raw, zcol=12, vname="prcomp1_0607", overwrite=TRUE)



More information about the grass-stats mailing list