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@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@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@data$pc1 <- NA
mod_pre_b267.raw@data$pc2 <- NA
mod_pre_b267.raw@data$pc3 <- NA

# fill-in with PCs
mod_pre_b267.raw@data$pc1[mod_pre_b267.na_mask] <- prcomp.mod_pre_b267
$x[,"PC1"]
mod_pre_b267.raw@data$pc2[mod_pre_b267.na_mask] <- prcomp.mod_pre_b267
$x[,"PC2"]
mod_pre_b267.raw@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

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

[Roger Bivand]

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.raw@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@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@data$pc1 <- NA
> mod_pre_b267.raw@data$pc2 <- NA
> mod_pre_b267.raw@data$pc3 <- NA
>
> # fill-in with PCs
> mod_pre_b267.raw@data$pc1[mod_pre_b267.na_mask] <- prcomp.mod_pre_b267
> $x[,"PC1"]
> mod_pre_b267.raw@data$pc2[mod_pre_b267.na_mask] <- prcomp.mod_pre_b267
> $x[,"PC2"]
> mod_pre_b267.raw@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: [phone]; fax [phone]
e-mail: [hidden email]



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

[Nikos Alexandris]

Thank you very much for your time Roger :-)

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


Roger:
> The NAs you show in the fake2 image are in band 7 (b7) and propagate
> correctly, b7 shown after import to R;

Right: _after_ import to R. They are just _999_ values prior importing
in R.


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

But why is "NODATA=" set to 999? It's a bug, right?

# reading the help page of readRAST6()
?readRAST6

---%<---
NODATA: by default NULL, in which case it is set to one less than
          'floor()'of the data values, otherwise an integer NODATA
          value (required to be integer by GRASS r.out.bin)
---%<---


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

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

[Nikos Alexandris]

Roger:
> I've committed a fix to sourceforge CVS, but setting
> > NODATA= manually is generally a good idea.


Nikos:
> But why is "NODATA=" set to 999? It's a bug, right?

Hmmm... one last thing: I see that it works properly only if some
out-of-range value is assigned to "NODATA=" (e.g. NODATA=-9999), right?
Setting NODATA=NULL is nonsense I guess.

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

[Roger Bivand]

On Wed, 22 Apr 2009, Nikos Alexandris wrote:

>
> Roger:
>> I've committed a fix to sourceforge CVS, but setting
>>> NODATA= manually is generally a good idea.
>
>
> Nikos:
>> But why is "NODATA=" set to 999? It's a bug, right?
>
> Hmmm... one last thing: I see that it works properly only if some
> out-of-range value is assigned to "NODATA=" (e.g. NODATA=-9999), right?
> Setting NODATA=NULL is nonsense I guess.
>

No, not nonsense. It is a default like any other. Working out what is in
range is not easy - it is based on r.info -r, but has to take edge cases
into account. The current "buglet" - the bug is the user not specifying
NODATA manually - was that in using parseGRASS in very recent releases,
one further layer of string to numeric conversion was needed to avoid the
value being set to 999 when the min and max values from r.info -r were not
"seen" as numeric.

There have been several similar 999 problems before, all can be avoided by
setting the NODATA= argument manually. For these bands, setting it
negative will always be OK, and setting manually also saves time, because
the heuristic doesn't have to be used. In general, NODATA= will simply
handle the NULL values correctly - the problem with any non-null value
given is that it may actually occur in the data in the current region.

Check out the CVS source from sourceforge - project r-spatial, module
spgrass6, it will work with your data.

Roger

>

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

Nikos:
> >> But why is "NODATA=" set to 999? It's a bug, right?
> > Hmmm... one last thing: I see that it works properly only if some
> > out-of-range value is assigned to "NODATA=" (e.g. NODATA=-9999), right?
> > Setting NODATA=NULL is nonsense I guess.


Roger:
> No, not nonsense. It is a default like any other. Working out what is in
> range is not easy - it is based on r.info -r, but has to take edge cases
> into account. The current "buglet" - the bug is the user not specifying
> NODATA manually -

