Mask a map using statistical significance

Dear all

I am working on windows 7.

I have update my version of GRASS GIS 7.0.4 and R-3.2.5.

I try to work with my old tools and I meet a problem with function 
InitGrass. i just want to share this problem and if possible understand 
if the problem depend on my installation of grass and r or release of 
InitGrass.

Before I can use:
loc=initGRASS(gisBase =
Chemin_de_l_exe_de_Grass,home=Chemin_des_donnees_Grass,gisDbase=Chemin_des_donnees_Grass,location=Localisation,mapset=Jeu_de_donnees,override=T)
I have this problem:
Error in if (class(t0) != "try-error" && is.character(t0) && nchar(t0) > 
: missing value where TRUE/FALSE needed

I thonk it is in this osurce code:
https://github.com/cran/rgrass7/blob/master/R/xml1.R

I solve the problem uisng a false folder
loc=initGRASS(gisBase = 
Chemin_de_l_exe_de_Grass,home=Chemin_des_donnees_Grass,gisDbase=Chemin_des_donnees_Grass,addon_base=Chemin_de_l_exe_de_Grass,location=Localisation,mapset=Jeu_de_donnees,override=T)

Best regards

File: C:\Program Files\R\R-3.2.5\Cerema
#Chemin de l'exe de Grass
C:\GRASS GIS 7.0.4
#Monde Grass temporaire (1=oui, 0=non)
0
#Chemin des donnees Grass
C:\GRASSDATA
#Localisation
Cerema
#Jeu de donnees
test
#Projection (EPSG)
2154
#Chemin des routines R
C:/PROGRA~1/R/R-32~1.5/cerema


Code:
# DICARTO_00_LectureMondeGrass
# Version 2.0 06/10/2015
# Cerema
# Frederic Pons, Celine Trmal

DICARTO_00_LectureMondeGrass = function(chem_routine,Numero_de_calcul)
{
   # Lecture de l'initialisation et cr??ation du monde Grass
fichier_init=paste(chem_routine,"\\Init_Routine_Cerema",as.character(Numero_de_calcul),".txt",sep="")
   fid=file(fichier_init, open = "r")
   lignes <- readLines(fid)
   close(fid)

   Chemin_de_l_exe_de_Grass=as.character(lignes[2])
   Chemin_des_donnees_Grass=as.character(lignes[6])
   Localisation=as.character(lignes[8])
   Jeu_de_donnees=as.character(lignes[10])
   Projection=as.numeric(lignes[12])
   Monde_grass_temp=as.numeric(lignes[4])

   if (Monde_grass_temp==0) {
     loc=initGRASS(gisBase = 
Chemin_de_l_exe_de_Grass,home=Chemin_des_donnees_Grass,gisDbase=Chemin_des_donnees_Grass,addon_base=Chemin_de_l_exe_de_Grass,location=Localisation,mapset=Jeu_de_donnees,override=T)
     ## convert default coordinate system x,y to epsg code 3347 (can be 
altered later)
#execGRASS("g.mapset",flags=c("c"),parameters=list(mapset="PERMANENT"))
     #execGRASS("g.proj",flags=c("c"),epsg=as.integer(Projection))
#execGRASS("g.mapset",flags=c("c"),parameters=list(mapset=Jeu_de_donnees))
   } else {
     loc=initGRASS(Chemin_de_l_exe_de_Grass, 
home=tempdir(),addon_base=Chemin_de_l_exe_de_Grass,override=T)
   }

return(cbind(Chemin_de_l_exe_de_Grass,loc$GISDBASE,loc$LOCATION,loc$MAPSET,Projection))
}
*Fr?d?ric Pons *
*Expert hydraulique sur les inondations et al?as c?tiers
**DREC/Service Risques Inondations Littoraux et Hydraulique **- T?l.: 
(33)4 42 24 76 68 *
*Direction Territoriale M?diterran?e
*
Centre d??tudes et d?expertise sur les risques, l?environnement, la 
mobilit? et l?am?nagement
www.cerema.fr <http://www.cerema.fr>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20160504/a605343e/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: cerema.png
Type: image/png
Size: 6094 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20160504/a605343e/attachment.png>

