绘制跨越 180 度国际日期变更线的凸包并计算面积

Question

我正在尝试使用凸包绘制物种范围区域，然后计算面积并创建一个图形。

180 度国际日期变更线有一个众所周知的问题，我一直在尝试解决 SE 上的许多示例，例如：

How to remedy a path that crosses the international dateline with R

这接近我的目标，但在 mapview 中绘制的不是 ggplot2： How to construct/plot convex hulls of polygons from points by factor using sf?

这是我的尝试：

library(tidyverse)
library(maps)
library(ggmap)
library(sf)
library(sp)
library(rnaturalearth)
library(rnaturalearthdata)
library(ggspatial)
library(mapproj)

生成物种出现数据，其中一些点跨越 180 经度

df <- data.frame(species = rep("sp1",8),
                 longitude = as.double(c(-170.2, -179.5, 55.9, 167.6, 154.3, 101.7, 70.54, -165.94)),
                  latitude = as.double(c(8.25, -24.75, 24.25,19.25, 33.45, -15.5, 5.56, 4.6)))

来自 map_data 和 plot 的以太平洋为中心的世界地图

world <- map_data("world2") 

map<-ggplot() +
geom_polygon(data = world, aes(x = long, y = lat, group = group),
col = "#78909C", fill = "#78909C", lwd = 0)+
coord_map(orientation = c(90,0, 150), ylim = c(-40, 40), xlim = c(20,210))

my map

向地图添加发生点

map +
geom_point(data = df, mapping = aes(x = longitude, y = latitude))

map with points

根据物种出现数据构造最小凸包。

species.sf <- df %>%
  st_as_sf( coords = c( "longitude", "latitude" ))

创建船体并环绕日期线

hull<- species.sf %>%
  summarise( geometry = st_combine( geometry ) ) %>%
  st_convex_hull()

hull<-st_wrap_dateline(hull,options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180"),
   quiet = TRUE)

绘制船体 - 在 180 处切割，但显然不包括所有出现点

map +
geom_point(data = df, mapping = aes(x = longitude, y = latitude))+
geom_sf(data=hull, inherit.aes = TRUE)

incorrect hull

计算船体面积 - 必须根据船体形状不正确

st_area(hull)

我还尝试将太平洋中心 CRS 应用于地图、点和船体，但怀疑我以错误的顺序或错误的位置应用这些？我对使用 R 进行空间分析非常陌生，所以任何帮助都非常有用。谢谢。

Answer 1

请在下面找到您的问题的解决方案。我使用了 sf 包中的函数 st_shift_longitude()。

Reprex

• 您的数据（无更改）

df <- data.frame(species = rep("sp1",8),
                 longitude = as.double(c(-170.2, -179.5, 55.9, 167.6, 154.3, 101.7, 70.54, -165.94)),
                 latitude = as.double(c(8.25, -24.75, 24.25,19.25, 33.45, -15.5, 5.56, 4.6)))



world <- map_data("world2") 

map<-ggplot() +
  geom_polygon(data = world, aes(x = long, y = lat, group = group),
               col = "#78909C", fill = "#78909C", lwd = 0) + 
  coord_map(orientation = c(90,0, 150), ylim = c(-40, 40), xlim = c(20,210))
map

• 将dataframe“df”转换为sf对象“species.sf”并用st_shift_longitude()移动经度

species.sf <- df %>% 
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>% 
  st_shift_longitude()


map + 
  geom_sf(data = species.sf, inherit.aes = TRUE) + 
  coord_sf(xlim = c(40, 210), ylim = c(-40, 40))

• 基于sf 对象“species.sf”和group_by(species) 计算凸包多边形（适用于您的一般情况）

hull <- species.sf %>%
  group_by(species) %>% 
  summarise( geometry = st_combine( geometry ) ) %>%
  st_convex_hull()

• 将sf对象“hull”转换回dataframe对象“hullDF”

hullDF <- hull %>% 
  st_geometry() %>% 
  st_coordinates() %>% 
  as.data.frame() 

• 最终结果的可视化

map + 
  geom_point(data = df, mapping = aes(x = longitude, y = latitude)) +
  geom_polygon(data = hullDF, mapping = aes(x = X, y = Y), fill = "lightgreen", alpha = 0.5)

• 计算船体多边形的面积（需要units 库将结果转换为平方千米）

library(units)

hull_area <- hull %>%  
  st_area() %>% 
  set_units(km^2)

hull_area
#> 74714882 [km^2]

^{由reprex package (v2.0.1) 于 2021 年 11 月 12 日创建}