:D

> was that in using parseGRASS in very recent releases,
> one further layer of string to numeric conversion was needed to avoid the
> value being set to 999 when the min and max values from r.info -r were not
> "seen" as numeric.

> There have been several similar 999 problems before, all can be avoided by
> setting the NODATA= argument manually. For these bands, setting it
> negative will always be OK, and setting manually also saves time, because
> the heuristic doesn't have to be used. In general, NODATA= will simply
> handle the NULL values correctly - the problem with any non-null value
> given is that it may actually occur in the data in the current region.

> Check out the CVS source from sourceforge - project r-spatial, module
> spgrass6, it will work with your data.

Thanks for your efforts Roger. I think I have to postpone for now
trying-out the latest source code (or maybe I'll try it out tonight).

I would like to note that I skimmed several times through "Applied
Spatial Data Analysis with R" [1] but I did not find any concrete
comment about the "NODATA=" option. Maybe I wasn't careful enough...

Kindest regards, Nikos
---
[1] Applied Spatial Data Analysis with R, Series: Use R
Bivand, Roger S., Pebesma, Edzer J., Gómez-Rubio, Virgilio

2008, XIV, 378 p., Softcover
ISBN: 978-0-387-78170-9


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

[hamish-2]

I posted to grass bug #73 what I should have posted here, (how to deal
with nodata is a common pain apparently)

http://trac.osgeo.org/grass/ticket/73#comment:47

* 999 is dangerous, it can easily be an elevation value. -9999 is bad
too, you just hit it less often (data doesn't have to be elevation in
meters). But even though you hit it less often doesn't mean it's still
not going to bite you one day, so it's a bug. but alas, probably for
r.out.arc we are stuck with what Arc allows for input values.


Hamish



     

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

[Roger Bivand]

On Wed, 22 Apr 2009, Nikos Alexandris wrote:

>
> Nikos:
>>>> But why is "NODATA=" set to 999? It's a bug, right?
>>> Hmmm... one last thing: I see that it works properly only if some
>>> out-of-range value is assigned to "NODATA=" (e.g. NODATA=-9999), right?
>>> Setting NODATA=NULL is nonsense I guess.
>
>
> Roger:
>> No, not nonsense. It is a default like any other. Working out what is in
>> range is not easy - it is based on r.info -r, but has to take edge cases
>> into account. The current "buglet" - the bug is the user not specifying
>> NODATA manually -
>
> :D
>
>> was that in using parseGRASS in very recent releases,
>> one further layer of string to numeric conversion was needed to avoid the
>> value being set to 999 when the min and max values from r.info -r were not
>> "seen" as numeric.
>
>> There have been several similar 999 problems before, all can be avoided by
>> setting the NODATA= argument manually. For these bands, setting it
>> negative will always be OK, and setting manually also saves time, because
>> the heuristic doesn't have to be used. In general, NODATA= will simply
>> handle the NULL values correctly - the problem with any non-null value
>> given is that it may actually occur in the data in the current region.
>
>> Check out the CVS source from sourceforge - project r-spatial, module
>> spgrass6, it will work with your data.
>
> Thanks for your efforts Roger. I think I have to postpone for now
> trying-out the latest source code (or maybe I'll try it out tonight).
>
> I would like to note that I skimmed several times through "Applied
> Spatial Data Analysis with R" [1] but I did not find any concrete
> comment about the "NODATA=" option. Maybe I wasn't careful enough...

Right, there isn't anything in the GRASS book either, where the interface
is described in more detail. In fact, the use of rgdal is quite recent, so
most of this implementation post-dates both sources.

Roger

>
> Kindest regards, Nikos
> ---
> [1] Applied Spatial Data Analysis with R, Series: Use R
> Bivand, Roger S., Pebesma, Edzer J., Gómez-Rubio, Virgilio
>
> 2008, XIV, 378 p., Softcover
> ISBN: 978-0-387-78170-9
>
>
>

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