Raster maps from GRASS to R and back to GRASS?

[Nikos Alexandris]

Hi GRASS&R experts,

After playing around with R's functions for PCA, I need to

1. load GRASS raster maps in R
2. Do something with the data
   -involved workaround for NA's
3. Get back the NA's
4. Write back the transformed data in GRASS as new raster maps.


So, this is rather a (general?) GRASS - R question: how is it done,
after having calculated the principal components for a given set of
GRASS raster maps in R, to create new GRASS raster maps with them? I
mean how is it done to write the transformed values back to the right
coordinates?

Is there a clear example on the web using the writeRAST6()function?


Below an example.
Regards, Nikos
-------------------------------------------------------------------------

 # launch R from within GRASS
R


 # load spgrass6, and projection parameters
library(spgrass6)
G <- gmeta6()


 # load data
modis_2006_250m <- readRAST6 (c('MOD2006_239_250_sur_refl_b01',
'MOD2006_239_250_sur_refl_b02'))


 # workaround NA's as explained by Dylan [1]
modis_2006_250m.values <- which(!is.na(modis_2006_250m@data
$MOD2006_239_250_sur_refl_b01) & !is.na(modis_2006_250m@data
$MOD2006_239_250_sur_refl_b02))


 # create object without NA's
mod06_250m <- modis_2006_250m@data[modis_2006_250m.values, ]


 # perform pca based on SVD (default with center=TRUE, scale=FALSE)
prcomp.mod06_250m <- prcomp(mod06_250m)


 # check eigenvalues
summary(prcomp.mod06_250m)

Importance of components:
                           PC1      PC2
Standard deviation     3183.53 330.0584
Proportion of Variance    0.99   0.0106
Cumulative Proportion     0.99   1.0000


 # check eigenvectors (=rotation)
prcomp.mod06_250m$rotation

                                    PC1        PC2
MOD2006_239_250_sur_refl_b01 -0.3640966 -0.9313612
MOD2006_239_250_sur_refl_b02 -0.9313612  0.3640966


 # write back to GRASS the transformed data as new RASTER maps?
--> Any help here on how to get the values to the right coordinates? Is
it enough to just use some <--




=======================================================================
Extra info: comparing the resulted PC's from both R and GRASS
=======================================================================
 # compare with i.pca
i.pca input="MOD06_250_b01,MOD06_250_b02" output="TEST06_250m_b12"
rescale=0,0

Eigen values, (vectors), and [percent importance]:
   PC1 3481240.24 (  0.36  0.93 ) [ 98.92% ]
   PC2  37875.46 (  0.93 -0.36 ) [  1.08% ]


   ### Note ###
The signs of the eigenvectors exactly the opposite. Therefore we expect
that the Principal Components deriving from R's "prcomp()" will be the
opposite of the Components deriving from "i.pca"
   ### #### ###


 # get summary for PC's produced by "prcomp()"
summary(prcomp.mod06_250m$x)
      PC1                 PC2          
 Min.   :-10167.54   Min.   :-3919.99  
 1st Qu.: -3422.88   1st Qu.: -167.49  
 Median : -3119.27   Median :   45.32  
 Mean   : -3135.72   Mean   :   21.60  
 3rd Qu.: -2851.33   3rd Qu.:  239.13  
 Max.   :   -64.22   Max.   : 1611.15


 # get variances for PC's produce by "prcomp()"
var(prcomp.mod06_250m$x[, 1])
[1] 302107.8

var(prcomp.mod06_250m$x[, 2])
[1] 108471.9


 # compare with PC's from GRASS

### kept are only commonly reported parameters by GRASS' "r.univar" and
R's "summary()" ###


 # for GRASS' PC1
r.univar TEST2006_250m_b12.1

[...]
n: 350879
minimum: 64.2261
maximum: 10166.7
...
mean: 3139.65
mean of absolute values: 3139.65
standard deviation: 543.404
variance: 295287


 # for GRASS' PC2
r.univar TEST2006_250m_b12.2

 100%
total null and non-null cells: 1023328
total null cells: 672449

Of the non-null cells:
----------------------
n: 350879
minimum: -1609.93
maximum: 3922.17

mean: -22.0414
...
variance: 108444

