在 R 中创建几种植物物种的 Robinson 分布图 - 多边形错误和其他问题

Question

这是我第一次在 ggplot 中使用地图。我需要创建 90 个植物物种分布图。

因此，我设置了 API 以从 Plant of the World Online 获取数据，并从 World Geographical Scheme for Recording Plant Distributions https://github.com/tdwg/wgsrpd 获取 shp

require(kewr)#to connect to POWO APIs
 


#getting data from POWO

id.powo <- search_wcvp("Camellia japonica")
id.powo<- id.powo$results[[1]]$id

r <- lookup_powo(id.powo, distribution = TRUE)

native <- r$distribution$natives
introduced <- r$distribution$introduced

tdwg.native.name <- list()
for (i in 1:length(native)){
  tdwg.native.name[[i]]<- native[[i]]$name
}
tdwg.native.name <- unlist(tdwg.native.name)

col.native <- c("#B8DE95") #pastel green for native 


tdwg.introduced.name <- list()
for (i in 1:length(introduced)){
  tdwg.introduced.code[[i]]<- introduced[[i]]$name
}
tdwg.introduced.name <- unlist(tdwg.introduced.code)

#col.introduced <- c("#9F6CCC") #pastel violet  for introduced


#creating the map

install.packages("remotes")
remotes::install_github("barnabywalker/bazr")
library(bazr)

I downloaded the shp from https://github.com/tdwg/wgsrpd and then loaded in r

tdwg_level3.shp<- read_sf("/Users/...d/wgsrpd-master/level3/level3.shp")

tdwg_level3.shp.robinson<- st_transform(tdwg_level3.shp,  crs ="+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs") 


ggplot() + 
  geom_sf(data = tdwg_level3.shp.robinson)

问题从这里开始，多边形形成北极重叠并且没有匹配

此外，地图看起来很拥挤，因为某些区域被分区，并且有小岛使地图看起来很脏。最后，没有必要绘制南极洲......

因此，所需的输出应如下所示

此地图是使用此脚本构建的https://gist.github.com/valentinitnelav/065af8eba2d9455d9407e5d3890f6f86

我需要制作大约 100 张地图，每个物种一张。我将在数据框列中列出代码 id.powo 然后循环代码。

id.powo.list<- c("17414550-1", "296290-1", "263221-1")

谢谢

Answer 1

你很亲近！我会依次回答你的问题，但首先要说我已经在 terra 包中完成了这一切，因为我发现它是迄今为止最直观的空间内容。

require(terra)
require(tidyterra)
tdwg_level3.shp <- vect("/Users/...d/wgsrpd-master/level3/level3.shp")

1. 关于线条，这是因为有时当您重新投影世界的多边形时，新投影会在稍微不同的位置切割多边形（在这种情况下，俄罗斯/阿拉斯加之间的确切点，俄罗斯多边形的边缘被切割离开）。最简单的解决方法是在重新投影之前稍微修剪第一个多边形：

   tdwg_level3.shp <- crop(tdwg_level3.shp, extent(-179.9, 179.9, -90, 84))
   

2. 关于看起来很拥挤的地图，这很容易通过在 ggplot 中将“颜色”和“填充”设置为相同来解决，然后线条与一般颜色混合：

   ggplot() + 
     geom_spatvector(data = tdwg_level3.shp.robinson, fill="grey70", colour="grey70")
   

   World map with no boundary lines

   如果你想一些行（例如级别 1 的区别），那么首先您可以聚合 tdwg shapefile：

   tdwg_level3.shp.robinson,simple <- aggregate(tdwg_level3.shp.robinson, by="LEVEL2_COD")
   
   ggplot() + geom_spatvector(data = tdwg_level3.shp.robinson.simple)
   

   World map with aggregate boundary lines

3. 删除南极洲很容易，只需将 shapefile 子集化即可：

   tdwg_level3_rob <-tdwg_level3_rob[tdwg_level3_rob$LEVEL3_NAM != "Antarctica",]
   

