防止海岸線的局部映射

Question

如何防止國家多邊形在不同預測下被切斷？

在下面的例子中，我想做一個南極洲的立體投影地圖，包括緯度< -45°S。通過將我的y限制設置爲這個範圍，繪圖區域是正確的，但是國家多邊形也在這些極限處裁剪。有沒有辦法將海岸線繪製到繪圖區域的邊緣？

感謝您的任何建議。

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() 

Answer 1

的問題是，因爲你是一個specifing最大的-45 ylim，數據正在以緯度精確地切割。很明顯，情節的實際界限是以一種單獨的方式計算的，這可能是基於由此產生的預計限制（但是我沒有探索源代碼我一無所知）進行某種保守的假設。

您可以避免通過建立一個虛擬的情節不繪製海岸線，然後添加數據的頂部增加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() 

順便說一句，這個策略不會爲工作所有的預測，因爲一些意味着重疊的一些數據限制，但極地立體攝影只是從中心延伸出來，所以這裏沒有這樣的問題。

此外，將會有一個更新的方式來做到這一點與maptools，rgdal和rgeos，這將允許您使用數據共享一個合適的PROJ.4立體投影（但我還沒有弄清楚削減步）。 mapproj中的預測是「僅用於形狀」，將其他數據轉化爲afaik將是一項工作。

Answer 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() 

Answer 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()