R：水平图上的叠加图答案

【问题标题】：R: overlay plot on levelplotR：水平图上的叠加图
【发布时间】：2013-07-10 23:15:11
【问题描述】：

我有一个光栅文件 'airtemp' 和一个多边形 shapefile 'continents'。我想将“大陆”叠加在“airtemp”上，因此“大陆”的边界在“airtemp”上可见。我通过levelplot（格子）绘制光栅文件。我首先通过readShapeSpatial (maptools) 读取多边形，然后是plot。

问题是levelplot 和plot 有不同的比例。 Plot 往往有更小的框架。抱歉，我没有可重复的样本，但我觉得这对地球物理学家来说是一个相当普遍的问题。我在这里发现了一个类似的问题：

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

但我不太明白解决方案。

【问题讨论】：

答案是，levelplot 是一个格函数，plot 是一个基函数，很难混合基和网格图形。

标签： r plot levelplot

【解决方案1】：

您可以使用 +.trellis 和 layer 覆盖 shapefile latticeExtra 包中的函数（自动加载了rasterVis)。

library(raster)
library(rasterVis)

让我们构建一些数据来玩。如果你已经可以跳过这部分有一个光栅文件和一个 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)))

现在您使用rasterVis::levelplot 绘制光栅文件，即带有sp::sp.polygons的shapefile，就生成了整体图形 +.trellis 和 layer。

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

sp.polygons 使用透明颜色作为fill 的默认颜色，但您可以更改它：

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

【讨论】：

在此示例中使用包 tidyverse 可能会导致 map 和 layer 发生冲突。如果出现错误：尝试创建没有几何图形的图层，请使用 latticeExtra::layer。
什么是 X 轴和 Y 轴的“直方图”？我怎样才能删除它们？
在 rasterVis::levelplot 的帮助页面中记录了“margin”参数。
@OscarPerpiñán 是否可以用线条填充多边形（例如必须为多边形着色）？有特定的fill 选项吗？
@Nemesi 那个选项没有实现，但是你可以定义自己的面板函数来获取它（例如：stackoverflow.com/a/9422480/964866）

【解决方案2】：

根据this discussion，这是执行此操作的一种方法：它将 SpatialPolygonsDataFrame 分解为一个由 NA 分隔的多边形坐标矩阵。然后使用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.

这就是乐趣的开始：

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

然后是绘图：

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

lattice 的想法是在参数panel 中定义绘图函数（有关该主题的完整说明，请参见?xyplot）。 levelplot 本身的函数是levelplot。

当然，在你的情况下，使用base 图形绘制它似乎更简单：

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

【讨论】：