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

Question

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

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

这接近我的目标，但在 mapview 中绘图而不是 ggplot2：

这是我的尝试：

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 和地块

的以太平洋为中心的世界地图

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 创建

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

Plotting convex hulls crossing 180 degree international date line and calculating area

mapping

maps

r

ggplot2

sf