library(tidyverse)
library(Symobio)
library(sf)
library(rnaturalearth)
library(Rfast)

sf::sf_use_s2(use_s2 = FALSE)

# --------------------------------------#
#       Process range maps           ####
# --------------------------------------#
# Load range maps
range_maps = st_read(
  "~/share/groups/mas_data/Saved_Data_Dropbox_Business/Datensätze/Range Maps/IUCN_range_maps_mammals_version2016/TERRESTRIAL_MAMMALS.shp",
  geometry_column = "geometry",
  promote_to_multi = T
) %>% 
  dplyr::filter(!legend %in% c("Extinct", "Not Mapped"))

# Find species that intersect with Amazon
amazon = st_read(
  "data/geospatial/amazonia_polygons/amazonia_polygons.shp",
  geometry_column = "geometry",
  promote_to_multi = T
) 
range_maps = st_filter(range_maps, amazon)

# Match names against GBIF backbone
range_maps_names_matched = lapply(unique(range_maps$binomial), function(name){
  tryCatch({
    match_result = Symobio::gbif_match_name(name = name)
    if(match_result$status != "ACCEPTED"){
      match_result = gbif_match_name(usageKey = match_result$acceptedUsageKey)
    }
    
    name_matched = if("species" %in% names(match_result)) match_result$species else NA
    data.frame(name_orig = name, name_matched = name_matched)
  }, error = function(e){
    return(NULL)
  })
}) %>% 
  bind_rows()

save(range_maps_names_matched, file = "data/r_objects/range_maps_names_matched.RData")

# Subset range maps to target species and focal region
range_maps = range_maps %>% 
  inner_join(range_maps_names_matched, by = c("binomial" = "name_orig")) %>% 
  group_by(name_matched) %>% 
  summarize(geometry = suppressMessages(st_union(geometry))) 

save(range_maps, file = "data/r_objects/range_maps.RData")

# -------------------------------------------#
#       Process gridded range maps        ####
# -------------------------------------------#
range_maps_gridded = rnaturalearth::ne_countries() %>% 
  dplyr::filter(continent == "South America") %>% 
  sf::st_union() %>% 
  st_make_grid(square = FALSE, cellsize = 1) %>% 
  st_sf() %>% 
  st_join(range_maps, st_intersects, left = F) %>% 
  na.omit()

save(range_maps_gridded, file = "data/r_objects/range_maps_gridded.RData")

# ----------------------------------------------#
#       Calculate range dissimilarity        ####
# ----------------------------------------------#
load("data/r_objects/range_maps_gridded.RData")

range_maps_gridded_id = range_maps_gridded %>% 
  dplyr::group_by(geometry) %>%
  dplyr::mutate(geom_id = cur_group_id()) %>% 
  ungroup()

geometries_unique = range_maps_gridded_id %>% 
  dplyr::select(-name_matched) %>% 
  group_by(geom_id) %>% 
  slice_head(n = 1)

geom_dist = sf::st_distance(geometries_unique, geometries_unique)  %>%  # Takes ~ 10 mins
  as.matrix()

range_maps_split = range_maps_gridded_id %>% 
  group_by(name_matched) %>% 
  group_split()

mean_minimum_distances = lapply(range_maps_split, function(df1){  # Takes ~ 30 mins
  df_dist = lapply(range_maps_split, function(df2){
    dists_subset = geom_dist[df1$geom_id, df2$geom_id, drop = F]
    dist = mean(Rfast::rowMins(dists_subset, value = T)) # Mean minimum distance
    
    df_result = data.frame(
      species1 = df1$name_matched[1],
      species2 = df2$name_matched[1],
      dist = dist
    )
    
    return(df_result)
  }) %>% bind_rows()
  
  return(df_dist)
}) %>% bind_rows()

range_dist = mean_minimum_distances %>% 
  pivot_wider(names_from = species2, values_from = dist) %>% 
  column_to_rownames("species1") %>% 
  as.matrix()

range_dist = range_dist / max(range_dist) # Scale to (0,1)

save(range_dist, file = "data/r_objects/range_dist.RData")