Your example looks a bit rubbish because you have mostly isolated
pixels. Let's make an example with an area so we can see the shading.

r.pvalue = raster(outer(-24:25,-24:25, function(a,b){a^2+b^2})/120000)
range(r.pvalue[])

[1] 0.00000000 0.01041667

now convert your raster to polygons at the threshold you are
interested in. Let's imagine the area less than 0.001 (the central
circle here):

p_poly = rasterToPolygons(r.pvalue<0.001, dissolve=TRUE)

and overlay hatched polygons:

plot(r.pvalue) # or whatever your underlying data is
plot(p_poly[p_poly$layer==1,],density=5,add=TRUE, border=NA)

see help(polygon) for parameters to control the shading angle and
density. The "border=NA" above hides the outline which looks like
those maps you linked to.

Barry



On Wed, May 4, 2016 at 12:11 AM, Thiago V. dos Santos via R-sig-Geo

Dear all,

In climate studies, it is a common practice to perform some statistical test between two fields (maps), and then plot the resulting map using a significance mask. This masking is usually done by adding some kind of pattern (shading, crosshatching etc) on top of the actual color palette.

Examples can be seen here in this image https://www3.epa.gov/climatechange/images/science/GlobalPrecipMap-large.png and in the left images of this panel:
http://www.nature.com/nclimate/journal/vaop/ncurrent/images_article/nclimate2996-f3.jpg
In my case, I ran a statistical test for detecting trend on a time-series raster and I now have one raster with the trend for rainfall (in degree C per year) and one with the p-values associated to the test.

My data looks roughly like this:

require(raster)

## scratch raster objects and fill them with some random values
r.trend <- raster(nrows=50, ncols=50, xmn=-58, xmx=-48, ymn=-33, ymx=-22)
r.trend[] <- round(runif(50 * 50, min=0, max=3), digits=2)

r.pvalue <- raster(nrows=50, ncols=50, xmn=-58, xmx=-48, ymn=-33, ymx=-22)
r.pvalue[] <- round(runif(50 * 50, min=0.0001, max=0.1), digits=5)


What I would like to do is to plot r.trend, and on top of it plot r.pvalue (as a pattern) where r.pvalues < 0.01 (corresponding to a significance level of 99%).

Has anyone here ever tried to do a similar plot and can point me to any direction?

I usually use rasterVis and ggplot2 to plot maps, but I would be open to try some other package and even other SIG software other than R.

Many thanks in advance, -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Dear all

I am working on windows 7.

I have update my version of GRASS GIS 7.0.4 and R-3.2.5.

I try to work with my old tools and I meet a problem with function 
InitGrass. i just want to share this problem and if possible understand 
if the problem depend on my installation of grass and r or release of 
InitGrass.

You have not said which version of rgrass7 you are using.

I'm afraid this is very unclear. You are almost certainly using a very 
suboptimal approach. Have you thought of starting GRASS properly once, 
and doing whatever your function does in new mapsets? Or are the locations 
different?

Please move this to the grass-stats list.

After an error, do run traceback(). And establish first that initGRASS() 
works run in the R session itself rather than in a function (I doubt that 
this is easy to debug).

Roger

PS: Please note that the link to github is *not* the source, it is an 
unauthorised copy, and as maintainer I refuse to consider any references 
to it. The released source is on CRAN, and the development trunk may be 
browsed on R-Forge. Just "thinking" is no use, if you run under debug(), 
give the line number of specific issues.