######## Result is #########
### almost identical :-) ###
############################

_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Dylan Beaudette]

On Tuesday 03 March 2009, Nikos Alexandris wrote:

> Hi GRASS&R experts,
>
> After playing around with R's functions for PCA, I need to
>
> 1. load GRASS raster maps in R
> 2. Do something with the data
>    -involved workaround for NA's
> 3. Get back the NA's
> 4. Write back the transformed data in GRASS as new raster maps.
>
>
> So, this is rather a (general?) GRASS - R question: how is it done,
> after having calculated the principal components for a given set of
> GRASS raster maps in R, to create new GRASS raster maps with them? I
> mean how is it done to write the transformed values back to the right
> coordinates?
>
> Is there a clear example on the web using the writeRAST6()function?
>

Does this help:
http://casoilresource.lawr.ucdavis.edu/drupal/node/664

Cheers,

Dylan

> Below an example.
> Regards, Nikos
> -------------------------------------------------------------------------
>
>  # launch R from within GRASS
> R
>
>
>  # load spgrass6, and projection parameters
> library(spgrass6)
> G <- gmeta6()
>
>
>  # load data
> modis_2006_250m <- readRAST6 (c('MOD2006_239_250_sur_refl_b01',
> 'MOD2006_239_250_sur_refl_b02'))
>
>
>  # workaround NA's as explained by Dylan [1]
> modis_2006_250m.values <- which(!is.na(modis_2006_250m@data
> $MOD2006_239_250_sur_refl_b01) & !is.na(modis_2006_250m@data
> $MOD2006_239_250_sur_refl_b02))
>
>
>  # create object without NA's
> mod06_250m <- modis_2006_250m@data[modis_2006_250m.values, ]
>
>
>  # perform pca based on SVD (default with center=TRUE, scale=FALSE)
> prcomp.mod06_250m <- prcomp(mod06_250m)
>
>
>  # check eigenvalues
> summary(prcomp.mod06_250m)
>
> Importance of components:
>                            PC1      PC2
> Standard deviation     3183.53 330.0584
> Proportion of Variance    0.99   0.0106
> Cumulative Proportion     0.99   1.0000
>
>
>  # check eigenvectors (=rotation)
> prcomp.mod06_250m$rotation
>
>                                     PC1        PC2
> MOD2006_239_250_sur_refl_b01 -0.3640966 -0.9313612
> MOD2006_239_250_sur_refl_b02 -0.9313612  0.3640966
>
>
>  # write back to GRASS the transformed data as new RASTER maps?
> --> Any help here on how to get the values to the right coordinates? Is
> it enough to just use some <--
>
>
>
>
> =======================================================================
> Extra info: comparing the resulted PC's from both R and GRASS
> =======================================================================
>  # compare with i.pca
> i.pca input="MOD06_250_b01,MOD06_250_b02" output="TEST06_250m_b12"
> rescale=0,0
>
> Eigen values, (vectors), and [percent importance]:
>    PC1 3481240.24 (  0.36  0.93 ) [ 98.92% ]
>    PC2  37875.46 (  0.93 -0.36 ) [  1.08% ]
>
>
>    ### Note ###
> The signs of the eigenvectors exactly the opposite. Therefore we expect
> that the Principal Components deriving from R's "prcomp()" will be the
> opposite of the Components deriving from "i.pca"
>    ### #### ###
>
>
>  # get summary for PC's produced by "prcomp()"
> summary(prcomp.mod06_250m$x)
>       PC1                 PC2
>  Min.   :-10167.54   Min.   :-3919.99
>  1st Qu.: -3422.88   1st Qu.: -167.49
>  Median : -3119.27   Median :   45.32
>  Mean   : -3135.72   Mean   :   21.60
>  3rd Qu.: -2851.33   3rd Qu.:  239.13
>  Max.   :   -64.22   Max.   : 1611.15
>
>
>  # get variances for PC's produce by "prcomp()"
> var(prcomp.mod06_250m$x[, 1])
> [1] 302107.8
>
> var(prcomp.mod06_250m$x[, 2])
> [1] 108471.9
>
>
>  # compare with PC's from GRASS
>
> ### kept are only commonly reported parameters by GRASS' "r.univar" and
> R's "summary()" ###
>
>
>  # for GRASS' PC1
> r.univar TEST2006_250m_b12.1
>
> [...]
> n: 350879
> minimum: 64.2261
> maximum: 10166.7
> ...
> mean: 3139.65
> mean of absolute values: 3139.65
> standard deviation: 543.404
> variance: 295287
>
>
>  # for GRASS' PC2
> r.univar TEST2006_250m_b12.2
>
>  100%
> total null and non-null cells: 1023328
> total null cells: 672449
>
> Of the non-null cells:
> ----------------------
> n: 350879
> minimum: -1609.93
> maximum: 3922.17
>
> mean: -22.0414
> ...
> variance: 108444
>
> ######## Result is #########
> ### almost identical :-) ###
> ############################
>
> _______________________________________________
> grass-stats mailing list
> [hidden email]
> 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
_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Nikos Alexandris]

