R - Contour plot from raster dataset with country borders overlaid

Question

I have a fairly simple and probably common task, plotting a raster dataset with countour lines and adding country borders together in one plot, however I did not find a solution anywhere. There are a a few hints available (such as this one), but no raster dataset is used there and I can't get it to work. The dataset I am using is actually in netcdf format and available here (15mb in size) and contains about 40 years of gridded precipitation data. Here is my line of code:

setwd("...netcdf Data/GPCP")
library("raster") 
library("maps")

nc_brick79_17 <- brick("precip.mon.mean.nc") # load in the ncdf data as a 
raster brick
newextent <- c(85, 125, -20, 20) # specify region of interest
SEA_brick <- crop(nc_brick79_17, newextent) # crop the region
day1 <- SEA_brick[[1]] # select very first day as example
colfunc<-colorRampPalette(c("white","lightblue","yellow","red","purple")) # colorscale for plotting

So it works of course when I just plot the raster data together with a map overlaid:

plot(day1, col=(colfunc(100)), interpolate=F, main="day1",legend.args=list(text='mm/hr', side=4,font=1, line=2.5, cex=1.1))
map("world", add=TRUE, lwd=0.5, interior = FALSE, col = "black")

We get this plot (Raster Plot with country borders added)

Now the code I use to generate the contour plot is the following:

filledContour(day1,zlim=c(0,20),color=colorRampPalette(c("white","lightblue","yellow","red","purple")), 
              xlab = "Longitude (°)", ylab = "Latitude (°)")
map("world", add=TRUE, lwd=0.5, interior = FALSE, col = "black") # add map overlay 

I end up with a plot where obviously the country borders do not align and are even covering the colorbar.

Contour plot with map overlay shifted

In this last part I am trying to add the country boundaries to the contour plot, but it does not work, even though it should I assume. The map is simply not there, no error though:

filledContour(day1, zlim=c(0,20), 
               color.palette = colorRampPalette(c("white","lightblue","yellow","red","purple")), 
               xlab = "Longitude (°)", ylab = "Latitude (°)", 
               xlim = c(90, 120), ylim = c(-20, 20), nlevels = 25, 
               plot.axes = {axis(1); axis(2);    
               map('world', xlim = c(90, 120), ylim = c(-20, 20), add = TRUE, lwd=0.5, col = "black")})  

From that line of code I get this plot.

Contour plot but no country borders added

What could I improve or is there any mistake somewhere? Thank you!

Answer 1

I chose to use ggplot here. I leave two maps for you. The first one is the one you created. This is a replication with ggplot. The second one is the one you could not produce. There are many things to explain. But I am afraid I do not have enough time to write all. But I left some comments in my code below. Please check this question to learn more about the second graphic. Finally, I'd like to give credit to hrbrmstr who wrote a great answer in the linked question.

library(maptools)
library(akima)
library(raster)
library(ggplot2)

# This is a data set from the maptools package
data(wrld_simpl)

# Create a data.frame object for ggplot. ggplot requires a data frame.

mymap <- fortify(wrld_simpl)

# This part is your code.

nc_brick79_17 <- brick("precip.mon.mean.nc")
newextent <- c(85, 125, -20, 20)
SEA_brick <- crop(nc_brick79_17, newextent)
day1 <- SEA_brick[[1]]


# Create a data frame with a raster object. This is a spatial class
# data frame, not a regular data frame. Then, convert it to a data frame.

spdf <- as(day1, "SpatialPixelsDataFrame")
mydf <- as.data.frame(spdf)
colnames(mydf) <- c("value", "x", "y")


# This part creates the first graphic that you drew. You draw a map.
# Then, you add tiles on it. Then, you add colors as you wish.
# Since we have a world map data set, we trim it at the end.

ggplot() +
geom_map(data = mymap, map = mymap, aes(x = long, y = lat, map_id = id), fill = "white", color = "black") +
geom_tile(data = mydf, aes(x = x, y = y, fill = value), alpha = 0.4) + 
scale_fill_gradientn(colors = c("white", "lightblue", "yellow", "red", "purple")) +
scale_x_continuous(limits = c(85, 125), expand = c(0, 0)) +
scale_y_continuous(limits = c( -20, 20), expand = c(0, 0)) +
coord_equal()

ggplot version of filled.contour()

# As I mentioned above, you want to study the linked question for this part.

mydf2 <- with(mydf, interp(x = x, 
                           y = y, 
                           z = value,
                           xo = seq(min(x), max(x), length = 400),
                           duplicate = "mean"))

gdat <- interp2xyz(mydf2, data.frame = TRUE)

# You need to draw countries as lines here. You gotta do that after you draw
# the contours. Otherwise, you will not see the map.

ggplot(data = gdat, aes(x = x, y = y, z = z)) + 
geom_tile(aes(fill = z)) + 
stat_contour(aes(fill = ..level..), geom = "polygon", binwidth = 0.007) + 
geom_contour(color = "white") +
geom_path(data = mymap, aes(x = long, y = lat, group = group), inherit.aes = FALSE) +
scale_x_continuous(limits = c(85, 125), expand = c(0, 0)) +
scale_y_continuous(limits = c(-20, 20), expand = c(0, 0)) +
scale_fill_gradientn(colors = c("white", "lightblue", "yellow", "red", "purple")) +
coord_equal() +
theme_bw()