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

Dylan Beaudette debeaudette at ucdavis.edu
Mon Apr 13 11:18:19 EDT 2009


On Monday 13 April 2009, Nikos Alexandris wrote:
> GRASS_&_R-experts,
> I really need some hand-of-help here.

Hi Nilkos!

> 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.

I think that you are going to have to maintain a band-wise null mask in R, as 
opposed to a single null mask for the entire set of bands.

> 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

It could be. The way that sp objects and their 'bands' (columns in the 
attached dataframe) function makes for some tricky NULL handling.

One thing you might consider is generating a single NULL mask that represents  
NULL values that occur in any, or all bands -- akin to an OR operation. I 
think that na.omit(sp_object at data) may be able to return just the 'complete' 
data, along with a list of indices where NA values have been removed.

The approach listed on the link you posted:
x.vals <- which( !is.na(x at data$r1) & !is.na(x at data$r2) & !is.na(x at data$r3) 
& !is.na(x at data$r4))

will only generate the non-NULL maks correctly when all bands have the same 
pattern of NULL values. Maybe I should add that hint.

Hope that helps,

Dylan




> > 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)
>
> _______________________________________________
> grass-stats mailing list
> grass-stats at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-stats



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


More information about the grass-stats mailing list