0

I am trying to do a levelplot with a Pacific centered map.

This is the code I am using

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")}
)

and this is what I get. The polygons seem to be all messed up.

enter image description here

user3910073
  • 511
  • 1
  • 6
  • 23
  • what would you want to see? Show your data, etc – Amit Kohli Oct 02 '15 at 12:00
  • @AmitKohli the data are irrelevant I would like to be able to see the map entered in the pacific – user3910073 Oct 02 '15 at 12:43
  • Possible duplicate of [Plot Map of Pacific with filled countries](http://stackoverflow.com/questions/18994274/plot-map-of-pacific-with-filled-countries) – Cotton.Rockwood Feb 02 '16 at 20:23
  • Hi @user3910073... this is a problem with trying to draw polygons across the dateline. The question has already been addressed in several places: http://stackoverflow.com/q/18994274/3897439 and http://stackoverflow.com/q/10620862/3897439 – Cotton.Rockwood Feb 02 '16 at 20:25

1 Answers1

0

Took me way too long to find the answer to this question over here by Michael Sumner. I've only added an spTransform to ensure the 2 pieces go back together with 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'))

Pac centric map that doesn't have lines across it