在正投影中繪製世界地圖時剪切多邊形

我嘗試使用從here修改的代碼在R中繪製世界地圖，正如下面給出的。還顯示了產出數字 - 很明顯，在邊界附近的陸地地區被削減，在這種情況下，俄羅斯和南極洲。我相信這是由於多邊形上的某些點包裝到可見側的「後面」，它們通過映射函數轉換爲NAs。有什麼辦法可以解決這個問題嗎？在正投影中繪製世界地圖時剪切多邊形

我真的需要那些缺失的區域，因爲我的最終目標是繪製幾個這樣的地圖，每個地圖的中心點稍有不同。如果某些大地塊隨意進出，那看起來會很奇怪。

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)

來源

2016-03-10 David

附註：我不知道爲什麼，但我根本無法把在「映射」標籤 - 它someh ow自動更改爲「詞典」... – David

是的，它的工作。謝謝。 – David

我通過投影是剛出來的觀點地球外的點，並與這些「出」點填充可見多邊形一起解決它。然後用黑色多邊形覆蓋溢出的陸地。

#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) 
}

我也寫投影功能：

爲了做到這一點，我使用下面的函數（從地圖包修改）中檢索連續的陸塊的座標。「前」 = 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和岑隆是在可見地球的中心的座標）的代碼：

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

結果看起來是這樣的（應用在全球範圍內的黑色層前）：

來源

2016-03-12 23:19:43 David

你是怎麼繪製的？導致答案沒有繪圖功能 – user3617715

在正投影中繪製世界地圖時剪切多邊形

回答

相關問題