[R-sig-Geo] Mask a map using statistical significance
Barry Rowlingson
b.rowlingson at lancaster.ac.uk
Wed May 4 18:02:17 CEST 2016
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
More information about the R-sig-Geo
mailing list