5

I have a raster file 'airtemp' and a polygon shapefile 'continents'. I'd like to superimpose the 'continents' on 'airtemp', so the boundary of 'continents' is visible on top of 'airtemp'. I plot the raster file by levelplot (lattice). I read the polygon by readShapeSpatial (maptools) first and then plot.

The problem is levelplot and plot have different scales. Plot tends to have smaller frame. Sorry I don't have a reproducible sample, but I feel this is a rather common issue for geophysicists. I've found a similar question here:

http://r.789695.n4.nabble.com/overlaying-a-levelplot-on-a-map-plot-td2019419.html

but I don't quite understand the solution.

EDU
  • 267
  • 1
  • 5
  • 13
  • 1
    The answer say , that `levelplot` is a lattice function, `plot` is a base one, very hard to mix base and grid graphics. – agstudy Jul 11 '13 at 02:29

2 Answers2

12

You can overlay the shapefile using the +.trellis and layer functions from the latticeExtra package (which is automatically loaded with rasterVis).

library(raster)
library(rasterVis)

Let's build some data to play. You can skip this part if you already have a raster file and a shapefile.

library(maps)
library(mapdata)
library(maptools)

## raster
myRaster <- raster(xmn=-100, xmx=100, ymn=-60, ymx=60)
myRaster <- init(myRaster, runif)

## polygon shapefile
ext <- as.vector(extent(myRaster))

boundaries <- map('worldHires', fill=TRUE,
    xlim=ext[1:2], ylim=ext[3:4],
    plot=FALSE)

## read the map2SpatialPolygons help page for details
IDs <- sapply(strsplit(boundaries$names, ":"), function(x) x[1])
bPols <- map2SpatialPolygons(boundaries, IDs=IDs,
                              proj4string=CRS(projection(myRaster)))

Now you plot the raster file with rasterVis::levelplot, the shapefile with sp::sp.polygons, and the overall graphic is produced with +.trellis and layer.

levelplot(myRaster) + layer(sp.polygons(bPols))

overlay with transparent color

sp.polygons uses a transparent color as default for fill, but you can change it:

levelplot(myRaster) + layer(sp.polygons(bPols, fill='white', alpha=0.3))

overlay with white color

Oscar Perpiñán
  • 4,491
  • 17
  • 28
  • 1
    Using package tidyverse with this example might create conflict for `map` and `layer`. If you get Error: Attempted to create layer with no geom, use `latticeExtra::layer`. – Matifou Sep 29 '17 at 21:14
  • what are the "histograms" of the X and Y axis? how can I remove them? – MilloMarinE Nov 27 '18 at 01:01
  • It is documented in the help page of rasterVis::levelplot, the "margin" argument. – Oscar Perpiñán Nov 27 '18 at 09:26
  • @OscarPerpiñán would it be possible to fill the polygons with lines (like if you have to shade a polygon)? is there a specific `fill` option? – Nemesi Mar 28 '19 at 15:08
  • 1
    @Nemesi That option is not implemented, but you may define your own panel function to get it (for example: https://stackoverflow.com/a/9422480/964866) – Oscar Perpiñán Mar 29 '19 at 10:52
1

As per this discussion, here is one way of doing this: it consists in breaking the SpatialPolygonsDataFrame into one single matrix of polygons coordinates separated by NAs. Then plotting this on the levelplot using panel.polygon.

library(maptools)
a <- matrix(rnorm(360*180),nrow=360,ncol=180) #Some random data (=your airtemp)
b <- readShapeSpatial("110-m_land.shp") #I used here a world map from Natural Earth.

And that's where the fun begins:

lb <- as(b, "SpatialPolygons")
llb <- slot(lb, "polygons")
B <- lapply(llb, slot, "Polygons") #At this point we have a list of SpatialPolygons
coords <- matrix(nrow=0, ncol=2)
for (i in seq_along(B)){
    for (j in seq_along(B[[i]])) {
        crds <- rbind(slot(B[[i]][[j]], "coords"), c(NA, NA)) #the NAs are used to separate the lines
        coords <- rbind(coords, crds)
        }
    }
coords[,1] <- coords[,1]+180 # Because here your levelplot will be ranging from 0 to 360°
coords[,2] <- coords[,2]+90 # and 0 to 180° instead of -180 to 180 and -90 to 90

And then comes the plotting:

levelplot(a, panel=function(...){
                        panel.levelplot(...)
                        panel.polygon(coords)})

The idea in lattice is to define the plotting functions in argument panel(see ?xyplot for a complete explanation on the subject). The function for the levelplot itself is levelplot.

enter image description here

Of course, in your case, it seems way simpler to plot this using base graphics:

image(seq(-180,180,by=1),seq(-90,90,by=1),a)
plot(b, add=TRUE)

enter image description here

plannapus
  • 18,529
  • 4
  • 72
  • 94