On Tue, 2009-03-03 at 08:53 -0800, Dylan Beaudette wrote:
>
> Does this help:
> http://casoilresource.lawr.ucdavis.edu/drupal/node/664
>
> Cheers,
>
> Dylan

:DDDDDD  Amazing, I need to rest. How in the world did I forget that
when I wanted (but forgot) to cite the same page for the workaround??

:DDD

Cheers, Nikos

_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Nikos Alexandris]

Hi all!


Dylan:
> > Does this help:
> > http://casoilresource.lawr.ucdavis.edu/drupal/node/664

I need a confirmation that I am doing the NA-work-around the _right_
way. 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.


# 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@data$a) & !is.na(raw@data$b) & !
is.na(raw@data$c) & !is.na(raw@data$d) & !is.na(raw@data$e) & !
is.na(raw@data$f))

# get data
data <- raw@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@data$pc1 <- NA
raw@data$pc2 <- NA
raw@data$pc3 <- NA
raw@data$pc4 <- NA
raw@data$pc5 <- NA
raw@data$pc6 <- NA
str(raw)

# fill-in the PCs
raw@data$pc1[values] <- prcomp.data$x[,"PC1"]
raw@data$pc1[values] <- prcomp.data$x[,"PC1"]
raw@data$pc1[values] <- prcomp.data$x[,"PC1"]
raw@data$pc1[values] <- prcomp.data$x[,"PC1"]
raw@data$pc1[values] <- prcomp.data$x[,"PC1"]
mod_prepost_b267.raw@data$pc2[mod_prepost_b267.values] <-
prcomp.modis0607.nonas$x[,"PC2"]
mod_prepost_b267.raw@data$pc3[mod_prepost_b267.values] <-
prcomp.modis0607.nonas$x[,"PC3"]
mod_prepost_b267.raw@data$pc4[mod_prepost_b267.values] <-
prcomp.modis0607.nonas$x[,"PC4"]
mod_prepost_b267.raw@data$pc5[mod_prepost_b267.values] <-
prcomp.modis0607.nonas$x[,"PC5"]
mod_prepost_b267.raw@data$pc6[mod_prepost_b267.values] <-
prcomp.modis0607.nonas$x[,"PC6"]
str(mod_prepost_b267.raw)

# write grass raster maps
writeRAST6(mod_prepost_b267.raw, zcol=7, vname="prcomp1_0607",
overwrite=TRUE)
writeRAST6(mod_prepost_b267.raw, zcol=7, vname="prcomp2_0607",
overwrite=TRUE)
writeRAST6(mod_prepost_b267.raw, zcol=7, vname="prcomp3_0607",
overwrite=TRUE)
writeRAST6(mod_prepost_b267.raw, zcol=7, vname="prcomp4_0607",
overwrite=TRUE)
writeRAST6(mod_prepost_b267.raw, zcol=7, vname="prcomp5_0607",
overwrite=TRUE)
writeRAST6(mod_prepost_b267.raw, zcol=7, vname="prcomp6_0607",
overwrite=TRUE)

_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Nikos Alexandris]

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@data$a) & !is.na(raw@data$b) & !
is.na(raw@data$c) & !is.na(raw@data$d) & !is.na(raw@data$e) & !
is.na(raw@data$f))