4. 获取本地/引入的范围还有更多内容。首先，为每个创建新的 shapefile：

   tdwg_level3_nat <- tdwg_level3.shp.robinson[tdwg_level3.shp.robinson$LEVEL3_NAM %in% tdwg.native.name]
   
   tdwg_level3_int <- tdwg_level3.shp.robinson[tdwg_level3.shp.robinson$LEVEL3_NAM %in% tdwg.introduced.name]
   

   然后把它们放在一起（你会注意到必须做一些花哨的步法才能让图例起作用，但它们都放在一起）：

   ggplot() + 
     geom_spatvector(data = tdwg_level3_rob, fill="grey70", colour="grey70")+
     geom_spatvector(data = tdwg_level3_nat, aes(fill="Native", colour="Native"))+
     geom_spatvector(data = tdwg_level3_int, aes(fill="Introduced", colour="Introduced"))+
     scale_color_manual(values=colors)+
     scale_fill_manual(values=colors)+
     guides(color = "none")+ #This kills the "colour" legend, otherwise you get one for fill and one for colour
     theme(plot.background = element_blank(), panel.background = element_blank(), panel.grid=element_line(colour="grey90"),
       legend.position = c(0.2, 0.25), legend.background = element_rect(fill = NA, color = NA), 
       legend.title=element_blank())
   

   完整的分析（从下载形状时开始）最终是：

   tdwg_level3.shp <- vect(paste0(DataFP, "RandomCrap/tdwg wgsrpd master level3/level3.shp"))
   tdwg_level3.shp <- crop(tdwg_level3.shp, extent(-179.9, 179.9, -90, 84))
   tdwg_level3.shp.robinson <- project(tdwg_level3.shp,  "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
   tdwg_level3.shp.robinson <-tdwg_level3.shp.robinson[tdwg_level3.shp.robinson$LEVEL3_NAM != "Antarctica",]
   
   tdwg_level3_nat <- tdwg_level3.shp.robinson[tdwg_level3.shp.robinson$LEVEL3_NAM %in% tdwg.native.name]
   tdwg_level3_int <- tdwg_level3.shp.robinson[tdwg_level3.shp.robinson$LEVEL3_NAM %in% tdwg.introduced.name]
   
   colors <- c("Native" = col.native, "Introduced" = col.introduced) #Create a "colours" item so it can be used in the legend
   
   ggplot() + 
     geom_spatvector(data = tdwg_level3.shp.robinson, fill="grey70", colour="grey70")+
     geom_spatvector(data = tdwg_level3_nat, aes(fill="Native", colour="Native"))+
     geom_spatvector(data = tdwg_level3_int, aes(fill="Introduced", colour="Introduced"))+
     scale_color_manual(values=colors)+
     scale_fill_manual(values=colors)+
     guides(color = "none")+ #This kills the "colour" legend, otherwise you get one for fill and one for colour
     theme(plot.background = element_blank(), panel.background = element_blank(), panel.grid=element_line(colour="grey90"),
           legend.position = c(0.2, 0.25), legend.background = element_rect(fill = NA, color = NA), 
           legend.title=element_blank())
   

   （您可以使用 Valentin 链接来获取您想要的确切网格线等）

   Final Map

Answer 2

这是一个底图方法。这使用尚未在 CRAN 上的“terra 1.6.41”。你可以用install.packages('terra', repos='https://rspatial.r-universe.dev')安装这个版本

library(terra)
#terra 1.6.41
library(geodata)
w <- world(path=".")
w <- w[w$NAME_0 != "Antarctica", ]
w <- project(w, "+proj=robin")

w$value <- NA
w$value[w$NAME_0 == "China"]<- "Native"
w$value[w$NAME_0 == "India"]<- "Introduced"

g <- graticule(60, seq(-60,90,30), crs="+proj=robin")

plot(g, mar=0, background="azure", col="light gray", lab.cex=.4, off.lat=-.1, retro=T, lab.lon=2:7)
plot(w, "value", colNA=gray(.9), add=TRUE, col=c("light blue", "light green"), plg=list(x=-13800000, y=-3200000, cex=.7))