Before I can use:
loc=initGRASS(gisBase =
Chemin_de_l_exe_de_Grass,home=Chemin_des_donnees_Grass,gisDbase=Chemin_des_donnees_Grass,location=Localisation,mapset=Jeu_de_donnees,override=T)
I have this problem:
Error in if (class(t0) != "try-error" && is.character(t0) && nchar(t0) > 
:  missing value where TRUE/FALSE needed

I thonk it is in this osurce code:
https://github.com/cran/rgrass7/blob/master/R/xml1.R

I solve the problem uisng a false folder
loc=initGRASS(gisBase = 
Chemin_de_l_exe_de_Grass,home=Chemin_des_donnees_Grass,gisDbase=Chemin_des_donnees_Grass,addon_base=Chemin_de_l_exe_de_Grass,location=Localisation,mapset=Jeu_de_donnees,override=T)

Best regards

File: C:\Program Files\R\R-3.2.5\Cerema
#Chemin de l'exe de Grass
C:\GRASS GIS 7.0.4
#Monde Grass temporaire (1=oui, 0=non)
0
#Chemin des donnees Grass
C:\GRASSDATA
#Localisation
Cerema
#Jeu de donnees
test
#Projection (EPSG)
2154
#Chemin des routines R
C:/PROGRA~1/R/R-32~1.5/cerema


Code:
#  DICARTO_00_LectureMondeGrass
#  Version 2.0 06/10/2015
#  Cerema
#  Frederic Pons, Celine Trmal

DICARTO_00_LectureMondeGrass = function(chem_routine,Numero_de_calcul)
{
 # Lecture de l'initialisation et cr??ation du monde Grass
fichier_init=paste(chem_routine,"\\Init_Routine_Cerema",as.character(Numero_de_calcul),".txt",sep="")
  fid=file(fichier_init, open = "r")
  lignes <- readLines(fid)
  close(fid)

  Chemin_de_l_exe_de_Grass=as.character(lignes[2])
  Chemin_des_donnees_Grass=as.character(lignes[6])
  Localisation=as.character(lignes[8])
  Jeu_de_donnees=as.character(lignes[10])
  Projection=as.numeric(lignes[12])
  Monde_grass_temp=as.numeric(lignes[4])

  if (Monde_grass_temp==0) {
   loc=initGRASS(gisBase = 
Chemin_de_l_exe_de_Grass,home=Chemin_des_donnees_Grass,gisDbase=Chemin_des_donnees_Grass,addon_base=Chemin_de_l_exe_de_Grass,location=Localisation,mapset=Jeu_de_donnees,override=T)
   ## convert default coordinate system x,y to epsg code 3347 (can be 
altered later)
#execGRASS("g.mapset",flags=c("c"),parameters=list(mapset="PERMANENT"))
   #execGRASS("g.proj",flags=c("c"),epsg=as.integer(Projection))
#execGRASS("g.mapset",flags=c("c"),parameters=list(mapset=Jeu_de_donnees))
  } else {
    loc=initGRASS(Chemin_de_l_exe_de_Grass, 
home=tempdir(),addon_base=Chemin_de_l_exe_de_Grass,override=T)
  }

return(cbind(Chemin_de_l_exe_de_Grass,loc$GISDBASE,loc$LOCATION,loc$MAPSET,Projection))
}
*Fr?d?ric Pons *
*Expert hydraulique sur les inondations et al?as c?tiers
**DREC/Service Risques Inondations Littoraux et Hydraulique **- T?l.: (33)4 
42 24 76 68 *
*Direction Territoriale M?diterran?e
*
Centre d??tudes et d?expertise sur les risques, l?environnement, la mobilit? 
et l?am?nagement
www.cerema.fr <http://www.cerema.fr>

Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
http://depsy.org/person/434412

Thanks Barry and Arnaud (off-list) for the valuable hints.
I ended up managing to use rasterVis' levelplot to produce the map I was looking for.
The key step was to convert the raster I wanted to use as a mask to SpatialPoints, and then pass it as an additional layer to levelplot.
The final map, produced using my actual data, can be visualized here:?https://dl.dropboxusercontent.com/u/27700634/rainfall_trend.png.?Stippling indicates regions exceeding the 95% statistical significance.
?
| ? |
| ? |  | ? | ? | ? | ? | ? |
|  |
|  |
| View on dl.dropboxuserconten... | Preview by Yahoo |
|  |
| ? |


Greetings,
?-- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota

Your example looks a bit rubbish because you have mostly isolated
pixels. Let's make an example with an area so we can see the shading.

r.pvalue = raster(outer(-24:25,-24:25, function(a,b){a^2+b^2})/120000)
range(r.pvalue[])

[1] 0.00000000 0.01041667

now convert your raster to polygons at the threshold you are
interested in. Let's imagine the area less than 0.001 (the central
circle here):

