在正射投影中绘制世界地图时裁剪多边形

Question

我尝试在 R 中使用正交投影绘制世界地图，使用从 here 修改的代码，如下所示。还显示了输出图 - 很明显边界附近的陆地区域被裁剪了，在这种情况下是俄罗斯和南极洲。我相信这是由于多边形上的一些点环绕到可见侧的 "back"，这些点被映射函数转换为 NA。有什么办法可以解决这个问题吗？

我真的需要那些缺失的区域，因为我的最终目标是绘制其中几张地图，每张地图的中心点都略有不同。如果某些陆地随意进出，那看起来会很奇怪。

library(maps)
library(mapdata)

## start plot & extract coordinates from orthographic map
o <- c(-5,155,0) #orientation
xy <- map("world",proj="orthographic",orientation=o)
## draw a circle around the points for coloring the ocean 
polygon(sin(seq(0,2*pi,length.out=100)),cos(seq(0,2*pi,length.out=100)),col=rgb(0.6,0.6,0.9),border=rgb(1,1,1,0.5),lwd=2)
## overlay world map
map("worldHires",proj="orthographic",orientation=o,fill=TRUE,col=rgb(0.5,0.8,0.5),resolution=0,add=TRUE)

Answer 1

我通过将刚好在视野之外的点投影到地球之外，并用这些 "out" 点填充可见的多边形来解决这个问题。然后使用黑色多边形覆盖溢出的陆地。

为了做到这一点，我使用以下函数（从地图包修改）检索了连续大陆的坐标：

#Get contiguous country coordinates
contigcoord <- function (database = "world", regions = ".", exact = FALSE, boundary = TRUE, interior = TRUE, fill = FALSE, xlim = NULL, ylim = NULL){
    if (is.character(database))
        as.polygon = fill
    else as.polygon = TRUE
    coord <- maps:::map.poly(database, regions, exact, xlim, ylim, boundary, 
        interior, fill, as.polygon)
    return(coord)
}

我也写了投影函数。 "front" = 1 如果该点应该可见：

#Lat/Long to x-y (0,1 range) in orthographic projection. Cen is the centre point of map
mapproj <- function(lat,long,cenlat,cenlong){
    d2r=pi/180; lat=lat*d2r; long=long*d2r; cenlat=cenlat*d2r; cenlong=cenlong*d2r
    x=cos(lat)*sin(long-cenlong)
    y=cos(cenlat)*sin(lat)-sin(cenlat)*cos(lat)*cos(long-cenlong)
    front=sin(cenlat)*sin(lat)+cos(cenlat)*cos(lat)*cos(long-cenlong) > 0
    return(list(x=x,y=y,front=front))
}

用于填充陆地的代码如下（cenlat和cenlong是可见地球中心的坐标）：

xy <- contigcoord("world",fill=TRUE)
coord <- cbind(xy$x,xy$y)
coordtr <- mapproj(coord[,2],coord[,1],cenlat,cenlong)
coord <- cbind(coord,coordtr$x,coordtr$y,coordtr$front)

naloc <- (1:nrow(coord))[!complete.cases(coord)]
naloc <- c(0,naloc)
for(i in 2:length(naloc)){
    thispoly <- coord[(naloc[i-1]+1):(naloc[i]-1),3:5,drop=F]
    thispoly <- rbind(thispoly,thispoly[1,])
    unq <- unique(thispoly[,3])
    if(length(unq) == 1){ 
        if(unq == 1){ #Polygon is fully on front side
            polygon(thispoly[,1],thispoly[,2],col=rgb(0.5,0.8,0.5),border=NA)
        }
    } else { #front and back present
        ind <- thispoly[,3] == 0
        #project points "outside" the globe
        temdist <- pmax(sqrt(rowSums(as.matrix(thispoly[ind,1:2]^2))),1e-5)
        thispoly[ind,1:2] <- thispoly[ind,1:2]*(2-temdist)/temdist
        polygon(thispoly[,1],thispoly[,2],col=rgb(0.5,0.8,0.5),border=NA)
    }
}

结果如下所示（在全球范围内应用黑色图层之前）：

在正射投影中绘制世界地图时裁剪多边形

Clipping polygons when drawing world map in orthographic projection

maps

r

orthographic