# get data
data <- raw@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@data$pc1 <- NA
raw@data$pc2 <- NA
raw@data$pc3 <- NA
raw@data$pc4 <- NA
raw@data$pc5 <- NA
raw@data$pc6 <- NA
str(raw)

# fill-in the PCs
raw@data$pc1[values] <- prcomp.data$x[,"PC1"]
raw@data$pc2[values] <- prcomp.data$x[,"PC1"]
raw@data$pc3[values] <- prcomp.data$x[,"PC1"]
raw@data$pc4[values] <- prcomp.data$x[,"PC1"]
raw@data$pc5[values] <- prcomp.data$x[,"PC1"]
raw@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
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Nikos Alexandris]

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@data$a) & !is.na(raw@data$b) & !
> is.na(raw@data$c) & !is.na(raw@data$d) & !is.na(raw@data$e) & !
> is.na(raw@data$f))
>
> # get data
> data <- raw@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@data$pc1 <- NA
> raw@data$pc2 <- NA
> raw@data$pc3 <- NA
> raw@data$pc4 <- NA
> raw@data$pc5 <- NA
> raw@data$pc6 <- NA
> str(raw)
>
> # fill-in the PCs
> raw@data$pc1[values] <- prcomp.data$x[,"PC1"]
> raw@data$pc2[values] <- prcomp.data$x[,"PC1"]
> raw@data$pc3[values] <- prcomp.data$x[,"PC1"]
> raw@data$pc4[values] <- prcomp.data$x[,"PC1"]
> raw@data$pc5[values] <- prcomp.data$x[,"PC1"]
> raw@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
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Dylan Beaudette]

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@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@data$r1) & !is.na(x@data$r2) & !is.na(x@data$r3)
& !is.na(x@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@data$a) & !is.na(raw@data$b) & !
> > is.na(raw@data$c) & !is.na(raw@data$d) & !is.na(raw@data$e) & !
> > is.na(raw@data$f))
> >
> > # get data
> > data <- raw@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@data$pc1 <- NA
> > raw@data$pc2 <- NA
> > raw@data$pc3 <- NA
> > raw@data$pc4 <- NA
> > raw@data$pc5 <- NA
> > raw@data$pc6 <- NA
> > str(raw)
> >
> > # fill-in the PCs
> > raw@data$pc1[values] <- prcomp.data$x[,"PC1"]
> > raw@data$pc2[values] <- prcomp.data$x[,"PC1"]
> > raw@data$pc3[values] <- prcomp.data$x[,"PC1"]
> > raw@data$pc4[values] <- prcomp.data$x[,"PC1"]
> > raw@data$pc5[values] <- prcomp.data$x[,"PC1"]
> > raw@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
> [hidden email]
> 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
_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Roger Bivand]

On Mon, 13 Apr 2009, Dylan Beaudette wrote:

> 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@data) may be able to return just the 'complete'
> data, along with a list of indices where NA values have been removed.

Would simply using complete.cases(x@data) be what you are looking for?

Roger

>
> The approach listed on the link you posted:
> x.vals <- which( !is.na(x@data$r1) & !is.na(x@data$r2) & !is.na(x@data$r3)
> & !is.na(x@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@data$a) & !is.na(raw@data$b) & !
>>> is.na(raw@data$c) & !is.na(raw@data$d) & !is.na(raw@data$e) & !
>>> is.na(raw@data$f))
>>>
>>> # get data
>>> data <- raw@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@data$pc1 <- NA
>>> raw@data$pc2 <- NA
>>> raw@data$pc3 <- NA
>>> raw@data$pc4 <- NA
>>> raw@data$pc5 <- NA
>>> raw@data$pc6 <- NA
>>> str(raw)
>>>
>>> # fill-in the PCs
>>> raw@data$pc1[values] <- prcomp.data$x[,"PC1"]
>>> raw@data$pc2[values] <- prcomp.data$x[,"PC1"]
>>> raw@data$pc3[values] <- prcomp.data$x[,"PC1"]
>>> raw@data$pc4[values] <- prcomp.data$x[,"PC1"]
>>> raw@data$pc5[values] <- prcomp.data$x[,"PC1"]
>>> raw@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
>> [hidden email]
>> http://lists.osgeo.org/mailman/listinfo/grass-stats
>
>
>
>

