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:


# 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(poly, add=T)

# Mask and plot
rr <- mask(r, poly)

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[[1]]@Polygons[[1]]@coords
newHole <- Polygon(coordsHole,hole=TRUE)
listPol <- polyX@polygons[[1]]@Polygons
listPol[[length(listPol)+1]] <- newHole
punch <- Polygons(listPol,polyX@polygons[[1]]@ID)
new <- SpatialPolygons(list(punch),proj4string=poly@proj4string)

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.