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

[Nikos Alexandris]

Roger, Dylan and grass-stats-listers,

  I can't "solve" the random-NULL production of writeRAST6. I tried
again with only 3 bands which are identical with respect to their
extent/resolution/etc. I am not sending more data (e.g. as I promised to
send 6 modis bands) since the "problem" occurs with simpler dataset(s).

  The grass-location with the specific MODIS bands is the one I have
already uploaded in gregis [1].

(
I also repeated the test with landsat7 images from the spearfish
location but all works fine there!
)

  I am trying to present the problem as simple as possible (commands
follow after step-by-step description):

Step-by-step:
1. launch R from within a GRASS-location
2. load library spgrass6 and projection information (gmeta6)
3. load _three_ modis bands in one object
4. create an NA-mask using the function complete.cases()
5. extract "values" only form the initial 3-modis-bands object
6. perform pca on "values" based on svd [that is with function prcomp()]
7. create new columns in the initial 3-modis-bands object
8. fill-in with the PC's produced by prcomp()
9. write back to grass with writeRAST6

Result:
Viewing the produced raster maps there are many "random" NULLs where all
original raster maps are filled with some integer value.

I guess those NULLs are not "random" at all - but I have no idea what
the reason for this behavior is.

(
Check the "Value Tool", bottom left on the screen-shots [2][3][4]. Note
that the mouse-pointer is a bit "shifted" -- don't know why, but it
should be over the "black pixels"!
)

# # # # # # # # # # #
# using: grass65, R2.8.1, latest spgrass6 installed today, GDAL 1.6
# on: Ubuntu Intrepid 64-bit
# launch grass in location/mapset "peloponnese_hgrs87/PERMAMENT"
# code in R

# load spgrass6 and projection information
library(spgrass6) ; G <- gmeta6()

# load modis bands
mod_pre_b267.raw <- readRAST6 (c('mod_b2', 'mod_b6', 'mod_b7'))
# str(mod_pre_b267.raw)
# summary(mod_pre_b267.raw)

# create NA mask using complete.cases()
mod_pre_b267.na_mask <- complete.cases(mod_pre_b267.raw at data)
# str(mod_pre_b267.na_mask)
# summary(mod_pre_b267.na_mask)

# get values based on na_mask
mod_pre_b267 <- mod_pre_b267.raw at data[mod_pre_b267.na_mask, ]

# perform pca (svd)
prcomp.mod_pre_b267 <- prcomp(mod_pre_b267, scale=TRUE)
summary(prcomp.mod_pre_b267)
prcomp.mod_pre_b267$rotation

# create new columns in initial object
mod_pre_b267.raw at data$pc1 <- NA
mod_pre_b267.raw at data$pc2 <- NA
mod_pre_b267.raw at data$pc3 <- NA

# fill-in with PCs
mod_pre_b267.raw at data$pc1[mod_pre_b267.na_mask] <- prcomp.mod_pre_b267
$x[,"PC1"]
mod_pre_b267.raw at data$pc2[mod_pre_b267.na_mask] <- prcomp.mod_pre_b267
$x[,"PC2"]
mod_pre_b267.raw at data$pc3[mod_pre_b267.na_mask] <- prcomp.mod_pre_b267
$x[,"PC3"]

# write back to grass
writeRAST6(mod_pre_b267.raw, zcol=4, vname="fake1", overwrite=TRUE)
writeRAST6(mod_pre_b267.raw, zcol=5, vname="fake2", overwrite=TRUE)
writeRAST6(mod_pre_b267.raw, zcol=6, vname="fake3", overwrite=TRUE)

# check the "random" NULL's, e.g. with qgis
# # # # # # # # # # #


Hopefully you will have the time to have a look at it. Kindest regards,
Nikos
---

References
[1] grass location containint three MODIS surface reflectance bands:
https://git.berlios.de/cgi-bin/gitweb.cgi?p=gregis;a=blob_plain;f=peloponnesos/modis_peloponnese_postfire07.zip;hb=peloponnese

[2] http://imageshack.gr/view.php?file=owj2xu833u6678c3zvch.png
[3] http://imageshack.gr/view.php?file=n65vfvoo07t3njywfjjo.png
[4] http://imageshack.gr/view.php?file=sb9l315wazoeue0rgrj9.png