p_poly = rasterToPolygons(r.pvalue<0.001, dissolve=TRUE)

and overlay hatched polygons:

plot(r.pvalue) # or whatever your underlying data is
plot(p_poly[p_poly$layer==1,],density=5,add=TRUE, border=NA)

see help(polygon) for parameters to control the shading angle and
density. The "border=NA" above hides the outline which looks like
those maps you linked to.

Barry




On Wed, May 4, 2016 at 12:11 AM, Thiago V. dos Santos via R-sig-Geo

Dear all,

In climate studies, it is a common practice to perform some statistical test between two fields (maps), and then plot the resulting map using a significance mask. This masking is usually done by adding some kind of pattern (shading, crosshatching etc) on top of the actual color palette.

Examples can be seen here in this image https://www3.epa.gov/climatechange/images/science/GlobalPrecipMap-large.png and in the left images of this panel:
http://www.nature.com/nclimate/journal/vaop/ncurrent/images_article/nclimate2996-f3.jpg
In my case, I ran a statistical test for detecting trend on a time-series raster and I now have one raster with the trend for rainfall (in degree C per year) and one with the p-values associated to the test.

My data looks roughly like this:

require(raster)

## scratch raster objects and fill them with some random values
r.trend <- raster(nrows=50, ncols=50, xmn=-58, xmx=-48, ymn=-33, ymx=-22)
r.trend[] <- round(runif(50 * 50, min=0, max=3), digits=2)

r.pvalue <- raster(nrows=50, ncols=50, xmn=-58, xmx=-48, ymn=-33, ymx=-22)
r.pvalue[] <- round(runif(50 * 50, min=0.0001, max=0.1), digits=5)


What I would like to do is to plot r.trend, and on top of it plot r.pvalue (as a pattern) where r.pvalues < 0.01 (corresponding to a significance level of 99%).

Has anyone here ever tried to do a similar plot and can point me to any direction?

I usually use rasterVis and ggplot2 to plot maps, but I would be open to try some other package and even other SIG software other than R.

Many thanks in advance, -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

I will appreciate instructions to reproduce this map...

Thanks...

Ruben FC

El 05/05/2016 00:06, "Thiago V. dos Santos via R-sig-Geo" <
r-sig-geo at r-project.org> escribi?: Thanks Barry and Arnaud (off-list) for
the valuable hints. I ended up managing to use rasterVis' levelplot to
produce the map I was looking for. The key step was to convert the raster I
wanted to use as a mask to SpatialPoints, and then pass it as an additional
layer to levelplot. The final map, produced using my actual data, can be
visualized here: https://dl.
dropboxusercontent.com/u/27700634/rainfall_trend.png. Stippling indicates
regions exceeding the 95% statistical significance.

Ruben, 

This is a tentative reproducible example using "real" spatial data, rather than the synthetic data I proposed in the original question.

(I should have thought of this example earlier, in order to spatially please Barry lol)

require(raster)
require(rasterVis)

# Scratch raster objects
data(volcano)
r1 <- raster(volcano)

over <- ifelse(volcano >=160 & volcano <=180, 1, NA) # This is the "mask" raster
r2 <- raster(over)

# And this is the key step:
# To convert the "mask" raster to spatial points
r.mask <- rasterToPoints(r2, spatial=TRUE)

# Plot
levelplot(r1, margin=F) +
layer(sp.points(r.mask, pch=20, cex=0.3, alpha=0.8))


It looks a little better now. To control the parameters of the points (such as color, size, type etc), documentation is available in ?sp.points . 
 Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota

I will appreciate instructions to reproduce this map...
Thanks...
Ruben FC

El 05/05/2016 00:06, "Thiago V. dos Santos via R-sig-Geo" <r-sig-geo at r-project.org> escribi?: Thanks Barry and Arnaud (off-list) for the valuable hints. I ended up managing to use rasterVis' levelplot to produce the map I was looking for. The key step was to convert the raster I wanted to use as a mask to SpatialPoints, and then pass it as an additional layer to levelplot. The final map, produced using my actual data, can be visualized here: https://dl. dropboxusercontent.com/u/27700634/rainfall_trend.png. Stippling indicates regions exceeding the 95% statistical significance.