在正射投影下绘制世界地图出现了“非有限点”问题。

Question

我有一个世界各国的shapefile，可以从这里下载。我可以使用R语言绘制它.

countries <- readOGR("shp","TM_WORLD_BORDERS-0.3",encoding="UTF-8",stringsAsFactors=F)
par(mar=c(0,0,0,0),bg=rgb(0.3,0.4,1))
plot(countries,col=rgb(1,0.8,0.4))

现在我想在正射投影中绘制它（从外太空看到的地球），所以我正在尝试。

countries <- spTransform(countries,CRS("+proj=ortho +lat_0=-10 +lon_0=-60"))

我也尝试调整x_0和y_0参数（如此处所述），但总是出现错误：

non finite transformation detected:
[1] 45.08332 39.76804      Inf      Inf
Erro em .spTransform_Polygon(input[[i]], to_args = to_args, from_args = from_args,  : 
  failure in Polygons 3 Polygon 1 points 1
Além disso: Mensagens de aviso perdidas:
In .spTransform_Polygon(input[[i]], to_args = to_args, from_args = from_args,  :
  108 projected point(s) not finite

有时在第3个多边形中，有时在第7个中。这些“Inf”是从哪里来的？我需要更改任何参数吗？我想绘制出像这样的地图，但是要居中于南美洲之上。感谢您的帮助！

- Rodrigo

---

你正在尝试将所有坐标从-180到180，从-90到90进行投影。但是使用你尝试使用的投影方式不可能实现。首先尝试裁剪要看到投影的区域。我强烈建议你阅读一些有关地理投影的文献。 - user3710546

请参阅：http://en.wikipedia.org/wiki/Orthographic_projection_in_cartography#Mathematics，了解在何处裁剪的经验法则。 - plannapus

谢谢，Pascal。我想绘制整个半球世界，就像上面的图片一样。我认为我只需要中心纬度/经度，算法会自动裁剪所需部分。 - Rodrigo

在第一个例子中（http://matplotlib.org/basemap/users/examples.html），他们不需要做你所说的。这证明了有一种更简单的方法，不是吗？ - Rodrigo

另一个要考虑的问题是，裁剪被地球遮挡的点会破坏多边形的拓扑结构。 - Rodrigo

Answer 1

尝试使用maps包。它会警告那些无法投影的点，但不会出错并关闭进程。稍作调整，具体来说是设置海洋的填充颜色 (这个答案有助于解决这个问题)，用几行代码就能够模拟你附加的地图：

library(maps)

## start plot & extract coordinates from orthographic map
o <- c(-10,-60,0) # oreantation
xy <- map("world",proj="orthographic", orientation=o, bg="black")
xy <- na.omit(data.frame(do.call(cbind, xy[c("x","y")])))

## draw a circle around the points for coloring the ocean 
polygon(max(xy$x)*sin(seq(0,2*pi,length.out=100)),max(xy$y)*cos(seq(0,2*pi,length.out=100)), 
        col="blue4", border=rgb(1,1,1,0.5), lwd=2)

## overlay world map
colRamp <- colorRampPalette(c("lemonchiffon", "orangered"))
map("world",proj="orthographic", orientation=o, 
    fill=TRUE, col=colRamp(5), add=TRUE)

- Paul Regular

---

谢谢Paul，它正在工作。我正在为一个插图做这个，给出另一个更大的地图的位置。如何将另一个地图的边界框投影为此球体上的矩形？（它不会是一个矩形，而是一个弯曲的矩形...）我可以使用mapproject和rect来完成，但然后我得到一个非弯曲的矩形。 - Rodrigo

可能有更简单的方法，但我会手动构建矩形的坐标，使用mapproject进行投影，并使用polygon添加它：x <- seq(-60, -30, length.out=100); y <- seq(-30, -10, length.out=100); x <- c(x, rep(x[length(x)], length(x)), rev(x), rep(x[1], length(x))); y <- c(rep(y[1], length(y)), y, rep(y[length(y)], length(y)), rev(y)); xy <- mapproject(x, y, proj="orthographic", orientation=o); polygon(xy$x, xy$y) - Paul Regular

是的，当然，我怎么没想到呢？;) 再次感谢！ - Rodrigo