Overlay raster layer on map in ggplot2 in R?

Question

I am trying to overlay a raster layer onto a map in ggplot. The raster layer contains likelihood surfaces for each time point from a satellite tag. I also want to set cumulative probabilities(95%, 75%, 50%) on the raster layer.

I have figured out how to show the raster layer on the ggplot map, but the coordinates are not aligned with one another. I tried making each have the same projection but it does not seem to be working... I want them both to fit the boundaries of my model (xmin = 149, xmax = 154, ymin = -14, ymax = -8.75

Attached is my r code and the figure result:

#load data
ncname <- "152724-13-GPE3"
ncfname <- paste(ncname, ".nc", sep = "")
ncin <- nc_open(ncfname)

StackedObject<-stack("152724-13-GPE3.nc", varname = "monthly_residency_distributions")
MergedObject<-overlay(StackedObject,fun=mean ) 
MergedObject[is.na(MergedObject)]<-0 
Boundaries<-extent(c(149, 154, -14, -8.75)) 
ExtendedObject<-extend(MergedObject, Boundaries) 
Raster.big<-raster(ncol=1200,nrow=900,ext=Boundaries) 
Raster.HR<-resample(x=ExtendedObject, y=Raster.big, method="bilinear") 
Raster.HR@data@values<- Raster.HR@data@values/sum(Raster.HR@data@values) 
RasterVals<-sort(Raster.HR@data@values)
Raster.breaks <- c(RasterVals[max(which(cumsum(RasterVals)<= 0.05 ))], RasterVals[max(which(cumsum(RasterVals)<= 0.25 ))], RasterVals[max(which(cumsum(RasterVals)<= 0.50 ))], 1)
Raster.cols<-colorRampPalette(c("yellow","orange","red"))
RasterCols<- c(Raster.cols(3))

#Create Map
shape2 <- readOGR(dsn = "/Users/shannonmurphy/Desktop/PNG_adm/PNG_adm1.shp", layer = "PNG_adm1")
map<- crop(shape2, extent(149, 154, -14, -8.75))
projection(map)<- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

p <- ggplot() + geom_polygon(data = map, aes(x = long, y = lat, group = group), color = "black", size = 0.25) + coord_map()

projection(Raster.HR)<- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

#plot raster and ggplot 

par(mfrow=c(1,1))
plot(p)
par(mfrow=c(1,1), new = TRUE)
plot(Raster.HR, col=RasterCols, breaks=Raster.breaks, legend = NULL, bbox(map)) 

Please let me know if there is another package/line of code I should be using to do this! Appreciate any help

Answer 1

Ok I understand. You want to plot multiple raster layers on the ggplot or you want that the raster object is over your background polygon object. The problem with rasterVis::gplot is that it directly plot the raster and does not allow to add another one under or over. You remind me that I already had this need and modified function gplot to retrieve the data as a tibble so that you can then play with it as much as you want with dplyr and then ggplot2. Thanks for the reminder, I added it in my current github library for later use!
Let's use a reproducible example to show this function:

Create datasets

• Create a map of the world as a Raster to be use as background Raster map
• Create a raster of data, here a distance from a point (limited to a maximum distance)

The code:

library(raster)

# Get world map
library(maptools)
data(wrld_simpl)
# Transform World as raster
r <- raster(wrld_simpl, res = 1)
wrld_r <- rasterize(wrld_simpl, r)

# Lets create a raster of data
pt1 <- matrix(c(100,0), ncol = 2)
dist1 <- distanceFromPoints(r, pt1)
values(dist1)[values(dist1) > 5e6] <- NA
plot(dist1)

# Plot both
plot(wrld_r, col = "grey")
plot(dist1, add = TRUE)

Function to extract (part of) raster values and transform as a tibble

#' Transform raster as data.frame to be later used with ggplot
#' Modified from rasterVis::gplot
#'
#' @param x A Raster* object
#' @param maxpixels Maximum number of pixels to use
#'
#' @details rasterVis::gplot is nice to plot a raster in a ggplot but
#' if you want to plot different rasters on the same plot, you are stuck.
#' If you want to add other information or transform your raster as a
#' category raster, you can not do it. With `SDMSelect::gplot_data`, you retrieve your
#' raster as a tibble that can be modified as wanted using `dplyr` and
#' then plot in `ggplot` using `geom_tile`.
#' If Raster has levels, they will be joined to the final tibble.
#'
#' @export

gplot_data <- function(x, maxpixels = 50000)  {
  x <- raster::sampleRegular(x, maxpixels, asRaster = TRUE)
  coords <- raster::xyFromCell(x, seq_len(raster::ncell(x)))
  ## Extract values
  dat <- utils::stack(as.data.frame(raster::getValues(x))) 
  names(dat) <- c('value', 'variable')

  dat <- dplyr::as.tbl(data.frame(coords, dat))

  if (!is.null(levels(x))) {
    dat <- dplyr::left_join(dat, levels(x)[[1]], 
                            by = c("value" = "ID"))
  }
  dat
}

Plot multiple rasters in ggplot

You can use gplot_data to transform any raster as a tibble. You are then able to add any modification using dplyr and plot on ggplot with geom_tile. The interesting point is that you can use geom_tile as many time as you want with different raster data, provided that fill option is comparable. Otherwise, you can use the trick below to remove NA values in the background raster map and use a unique fill colour.

# With gplot_data
library(ggplot2)

# Transform rasters as data frame 
gplot_wrld_r <- gplot_data(wrld_r)
gplot_dist1 <- gplot_data(dist1)

# To define a unique fill colour for the world map,
# you need to remove NA values in gplot_wrld_r which 
# can be done with dplyr::filter
ggplot() +
  geom_tile(data = dplyr::filter(gplot_wrld_r, !is.na(value)), 
            aes(x = x, y = y), fill = "grey20") +
  geom_tile(data = gplot_dist1, 
            aes(x = x, y = y, fill = value)) +
  scale_fill_gradient("Distance",
                      low = 'yellow', high = 'blue',
                      na.value = NA) +
  coord_quickmap()

Plot raster over polygons

Of course, with a background map as a polygon object, this trick also let you add your raster over it:

wrld_simpl_sf <- sf::st_as_sf(wrld_simpl)

ggplot() +
  geom_sf(data = wrld_simpl_sf, fill = "grey20",
          colour = "white", size = 0.2) +
  geom_tile(data = gplot_dist1, 
            aes(x = x, y = y, fill = value)) +
  scale_fill_gradient("Distance",
                      low = 'yellow', high = 'blue',
                      na.value = NA)

EDIT: gplot_data is now in this simple R package: https://github.com/statnmap/cartomisc