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.

CodePudding user response：

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 <- st_as_sf(df,coords = c("longitude", "latitude"), crs = 4326)
species.sf <- st_shift_longitude(species.sf)


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" and shift the longitude with st_shift_longitude()

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


st_crs(hull) <- 4326
hull<-st_shift_longitude(st_geometry(hull))

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

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

• 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 <- st_area(hull)
set_units(hull_area, km^2)
#> 74714882 [km^2]

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