How to replicate ggplot2::geom_raster in base graphics?

Question

tl;dr

If you run both code snippets below (note: code was previously broken, but now is fixed and has been checked on multiple computers), you will see two plots of some raster data. One uses ggplot2 and produces a smoothed image with high-resolution coastlines that are somehow inherited from the polygons that I used to mask the raster.

Without using ggplot2, we can get the smoothed image using raster::plot(... , interpolate='bilinear') or rasterImage(interpolate = TRUE). But the coastlines in the map don't look good (they are jagged/pixellated, according to the resolution of the raster). You may have to zoom in quite a bit to see this. I have two questions:

1. How can I plot the raster in base graphics, but only where it overlaps a polygon whose edge has higher resolution than the raster (yielding something like the ggplot output)? I ask this question in more detail here, with an example of the output I want.
2. How does ggplot2 understand where the polygon edges are for correct plotting? I cannot see where/how I tell ggplot what the masking polygon is!

Gory details

I am plotting a raster in R. I initially did the plotting with ggplot2, using geom_raster(), but now I need to switch back to base R because I ultimately need to display multiple rasters with different color-scales on a single plot (this is a map, and color-scales will eventually differ by continent).

The code below will reproduce the plots (including downloading the necessary files). On my computer, the full code (both snippets below) takes about 1.5 minutes to run--most of this is drawing/rendering the plots.

library("rgdal")
library("sp")
library("raster")
library("ggplot2")
    
# Downloading and loading raster
download.file("https://dl.dropboxusercontent.com/u/49481174/Stack.zip", destfile="Stack.zip")
unzip("Stack.zip")
msrP2 <- raster("Stack/msrP2.grd")

plot(msrP2, interpolate='bilinear')

Note the pixellation along coastlines (you may need to zoom in to see this). This is expected based on the resolution of the raster.

But watch what happens in ggplot2!

warning, rendering the plot requires a decent chunk of memory on my machine

# Convert to data.frame for ggplot2
msr.p <- rasterToPoints(msrP2)
mdf <- data.frame(msr.p)
colnames(mdf) <- c("Longitude", "Latitude", "RichnessProp")

# plot using geom_raster
b_map <- ggplot(data=mdf, aes_string(y="Latitude", x="Longitude")) +
  geom_raster(aes(fill=RichnessProp)) +
  theme_bw() +
  coord_equal() +
  scale_fill_gradientn("DD", colours=c("gold1", "coral", "navy"), na.value="black") +
  theme(axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "right",
    legend.key = element_blank(),
    axis.ticks=element_blank()) +
  scale_y_continuous(limits=c(-5133051,6324167)) + 
  scale_x_continuous(limits=c(-20100000,20100000))

b_map

Now the coastlines are reproduced in high fidelity! Ultimately, the reason for this must have something to do with the fact that this raster was originally produced by masking a larger raster using high-resolution shapefiles of continents and lakes. So to re-state the questions, in order of importance:

Q1. How can I create this high-resolution plot in base R?

Q2. How in the heck does ggplot "know" where the coastline is based on the data.frame that I am passing to it??

Answer 1

Since you say that ggplot takes a huge amount of memory while plot doesn't, I'd think that the pixellation comes from the fact that raster::plot by default actually takes in only a subset of pixels to reduce memory footprint, depending on argument maxpixels:

maxpixels integer > 0. Maximum number of cells to use for the plot. If maxpixels < ncell(x), sampleRegular is used before plotting. If gridded=TRUE maxpixels may be ignored to get a larger sample

I'd try something like:

plot(msrP2, interpolate='bilinear', maxpixels=50000000)

to see if it solves your issue.

PS: Concernign Q2: ggplot2 can not know anything about the "coastline". If I'm right, the problem you are seeing is just related to the fact that raster::plot automatically "reduces" resolution to save you memory.