--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: [phone]; fax [phone]
e-mail: [hidden email]

_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Nikos Alexandris]

[
HiYo Dylan ;-)
]

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


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

If I understand your suggestion, you mean to:

# load the bands as usual (i.e. combine a, b, c, d, e, f in one
SpatialGridDataFrame)

# create NULL-"masks" (=vectors pointing to NULLs) band-wise

# extract data of interest in a normal "data.frame" and do something
with it

# add new "columns" in the original SpatialGridDataFrame

# fill-in the results by using the per-band NULL-masks

The question that immediately pops-up is: How should I decide which
NULL-mask to use for each PC since PCA redistributes the original data?

Or is this irrelevant?

:-O


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

Hmmm... (thinking), ok, I am close... I think (!?)


> I think that na.omit(sp_object@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@data$r1) & !is.na(x@data$r2) & !is.na(x@data$r3)
> & !is.na(x@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.

What happens if I just replace the & with an | (=OR) operator? Isn't it
the same? (I'll try it out).


> Hope that helps,
> Dylan

Of course it helps (me). Discussions like this one are (for me) the
essence of the MLs :D

Cheers, Nikos

_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Nikos Alexandris]

Nikos:
> >> 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.
[...]


Dylan:
> > 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@data) may be able to return just the 'complete'
> > data, along with a list of indices where NA values have been removed.


Roger:
> Would simply using complete.cases(x@data) be what you are looking for?


So, this is the _correct and easy_ way instead of using a couple of |
(=OR) operators in my case? Or am I misinterpreting it?

Kindest regards, Nikos

_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Dylan Beaudette]

On Monday 13 April 2009, Nikos Alexandris wrote:

> Nikos:
> > >> 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.
>
> [...]
>
> Dylan:
> > > 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@data) may be able to return
> > > just the 'complete' data, along with a list of indices where NA values
> > > have been removed.
>
> Roger:
> > Would simply using complete.cases(x@data) be what you are looking for?
>
> So, this is the _correct and easy_ way instead of using a couple of |
> (=OR) operators in my case? Or am I misinterpreting it?
>
> Kindest regards, Nikos

Hi.

As usualy Roger's suggestion seems to be the simplest, most direct approach.
It seems like a more general solution, as opposed to a sequence of
is.na() | ... one for each raster map.

Dylan

--
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341
_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Dylan Beaudette]

On Monday 13 April 2009, Nikos Alexandris wrote:

