Plotting polygons with rasters in base R, ggplot2 or levelplot

Question

I am trying to plot in R a raster layer with lines/polygon objects in R and each time I fail miserably with errors. I tried to do this in base R, ggplot2 and using levelplot but can't get the right result.

Source data can be found here.

What I need to do in the plot (all in one plot) is to: 1) zoom in a certain area defined as NIG. T 2) Display raster r values on a scale with cuts intervals. 3) Plot the country boundaries(shpAfr in base R and ggplot2 or world.outlines.spin levelplot). 4) Finally, include shpWater polygon layer (with col="blue" fill and contours).

library(raster)
library(maptools)
library(rasterVis)
library(viridis)
library(sf)
library(rgdal)
library(ggplot2)

r <- raster("raster_example.tif")
crs(r) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +to wgs84=0,0,0"
NIG <- c(2,14.5,4,14)
Reg_name <- "Nigeria"

shpAfr <- readOGR(dsn="Africa.shp")
proj4string(shpAfr) # describes data’s current coordinate reference system
#st_read(system.file("shape/nc.shp", package="sf"))

# Import water polygon
shpWater <- readOGR(dsn="waterbodies_africa.shp")
shpWater.count <- nrow(shpWater@data)
shpWater$id <- 1:shpWater.count
shpWater.fort <- fortify(shpWater, region='id')

# Import Africa admin map
shpAfr <- readOGR(dsn="Africa.shp")
shpAfr.count <- nrow(shpAfr@data)
shpAfr$id <- 1:shpAfr.count
shpAfr.fort <- fortify(shpAfr, region='id')

# Set colour intervals for plotting:
cuts=seq(0,1,0.1) #set breaks

Trying in base R, my problem is I can get the water shape fill in the right colour (fill and contour should be blue). If I try to plot both wrld_simpl and shpWater as polygon() I get into even bigger troubles.

plot(r, xlim = NIG[1:2], ylim = NIG[3:4],
     breaks=cuts, col = rev(plasma(11)))
lines(wrld_simpl,lwd = 1.5)
lines(shpWater, col="blue") # works but cannot fill the polygon
polygon(shpWater, col = "blue", border = "blue") # getting error here
Error in as.double(y) : 
  cannot coerce type 'S4' to vector of type 'double'

Ok, so now I try ggplot2, but I can't find a way to include a raster here without getting an error.

lon <- seq(r@extent@xmin,r@extent@xmax,
           (r@extent@xmax-r@extent@xmin)/r@ncols)
lat <- seq(r@extent@ymin,r@extent@ymax,
           (r@extent@ymax-r@extent@ymin)/r@nrows)

Plot1 <- ggplot()+
  geom_polygon(aes(x = long, y = lat, group=id),
               data = shpAfr.fort, color ="grey27", fill ="grey",
               alpha = .4, size = .2)+
  geom_raster(data = test, aes(fill=values))+ ## here it goes bad
  #geom_tile(data=test_df, aes(x=x, y=y, fill=value), alpha=0.8) + 
  #scale_fill_viridis() +
  geom_polygon(aes(x = long, y = lat, group=id),
               data = shpWater.fort, color ="lightskyblue2", fill ="lightskyblue2",
               size = .2)+coord_equal()+
  theme_minimal()+
  coord_map(xlim = Region[[3]][1:2],ylim = Region[[3]][3:4])

plot(Plot1)

Finally, I tried the levelplot and AGAIN failed.

mapTheme <- rasterTheme(region = rev(brewer.pal(10, "RdBu")))

# Get world outlines:
world.outlines <- map("world", plot=FALSE)
world.outlines.sp <- map2SpatialLines(world.outlines, proj4string = CRS("+proj=longlat"))

# Plot raster and polygon:
Plot2 <- levelplot(r,par.settings = mapTheme,pretty=TRUE,margin = F,
          xlim = NIG[1:2],ylim = NIG[3:4],
          col.regions=colorRampPalette(c("light blue","blue", "red")),
          main=paste0("test")) + layer(sp.lines(world.outlines.sp, col = "black", lwd = 0.5))
plot(Plot2 + layer(sp.lines(world.outlines.sp, col = "black", lwd = 0.5))
#Error: Attempted to create layer with no stat.

My results so far: 1) first image does not have the polygons filled with blue 2) second image has clearly world outlines not in the right location :

Answer 1

You would have probably have had answers a lot earlier if you had made a simple reprex, e.g. like this

library(raster)
r <- raster(res=1/12)
values(r) <- sample(100, ncell(r), replace=TRUE)
filename <- system.file("external/lux.shp", package="raster")
v <- shapefile(filename)

1. zoom in a certain area

One way to zoom is to use crop (alternatively use the ext argument in plot)

x <- crop(r, v)

2. Display raster r values on a scale with cuts intervals

cuts <- c(0,20,60,100)
plot(x, breaks=cuts, col=rainbow(3))

or

 y <- cut(x, cuts) 

3. Plot the country boundaries

lines(v)

4. Finally, include polygon layer (with col="blue" fill and contours).

plot(v[c(1,3),], col="blue", border="red", lwd=2, add=TRUE)

Answer 2

• 6 months later but I feel this question. My two thoughts are (1) I have had luck with plotting geom_sf and geom_stars together. You have to change your raster to a df before changing to a geom_stars. and (2) regardless of method, you need all datasets in the same projection - check with crs() and set all to the same with st_transform()

I didn't actually test this with your data but something like:

make raster into a df

test.df = as.data.frame (test, xy=TRUE) # Convert to data.frame, keeping the 
coordinates
class(test.df)

convert to geom_stars

test.stars = st_as_stars(test.df)

try your plot

    Plot1 <- ggplot()+

  geom_stars(data = test, aes(fill=values))+ #need to plot raster first I think?
scale_fill_identity( name = "", breaks = cuts,labels = "")+      

  geom_sf(data = shpAfr.fort, color ="grey27", size = .2)+
geom_sf(data = shpWater.fort, color ="lightskyblue2", fill 
  ="lightskyblue2", size = .2)+

  theme_minimal()+
  coord_sf( xlim = NIG[1:2], ylim = NIG[3:4]),expand = FALSE)

Plot1