[R-sig-Geo] Mask a map using statistical significance

[Thiago V. dos Santos]

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


On Wednesday, May 4, 2016 11:02 AM, Barry Rowlingson <b.rowlingson at lancaster.ac.uk> wrote:



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
<r-sig-geo at r-project.org> wrote:
> 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
	[[alternative HTML version deleted]]