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.

4

2 回答 2

12

+.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::levelplotshapefile，使用和sp::sp.polygons生成整体图形。+.trellislayer

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

用透明颜色覆盖

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

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

用白色覆盖

于 2013-08-21T11:35:07.530 回答
1

根据这个讨论，这是一种方法：它包括将 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.

在此处输入图像描述

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

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

在此处输入图像描述

于 2013-07-12T07:03:41.230 回答