geom_polygon + 地图投影 = 莫名其妙切成两种形状？
Posted
技术标签:
【中文标题】geom_polygon + 地图投影 = 莫名其妙切成两种形状？【英文标题】：geom_polygon + map projection = inexplicably cut into two shapes? 【发布时间】：2016-08-16 01:59:22 【问题描述】：我正在尝试使用 ggplot2 在 Winkel Tripel 投影中绘制世界地图；它最终将有一些数据。就我所知，ggplot 本身不能做 Winkel Tripel，所以我用手动预测来解决这个问题。我什么都有，除了海洋层，它不正确。代码：
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())
渲染：
您可以看到原本应该是单个多边形填充的内容已被分割成两个单独的多边形，它们仅覆盖地图的“端盖”，而不是中间。如何让它填满整个地图？我认为问题出在“边界”的定义上，但我在 geom_polygon 文档中没有看到任何内容来解释可能有什么问题。
【问题讨论】：
我认为这可能是一个蝴蝶结多边形，即 ll, ul, lr, ur, ll 而不是围绕 ll, ul, ur, lr, ll 的正确路径，它们代表上/下，左右。现在无法测试自己 @mdsumner 这就是第二次反转 lseq 的原因。 好吧，抱歉，顺便说一句，这是 scale_y_continuous 步骤，可能会以有问题的方式破坏多边形 - 您可以通过注释该行来看到它很好（但对于 y 范围） @mdsumner 我在发布之前尝试过，但并没有改变任何事情。奇怪的。这篇论文完成后，我可能会再玩一些 :) 【参考方案1】：从ggalt 开始，您可以直接进行那些基于 PROJ.4 的预测：
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
你原来的真正问题是你的多边形限制，因为你用你的边界框在南极洲的顶部划伤：
lseq = seq(-53, 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
尽管coord_proj() 无论如何都可以：
gg + coord_proj("+proj=wintri", ylim=c(-60, 85))
所以，稍微调整一下：
lseq = seq(-53, 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
而且不用直接与企鹅交往也能正常工作：
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", ylim=c(-53, 85))
gg <- gg + theme_map()
gg
【讨论】：
我怀疑子集 + 不使用 scale_y_limit 是真正的解决方法 更新了更多信息/示例（超级严重的感冒正在发生，所以我昨晚太简短了……抱歉）。【参考方案2】：虽然使用 'ggalt' 是“正确答案”，但这里有一种方法可以得到你想要的 - 首先裁剪出南极洲（虽然我不支持，但它真的很糟糕），然后使用 'scale_y_discrete'。我对 gg-verse 还不够熟悉，不知道为什么，但我收集到连续的作物切断了南部边界，并以某种方式再次绘制了围绕两半的路径。离散裁剪可能会选择“最近”的坐标，因此保持南部边界。
我们确实需要更好的工具来跨越 sp/ggplot2 的鸿沟，虽然我们对局部问题进行了部分修复，但两者并没有适当的统一。 'ggalt' 是最好的之一。
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")))
## not much of a "world" without Antarctica . .
badworld <- subset(spTransform(getMap(), CRS("+proj=wintri")),
!NAME == "Antarctica")
world <- fortify(badworld)
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_discrete(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())
【讨论】：
以上是关于geom_polygon + 地图投影 = 莫名其妙切成两种形状？的主要内容，如果未能解决你的问题，请参考以下文章