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

Question

I am trying to plot species range areas using convex hulls to then calculate the area and create a figure.

There is a well known issue with the 180 degree international dateline that I have been trying to remedy following many examples on SE, e.g:

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

This comes close to what I am aiming for but plots in mapview not ggplot2: How to construct/plot convex hulls of polygons from points by factor using sf?

Here is my attempt:

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

Generate species occurrence data, with some points crossing 180 longitude

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

Pacific centered world map from map_data and 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

Add occurrence points to the map

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

map with points

Construct minimum convex hulls from species occurrence data.

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

Create hulls and wrap around dateline

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

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

Plot hull - cuts at 180 but clearly not including all occurrence points

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

incorrect hull

Calculate Area of Hull - must be incorrect based on the hull shape

st_area(hull)

I have also tried applying a pacific centred CRS to the map, points and hulls but suspect that I am applying these either in the wrong order or wrong places? I am very new to using R for spatial analysis so any help hugely apriciated . Thanks.

Answer 1

Please find below a solution to your problem. I used the function st_shift_longitude() from the package sf.

Reprex

• Your data (no changes)

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

• Convert your dataframe "df" into the sf object "species.sf" and shift the longitude with 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))

• Compute the convex hull polygon based on the sf object "species.sf" with group_by(species) (for your general case)

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

• Convert back the sf object "hull" into the dataframe object "hullDF"

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

• Vizualisation of the final result

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)

• Compute area of the hull polygon (needs the units library to convert the result in squared km)

library(units)

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

hull_area
#> 74714882 [km^2]

^{Created on 2021-11-12 by the reprex package (v2.0.1)}