以太平洋为中心的Levelplot图

Question

我正在尝试做一个以太平洋为中心的levelplot地图。

这是我正在使用的代码

levelplot(x, col.regions=colorscale,              
          panel = function(x, y, ...) {
            panel.levelplot(x, y,  ...)
            mp <- map("world2", plot = FALSE, fill=TRUE,interior = FALSE,bg="white")
            lpolygon(mp$x, mp$y, fill="black", col="black")}
)

这就是我所得到的。多边形看起来都乱七八糟的。

​

​

Answer 1

我花了很长时间才通过Michael Sumner的here找到这个问题的答案。我只是添加了一个spTransform，以确保这两个片段可以与spRbind放在一起

library(maptools)
library(raster)


## first, do the jiggery pokery of wrld_simpl to
## convert from Atlantic to Pacific view
xrange <- c(0, 360)
yrange <- c(-90, 90)

require(maptools)
require(raster)
require(rgeos)

ext <- extent(xrange[1], xrange[2], yrange[1], yrange[2])
data(wrld_simpl, package = "maptools")

## this is a bit more general than needed for this example, since I
## wanted arbitrary crops originally

eastrange <- c(xrange[1], 180.0)
westrange <- c(-180.0, xrange[2] - 360)

ew <- extent(westrange[1], westrange[2], yrange[1], yrange[2])
ee <- extent(eastrange[1], eastrange[2], yrange[1], yrange[2])

geome <- as(ee, "SpatialPolygons")
geomw <- as(ew, "SpatialPolygons")

## why does this get dropped?
proj4string(geome) <- CRS(proj4string(wrld_simpl))
proj4string(geomw) <- CRS(proj4string(wrld_simpl))

worlde <- gIntersection(wrld_simpl, geome)
worldw <- gIntersection(wrld_simpl, geomw)
worldw <- elide(worldw, shift = c(360.0, 0.0))

proj4string(worldw) <- CRS(proj4string(wrld_simpl))

dfe <- data.frame(x = seq_len(length(worlde@polygons)))
dfw <- data.frame(x = seq_len(length(worldw@polygons)))
rownames(dfw) <- as.character(as.numeric(rownames(dfe)) + nrow(dfe))

worlde <- SpatialPolygonsDataFrame(worlde, dfe, match.ID = FALSE)
worldw <- SpatialPolygonsDataFrame(worldw, dfw, match.ID = FALSE)
worldw <- spChFIDs(worldw, rownames(dfw))

## I had to add this spTransform() call to stop the spRbind proj error
worldw_proj <- spTransform(worldw, CRSobj = worlde@proj4string)
world <- spRbind(worlde, worldw_proj)

## so, "world" is longitudes [0, 360], which is pretty much essential
## given the input data shown

r <- raster(nrows = 180, ncols = 360, xmn = 0, xmx = 360, ymn = -90,
            ymx = 90, crs = "+proj=longlat +ellps=WGS84 +over")

## fake a raster so we have something equivalent to your data (but simple/r)
rast <- rasterize(world, r)

## fill the ocean parts with something . . .
rast[is.na(rast)] <- 1:sum(is.na(getValues(rast)))

library(rasterVis)
levelplot(rast, margin = FALSE) + layer(sp.polygons(world, lwd=0.8,
                                                    col='darkgray'))

​

​