> [
> HiYo Dylan ;-)
> ]
>
> Nikos:
> > > 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.
>
> Dylan:
> > 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.
>
> If I understand your suggestion, you mean to:
>
> # load the bands as usual (i.e. combine a, b, c, d, e, f in one
> SpatialGridDataFrame)
>
> # create NULL-"masks" (=vectors pointing to NULLs) band-wise
>
> # extract data of interest in a normal "data.frame" and do something
> with it
>
> # add new "columns" in the original SpatialGridDataFrame
>
> # fill-in the results by using the per-band NULL-masks
>
> The question that immediately pops-up is: How should I decide which
> NULL-mask to use for each PC since PCA redistributes the original data?
>
> Or is this irrelevant?
>
> :-O
> :
> > > 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.
>
> Hmmm... (thinking), ok, I am close... I think (!?)
>
> > I think that na.omit(sp_object@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@data$r1) & !is.na(x@data$r2) &
> > !is.na(x@data$r3) & !is.na(x@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.
>
> What happens if I just replace the & with an | (=OR) operator? Isn't it
> the same? (I'll try it out).
>
> > Hope that helps,
> > Dylan
>
> Of course it helps (me). Discussions like this one are (for me) the
> essence of the MLs :D
>
> Cheers, Nikos

I guess I wasn't clear. I should have said something like: "for input maps
that have different NULL value (spatial distributions), you will need to
generate a new NULL mask that can be used to filter out all cells that
overlap with 1 or more NULL values in your stack of images". Using
complete.cases() as Roger suggested may accomplish this.

An idea on how to apply this function:

# generate some data; sprinkle in some NA ; convert to DF
x <- rnorm(100)
x[sample(1:100, 5)] <- NA
x.mat <- round(matrix(x, ncol=5), 2)
x.df <- as.data.frame(x.mat)

# extract an NA mask
x.na_mask <- complete.cases(x.df)
 
# just the 'rows' with complete data (across bands)
 x.df[x.na_mask,]

      V1    V2    V3    V4    V5
1   1.32 -1.27  0.00 -0.87  1.52
2  -0.44  1.33  0.78 -1.27  0.00
5   2.70 -0.08 -0.77  0.38  0.32
6  -0.52 -0.07  0.13 -0.43  1.12
7  -0.63  2.02  0.45  0.48 -0.59
8  -0.52 -0.78 -0.59 -0.07 -0.08
11 -1.37 -1.15  0.23  0.73  0.07
12 -0.11 -0.66  0.50 -1.14  0.71
13 -0.38 -0.08  1.00 -0.88 -0.19
14 -1.14 -0.45 -1.37 -0.43 -2.18
15  1.92  0.46 -0.72  0.12  0.27
16  0.51  0.22 -0.39 -1.54 -1.04
17  1.14 -1.11  1.04  0.00 -1.11
18  0.53  0.62 -0.10 -0.17  0.15
19  0.84  0.27 -1.01  0.87  0.31
20  0.82  0.87  0.46  1.35  0.13


# all the data
x.df
      V1    V2    V3    V4    V5
1   1.32 -1.27  0.00 -0.87  1.52
2  -0.44  1.33  0.78 -1.27  0.00
3  -0.31 -0.13 -1.14  2.08    NA
4   0.76    NA -1.20 -0.54 -0.85
5   2.70 -0.08 -0.77  0.38  0.32
6  -0.52 -0.07  0.13 -0.43  1.12
7  -0.63  2.02  0.45  0.48 -0.59
8  -0.52 -0.78 -0.59 -0.07 -0.08
9  -0.38 -0.51  0.44    NA -0.81
10 -1.38    NA -0.65 -0.50    NA
11 -1.37 -1.15  0.23  0.73  0.07
12 -0.11 -0.66  0.50 -1.14  0.71
13 -0.38 -0.08  1.00 -0.88 -0.19
14 -1.14 -0.45 -1.37 -0.43 -2.18
15  1.92  0.46 -0.72  0.12  0.27
16  0.51  0.22 -0.39 -1.54 -1.04
17  1.14 -1.11  1.04  0.00 -1.11
18  0.53  0.62 -0.10 -0.17  0.15
19  0.84  0.27 -1.01  0.87  0.31
20  0.82  0.87  0.46  1.35  0.13

# generate something cool from complete data
library(MASS)
p <- princomp(~ . , data=x.df[x.na_mask,])

# assign back to original data
x.df$pca_1[x.na_mask] <-  predict(p)[,1]

This would of course be done the the @data slot on an sp object. Thanks to
Roger for this slick approach.

How about that!?

Cheers,
Dylan



--
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341
_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Roger Bivand]

On Tue, 14 Apr 2009, Nikos Alexandris wrote:

>
> Nikos:
>>>> 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.
> [...]
>
>
> Dylan:
>>> 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@data) may be able to return just the 'complete'
>>> data, along with a list of indices where NA values have been removed.
>
>
> Roger:
>> Would simply using complete.cases(x@data) be what you are looking for?
>
>
> So, this is the _correct and easy_ way instead of using a couple of |
> (=OR) operators in my case? Or am I misinterpreting it?

Yes, this is the most direct way to choose the subset of a data.frame with
complete cases - used in lots of places that data.frame objects are used
(since 1992, White Book). Think data.frame objects, not GIS objects
(bands). This is about ontology, really.

Roger

>
> Kindest regards, Nikos
>
>

--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: [phone]; fax [phone]
e-mail: [hidden email]

_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Edzer Pebesma-2]

>>
>> So, this is the _correct and easy_ way instead of using a couple of |
>> (=OR) operators in my case? Or am I misinterpreting it?
>
> Yes, this is the most direct way to choose the subset of a data.frame
> with complete cases - used in lots of places that data.frame objects
> are used (since 1992, White Book). Think data.frame objects, not GIS
> objects (bands). This is about ontology, really.
Yes, or semantics at least.

