I am trying to produce a map (using R) that uses a color-scale to display raster data. I want to display the raster data only underneath a particular polygon. My first thought was to use
raster::mask to generate a new raster that is cropped to the polygon. The problem is that, because of the low resolution of the raster, I end up with jagged rather than smooth edges. Here's an example:
library(sp) library(raster) # Create Raster r <- raster(ncol=40, nrow=20) r <- rnorm(n=ncell(r)) # Create Polygon a <- seq(-150,150,1) b <- 50*a^2/(150^2) df <- data.frame(a=a,b=b) srl <- Polygon(df, hole=as.logical(NA)) Srl <- Polygons(list(srl), 1) poly <- SpatialPolygons(list(Srl), pO = 1:length(Srl), proj4string=CRS(as.character(NA))) # Plot raster and polygon plot(r) plot(poly, add=T) # Mask and plot rr <- mask(r, poly) plot(rr)
What I want is something that looks like this instead (this code is clumsy, but it shows what I want; basically it works by plotting a big white polygon on top of the raster with a hole over the area I want visible):
# The previous code snippet must be executed prior to this one e1 <- c(-180,180,180,-180) e2 <- c(90,90,-90,-90) df2 <- data.frame(e1=e1,e2=e2) srl2 <- Polygon(df2, hole=F) Srl2 <- Polygons(list(srl2), 2) poly2 <- SpatialPolygons(list(Srl2)) srl3 <- Polygon(df) Srl3 <- Polygons(list(srl3), 1) poly3 <- SpatialPolygons(list(Srl2,Srl3)) hole <- poly polyX <- poly2 coordsHole <- hole@polygons[]@Polygons[]@coords newHole <- Polygon(coordsHole,hole=TRUE) listPol <- polyX@polygons[]@Polygons listPol[[length(listPol)+1]] <- newHole punch <- Polygons(listPol,polyX@polygons[]@ID) new <- SpatialPolygons(list(punch),proj4string=poly@proj4string) plot(r) plot(new, col="white", border="NA", add=T)
In my case, I can't get away with this dumb hack, because the raster needs to be plotted on top of an existing plot. So how can I plot the raster, where pixels along the polygon edge are partially plotted, so only the parts inside the polygon show?
I don't want to just increase the raster resolution with
raster::disaggregate for a number of reasons. The raster is already big to begin with, and it could get pretty unwieldy. Also, for visualization purposes I want to interpolate the plotted raster (i.e. using
raster::plot(r, interpolate='bilinear'), and I'm pretty sure
disaggregate would mess with that.
Note that if I use ggplot2 to plot the raster with
geom_raster, it handles the edges exactly how I want. Unfortunately, I need to plot the raster as part of a larger plot for which I cannot use ggplot2.