zwol - 2 months ago 9
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:

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.

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
``````