Currently, a single gridded object in sp (SpatialGridDataFrame or
SpatialPixelsDataFrame) holds attributes in the @data slot, and this can
be a stack of bands that are spatially alligned, e.g. a set of bands
with different spectral characteristics. In addition bands (or variables
in the data.frame analogy) can be of different type, e.g. of type
logical or factor (think land use type).

The bands could also have a different time stamp, representing time
series of images, but the sp objects have no infrastructure to capture
the time stamp of a band.

--
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: [phone] http://ifgi.uni-muenster.de/
http://www.springer.com/978-0-387-78170-9 [hidden email]

_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Nikos Alexandris]

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


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

[...]

> An idea on how to apply this function:
>
> # load data
> # extract an NA mask using complete.cases
> # do something like prcomp()

> # assign back to original data
> x.df$pca_1[x.na_mask] <-  predict(p)[,1]
>
> This would of course be done the the @data slot on an sp object. Thanks to
> Roger for this slick approach.
>
> How about that!?

Very helpful Dylan. I am learning a lot, as usual. However, no-matter
what I try to do with my NA's I still get random NA's when I write back
to GRASS my data.

I tried to load the data, extract an NA_mask using complete.cases,
create "fake" columns and fill-in with the _untouched_ original data and
then write back to GRASS6. The result (or my process) is _noise_.

I'll try again tonight. If I am not going to make it I could try to
repeat this with spearfish data and see if something isn't ok with my
MODIS bands.

Cheers, Nikos

I don't understand what's the deal here :-(

_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Dylan Beaudette]

On Tuesday 14 April 2009, Nikos Alexandris wrote:

> Nikos:
> > > > > 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.
>
> Dylan:
> > > > 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.
>
> [...]
>
> > An idea on how to apply this function:
> >
> > # load data
> > # extract an NA mask using complete.cases
> > # do something like prcomp()
> >
> > # assign back to original data
> > x.df$pca_1[x.na_mask] <-  predict(p)[,1]
> >
> > This would of course be done the the @data slot on an sp object. Thanks
> > to Roger for this slick approach.
> >
> > How about that!?
>
> Very helpful Dylan. I am learning a lot, as usual. However, no-matter
> what I try to do with my NA's I still get random NA's when I write back
> to GRASS my data.
>
> I tried to load the data, extract an NA_mask using complete.cases,
> create "fake" columns and fill-in with the _untouched_ original data and
> then write back to GRASS6. The result (or my process) is _noise_.
>
> I'll try again tonight. If I am not going to make it I could try to
> repeat this with spearfish data and see if something isn't ok with my
> MODIS bands.
>
> Cheers, Nikos
>
> I don't understand what's the deal here :-(

I am pretty sure that readRAST6() respects the current region settings in
GRASS, however I am not sure how the 'plugin' option affects this. It could
be that your data are being read in at different resolutions, resulting in
different alignment of the 'bands' (columns) in the sp object. I would hope
that this is not possible, but I cannot think of anything else.

Can you reproduce this problem with another dataset? How about a subset of the
data?

Cheers,

Dylan



--
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341
_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Roger Bivand]

On Tue, 14 Apr 2009, Dylan Beaudette wrote:

