防止海岸线的部分测绘答案

【问题标题】：Preventing partial mapping of coastlines防止海岸线的部分测绘
【发布时间】：2012-07-15 08:22:39
【问题描述】：

如何防止国家多边形在不同的投影下被切断？

在下面的例子中，我想做一个南极洲的立体投影图，包括纬度

感谢您的建议。

library(maps)
library(mapproj)

ylim=c(-90,-45)
orientation=c(-90, 0, 0)

x11()
par(mar=c(1,1,1,1))
m <- map("world", plot=FALSE)
map("world",project="stereographic", orientation=orientation, ylim=ylim)
map.grid(m, nx=18,ny=18, col=8)
box()

【问题讨论】：

标签： r map-projections

【解决方案1】：

问题是因为您指定的最大 ylim 为 -45，并且数据在该纬度被准确切割。显然，绘图的实际边界是以单独的方式计算的，这可能与基于所产生的预计限制的某种保守假设有关（但如果不探索源代码，我对此一无所知）。

您可以通过设置一个不绘制海岸线的虚拟图来避免该问题，然后在顶部添加数据并增加ylim：

library(maps)
library(mapproj)

ylim = c(-90,-45)
orientation=c(-90, 0, 0)

x11()
par(mar=c(1,1,1,1))
m <- map("world", plot=FALSE)
map("world",project="stereographic", orientation=orientation, ylim=ylim, col = "transparent")
map("world",project="stereographic", orientation=orientation, ylim=ylim + c(-5, 5), add = TRUE)
map.grid(m, nx=18,ny=18, col=8)
box()

顺便说一句，此策略不适用于所有投影，因为有些投影意味着某些数据限制会重叠，但 Polar Stereographic 只是从中心向外延伸，因此这里不存在此类问题。

此外，使用 maptools、rgdal 和 rgeos 将有一种更新的方法来执行此操作，这将允许您使用共享正确 PROJ.4 立体投影的数据（但我还没有完全弄清楚剪裁步骤)。 mapproj 中的投影是“仅形状”的，将其他数据放入其中需要一些工作，afaik。

【讨论】：

【解决方案2】：

我意识到另一种选择是定义绘图区域的限制，而不是地图本身。这可能会在定义确切的绘图区域时提供一些灵活性（即 xaxs="i" 和 yaxs="i"）。您还保证将绘制所有多边形 - 在 @mdsumner 的示例中，缺少澳大利亚并且需要扩展 y 限制才能正确绘制它的多边形。

orientation=c(-90, 0, 0)

ylim <- c(mapproject(x=-180,y=-45, project="stereographic", orientation=orientation)$y, mapproject(x=0,y=-45, project="stereographic", orientation=orientation)$y)
xlim <- c(mapproject(x=-90,y=-45, project="stereographic", orientation=orientation)$x, mapproject(x=90,y=-45, project="stereographic", orientation=orientation)$x)

x11(width=6,height=6)
par(mar=c(1,1,1,1))
plot(0,0, t="n", ylim=ylim, xlim=xlim, xaxs="i", yaxs="i", xlab="", ylab="", xaxt="n", yaxt="n")
map("world",project="stereographic", orientation=orientation, add=TRUE)
map.grid(nx=18,ny=18, col=8)
box()

【讨论】：

【解决方案3】：

另一种方法是使用实际的 PROJ.4 投影并使用 maptools 包中的数据集。这是代码，在南极大陆仍然存在问题，海洋海岸线在 -90 处与极点相交（这有点烦人，并且需要找到该坐标并将其从转换结果中删除，或者通过拓扑工具运行来查找条子）。

library(maptools)
data(wrld_simpl)
xrange <- c(-180, 180, 0, 0)
yrange <- c(-90, -45, -90, -45)

stere <- "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

library(rgdal)

w <- spTransform(wrld_simpl, CRS(stere))

xy <- project(cbind(xrange, yrange), stere)
plot(w, xlim = range(xy[,1]), ylim = range(xy[,2]))
box()

【讨论】：