zwol zwol - 3 months ago 17
R Question

geom_polygon + map projection = inexplicably cut into two shapes?

I'm trying to draw a world map in Winkel Tripel projection, using ggplot2; it will ultimately have some data on top of it. Natively, as far as I know, ggplot can't do Winkel Tripel, so I have kludged around this with manual projections. I've got everything except the ocean layer, which doesn't come out right. Code:

suppressPackageStartupMessages({
library(ggplot2)
library(sp)
library(rworldmap)
library(rgdal)
})
ll.to.wt <- function (points)
as.data.frame(spTransform(SpatialPoints(points, CRS("+proj=longlat")),
CRS("+proj=wintri")))

world <- fortify(spTransform(getMap(), CRS("+proj=wintri")))
xlimits <- ll.to.wt(matrix(c(-180,180,0,0), nrow=2))$coords.x1
ylimits <- ll.to.wt(matrix(c(0,0,-60,85), nrow=2))$coords.x2
lseq = seq(-60, 85, by=.25)
boundary <- ll.to.wt(data.frame(
long = c(rep(-180, length(lseq)), rep(180, length(lseq)), -180),
lat = c(lseq, rev(lseq), lseq[1])))

ggplot() +
geom_polygon(data=boundary, aes(x=long, y=lat), fill="#9AC5E3") +
geom_map(data=world, map=world, aes(x=long, y=lat, map_id=id),
color="#888888", fill="#f2caae", size=0.25) +
scale_x_continuous(limits=xlimits, expand=c(0,0)) +
scale_y_continuous(limits=ylimits, expand=c(0,0)) +
coord_equal() +
theme(
axis.line=element_blank(),
axis.text.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
legend.justification = c(0,0), # bottom of box
legend.position = c(0,0), # bottom of picture
panel.background=element_blank(),
panel.border=element_blank(),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.margin=unit(0, "lines"),
plot.background=element_blank())


Rendered:

world map drawn correctly except that the blue ocean fill only covers the "end caps" of the map, not the region in the middle

You can see that what was meant to be a single polygon fill has been chopped into two separate polygons, and they only cover the "end caps" of the map, not the middle. How do I make it fill under the entire map? I presume that the problem is with the definition of "boundary", but I don't see anything in the
geom_polygon
documentation to explain what might be wrong.

Answer

Ever since ggalt you can do those PROJ.4-based projections directly:

library(ggplot2)
library(sp)
library(rworldmap)
library(rgdal)
library(ggalt)
library(ggthemes)

world <- fortify(getMap())
world <- subset(world, id != "Antarctica")
lseq = seq(-60, 85, by=.25)
boundary <- data.frame(
    long = c(rep(-180, length(lseq)), rep(180, length(lseq)), -180),
    lat  = c(lseq,                    rev(lseq),          lseq[1]))

gg <- ggplot()
gg <- gg + geom_polygon(data=boundary, aes(x=long, y=lat), fill="#9AC5E3")
gg <- gg + geom_map(data=world, map=world, 
                    aes(long, lat, map_id=id),
                    color="#888888", fill="#f2caae", size=0.25)
gg <- gg + coord_proj("+proj=wintri")
gg <- gg + theme_map()
gg

enter image description here