Mask a map using statistical significance
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