On Tue, 21 Apr 2009, Nikos Alexandris wrote:
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 NAs you show in the fake2 image are in band 7 (b7) and propagate correctly, b7 shown after import to R; the grass_nas image is from within GRASS showing the vector returned by complete.cases().
They are getting there because the NODATA= argument is not given, and is being set to 999. I've committed a fix to sourceforge CVS, but setting NODATA= manually is generally a good idea.
Roger
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....@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....@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....@data$pc1 <- NA
mod_pre_b267....@data$pc2 <- NA
mod_pre_b267....@data$pc3 <- NA
# fill-in with PCs
mod_pre_b267....@data$pc1[mod_pre_b267.na_mask] <- prcomp.mod_pre_b267
$x[,"PC1"]
mod_pre_b267....@data$pc2[mod_pre_b267.na_mask] <- prcomp.mod_pre_b267
$x[,"PC2"]
mod_pre_b267....@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
-- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: roger.biv...@nhh.no
<<attachment: grass_nas.png>>
<<attachment: b7nas.png>>
_______________________________________________ grass-stats mailing list grass-stats@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-stats