> On Tuesday 14 April 2009, Nikos Alexandris wrote:
>> Nikos:
>>>>>> 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.
>>
>> Dylan:
>>>>> 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.
>>
>> [...]
>>
>>> An idea on how to apply this function:
>>>
>>> # load data
>>> # extract an NA mask using complete.cases
>>> # do something like prcomp()
>>>
>>> # assign back to original data
>>> x.df$pca_1[x.na_mask] <-  predict(p)[,1]
>>>
>>> This would of course be done the the @data slot on an sp object. Thanks
>>> to Roger for this slick approach.
>>>
>>> How about that!?
>>
>> Very helpful Dylan. I am learning a lot, as usual. However, no-matter
>> what I try to do with my NA's I still get random NA's when I write back
>> to GRASS my data.
>>
>> I tried to load the data, extract an NA_mask using complete.cases,
>> create "fake" columns and fill-in with the _untouched_ original data and
>> then write back to GRASS6. The result (or my process) is _noise_.
>>
>> I'll try again tonight. If I am not going to make it I could try to
>> repeat this with spearfish data and see if something isn't ok with my
>> MODIS bands.
>>
>> Cheers, Nikos
>>
>> I don't understand what's the deal here :-(
>
> I am pretty sure that readRAST6() respects the current region settings in
> GRASS, however I am not sure how the 'plugin' option affects this. It could
> be that your data are being read in at different resolutions, resulting in
> different alignment of the 'bands' (columns) in the sp object. I would hope
> that this is not possible, but I cannot think of anything else.

This looks possible, but if the data are all read in a single call to
readRAST6() is not, because plugin is set to FALSE even if present when
more than one raster is being transferred (and when the plugin raster
doesn't seem to be in the current region settings).

Without a reproducible example, we are not going to get anywhere, I'm
afraid. We must have a tarball of the location, and the exact R code being
used.

Roger

>
> Can you reproduce this problem with another dataset? How about a subset of the
> data?
>
> Cheers,
>
> Dylan
>
>
>
>

--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: [phone]; fax [phone]
e-mail: [hidden email]

_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Nikos Alexandris]

Roger:
> Without a reproducible example, we are not going to get anywhere, I'm
> afraid. We must have a tarball of the location, and the exact R code being
> used.

Dylan:
> > Can you reproduce this problem with another dataset? How about a subset of the
> > data?

Sure! After all it's MODIS data I am dealing with (...public domain) and
"my" location is not that big. I'll do my best tomorrow or the day after
(?). I'll let the list know :-)

Cheers, Nikos

_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Nikos Alexandris]

On Tue, 2009-04-14 at 23:36 +0200, Nikos Alexandris wrote:

> Roger:
> > Without a reproducible example, we are not going to get anywhere, I'm
> > afraid. We must have a tarball of the location, and the exact R code being
> > used.
>
> Dylan:
> > > Can you reproduce this problem with another dataset? How about a subset of the
> > > data?
>
> Sure! After all it's MODIS data I am dealing with (...public domain) and
> "my" location is not that big. I'll do my best tomorrow or the day after
> (?). I'll let the list know :-)
>
> Cheers, Nikos

Roger and Dylan,
apologies for not being on-time. I was _forced_ to spend time on other
things.

I will try to make a clean-check of my attempt one last time and if I
still can't make it I'll post info + grass-location with the data.

Last time I tried to repeat my test I received the same "random NULL's"
when using *only* one of the data-sets. I suspect that something might
be wrong with the data(?!). Not sure though (=I can't find _any_ NULL's
in the NULL's that are _produced_ when I write the data back from R to
GRASS).

Kindest regards, Nikos

_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats

[Dylan Beaudette-2]

On Fri, Apr 17, 2009 at 6:52 AM, Nikos Alexandris
<[hidden email]> wrote:

> On Tue, 2009-04-14 at 23:36 +0200, Nikos Alexandris wrote:
>> Roger:
>> > Without a reproducible example, we are not going to get anywhere, I'm
>> > afraid. We must have a tarball of the location, and the exact R code being
>> > used.
>>
>> Dylan:
>> > > Can you reproduce this problem with another dataset? How about a subset of the
>> > > data?
>>
>> Sure! After all it's MODIS data I am dealing with (...public domain) and
>> "my" location is not that big. I'll do my best tomorrow or the day after
>> (?). I'll let the list know :-)
>>
>> Cheers, Nikos
>
> Roger and Dylan,
> apologies for not being on-time. I was _forced_ to spend time on other
> things.
>
> I will try to make a clean-check of my attempt one last time and if I
> still can't make it I'll post info + grass-location with the data.
>
> Last time I tried to repeat my test I received the same "random NULL's"
> when using *only* one of the data-sets. I suspect that something might
> be wrong with the data(?!). Not sure though (=I can't find _any_ NULL's
> in the NULL's that are _produced_ when I write the data back from R to
> GRASS).
>
> Kindest regards, Nikos
>

Hi Nikos. I am not sure if this is still a possibility in GRASS, but I
seem to recall that some operations can result in values of 'Nan'...
which isn't properly encoded as NULL. Glynn would know for sure about
this. It could be you have NULL or NULL-like cells that you don't know
about...? r.report, and r.null may be useful in identification of NULL
cells.

Good luck!

Dylan
_______________________________________________
grass-stats mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-stats