Modeling the Species Distribution of Sea Urchins in the Western Atlantic
+
+## Introduction
+
+This project assessed how climate change will impact the distribution of
+three key species of sea urchins in the Western Atlantic. This was done
+through SDMs produced using LGCP models under a Bayesian framework. For
+more details, reader is referred to the publication “A thorny future for
+sea urchins in the Tropical Western Atlantic” (under review). All codes
+and data used in the study are available in this folder, except original
+climatic layers, which are available online.
+
+## Interactive results
+
+An interactive version of the results of this work is availabe at
+[silasprincipe.github.io/herbinv](http://silasprincipe.github.io/herbinv)
+
+## Files included
+
+**data**: all data used in the study
+– **env**: environmental data
+– **lyva**: *Lytechinus variegatus* data
+– **eclu**: *Echinometra lucunter* data
+– **trve**: *Tripneustes ventricosus* data
+– **sst_limits**: SST data for the thermal limits part **codes**: all
+codes used in the study
+**functions**: all functions that were used either directly or as
+support in the study **gis**: shapefiles
+
+### Codes
+
+- data\_\*: data download and preparing
+- prep\_\*: prepare data (environmental variables)
+- lgcp_prepare_data: prepare all data necessary for fitting the LGCP
+ model through INLA
+- lgcp_model\_\*: models for each species
+- plot\_\*: produce figures
+- sst\_\*: codes used for producing the mechanistic model or for
+ estimating thermal limits
+- supp\_\*: supplementary analysis
+
+## Related repositories
+
+- Species Distribution Modeling of reef building coral species of the
+ Western Atlantic:
+ [silasprincipe/reefbuilders](https://github.com/silasprincipe/reefbuilders)
+
+------------------------------------------------------------------------
+
+Please, contact us if you have any questions on how to use the codes.
diff --git a/README.qmd b/README.qmd
new file mode 100644
index 0000000..f6d36bc
--- /dev/null
+++ b/README.qmd
@@ -0,0 +1,54 @@
+---
+format: gfm
+---
+
+
+
+```{r, echo = FALSE}
+knitr::opts_chunk$set(
+ collapse = TRUE,
+ comment = "#>",
+ fig.path = "figures/"
+)
+```
+
+#
+
+
+
+
+
+
+
+
+
+Echinometra lucunter
+
+
+
+Two species of Echinometra occur in the western Atlantic: E. viridis and E. lucunter. Echinometra lucunter is one of the most common sea urchins on Brazilian reefs, where it can attain great abundance in shallower areas. Although its importance in herbivory can be variable, species of this genus participate actively in the reef carbonate cycle.
+
+
+
+Alvaro E. Migotto. Black sea urchin. Cifonauta image database. Available at: http://cifonauta.cebimar.usp.br/media/4032/ Accessed: 2024-08-12. +
+Current and future distribution (SDMs)
+Both E. lucunter and T. ventricosus showed areas of higher suitability in the Caribbean and south of the Gulf of Mexico, especially along the coast of Campeche to Quintana Roo (Mexico). Both species also show areas of high suitability along the Antilles. Echinometra lucunter has a higher range of distribution to the south.
+
+
+
+
+
+Code
+suppressPackageStartupMessages(library(terra))
+suppressPackageStartupMessages(library(sf))
+library(leaflet)
+library(leaflet.providers)
+library(leafem)
+
+sp <- "eclu"
+
+basedir <- paste0("../results/", sp, "/predictions/")
+
+sdm_proj <- list.files(basedir)
+sdm_proj <- sdm_proj[grepl("mean", sdm_proj)]
+sdm_proj_cont <- sdm_proj[grepl("cont", sdm_proj)]
+
+proj_lays <- rast(paste0(basedir, sdm_proj_cont))
+proj_lays <- project(proj_lays, "EPSG:3857")
+
+# Normalize to 0-1
+proj_lays <- (proj_lays - min(terra::minmax(proj_lays$eclu_mean_m4_cont_current)[1,])) / (terra::minmax(proj_lays$eclu_mean_m4_cont_current)[2,] - terra::minmax(proj_lays$eclu_mean_m4_cont_current)[1,])
+
+# Get areas of extrapolation
+extrap_lays <- proj_lays[[2:4]]
+extrap_lays[extrap_lays <= terra::minmax(proj_lays$eclu_mean_m4_cont_current)[2,]] <- NA
+extrap_lays[!is.na(extrap_lays)] <- 1
+
+extrap_shape <- lapply(1:3, function(id){
+ terra::project(terra::as.polygons(extrap_lays[[id]]), "EPSG:4326")
+})
+
+# Set maximum to the maximum of current layer
+proj_lays[proj_lays > 1] <- 1
+
+# Load points
+pts <- read.csv(paste0("../data/", sp, "/", sp, "_filt.csv"))
+pts <- vect(pts, geom = c("x", "y"), crs = crs(rast(paste0(basedir, sdm_proj_cont[1]))))
+pts <- project(pts, "EPSG:4326")
+pts <- as.data.frame(geom(pts))
+
+# Plot maps
+leaflet() %>%
+ #addProviderTiles("OpenStreetMap.Mapnik", group = "OSM") %>%
+ addProviderTiles("Esri.WorldGrayCanvas", group = "ESRI Gray") %>%
+ addRasterImage(
+ proj_lays[[1]],
+ project = F,
+ colors = colorNumeric(
+ palette = rev(c("#A84C00", "#D97D27", "#F5BD44", "#FFD561", "#FFF291",
+ "#FFFFBF", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4", "#313695")),
+ domain = c(0,1),
+ na.color = "#00000000",
+ alpha = FALSE,
+ reverse = FALSE
+ ),
+ group = "Current"
+ ) %>%
+ addRasterImage(
+ proj_lays[[2]],
+ project = F,
+ colors = colorNumeric(
+ palette = rev(c("#A84C00", "#D97D27", "#F5BD44", "#FFD561", "#FFF291",
+ "#FFFFBF", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4", "#313695")),
+ domain = c(0,1),
+ na.color = "#00000000",
+ alpha = FALSE,
+ reverse = FALSE
+ ),
+ group = "SSP1 (RCP2.6)"
+ ) %>%
+ addPolygons(data = extrap_shape[[1]], group = "SSP1 (RCP2.6)") %>%
+ addRasterImage(
+ proj_lays[[3]],
+ project = F,
+ colors = colorNumeric(
+ palette = rev(c("#A84C00", "#D97D27", "#F5BD44", "#FFD561", "#FFF291",
+ "#FFFFBF", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4", "#313695")),
+ domain = c(0,1),
+ na.color = "#00000000",
+ alpha = FALSE,
+ reverse = FALSE
+ ),
+ group = "SSP2 (RCP4.5)"
+ ) %>%
+ addPolygons(data = extrap_shape[[2]], group = "SSP2 (RCP4.5)") %>%
+ addRasterImage(
+ proj_lays[[4]],
+ project = F,
+ colors = colorNumeric(
+ palette = rev(c("#A84C00", "#D97D27", "#F5BD44", "#FFD561", "#FFF291",
+ "#FFFFBF", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4", "#313695")),
+ domain = c(0,1),
+ na.color = "#00000000",
+ alpha = FALSE,
+ reverse = FALSE
+ ),
+ group = "SSP3 (RCP7.0)"
+ ) %>%
+ addPolygons(data = extrap_shape[[3]], group = "SSP3 (RCP7.0)") %>%
+ addLegend(pal = colorNumeric(
+ palette = c("#A84C00", "#D97D27", "#F5BD44", "#FFD561", "#FFF291",
+ "#FFFFBF", "#E0F3F8", "#ABD9E9", "#74ADD1", "#4575B4", "#313695"),
+ domain = c(0,1),
+ na.color = "#00000000",
+ alpha = FALSE,
+ reverse = FALSE
+ ), values = values(proj_lays[[1]]), title = "ROR", opacity = 1, position = "bottomright",
+ labFormat = labelFormat(transform = function(x) sort(x, decreasing = TRUE))) %>%
+ addCircleMarkers(lng = pts$x, lat = pts$y,
+ #clusterOptions = markerClusterOptions(),
+ radius = 5, weight = 2.5,
+ group = "Occurrence") %>%
+ addLayersControl(
+ baseGroups = c("Current", "SSP1 (RCP2.6)", "SSP2 (RCP4.5)", "SSP3 (RCP7.0)"),
+ overlayGroups = c("Occurrence"),
+ options = layersControlOptions(collapsed = FALSE)
+ ) %>%
+ setView(-60, 0, zoom=3)
+
+
+
+Changes in future distribution (SDMs)
+Echinometra lucunter do not present any apparent loss in its distribution range compared to the current scenario. This species would increase its range of suitable areas to the north and to the south.
+
+
+
+
+
+Code
+delta <- proj_lays[[2:4]] - proj_lays[[1]]
+
+# Plot maps
+leaflet() %>%
+ #addProviderTiles("OpenStreetMap.Mapnik", group = "OSM") %>%
+ addProviderTiles("Esri.WorldGrayCanvas", group = "ESRI Gray") %>%
+ addRasterImage(
+ delta[[1]],
+ project = F,
+ colors = colorNumeric(
+ palette = "BrBG",
+ domain = c(-1,1),
+ na.color = "#00000000",
+ alpha = FALSE,
+ reverse = FALSE
+ ),
+ group = "SSP1 (RCP2.6)"
+ ) %>%
+ addRasterImage(
+ delta[[2]],
+ project = F,
+ colors = colorNumeric(
+ palette = "BrBG",
+ domain = c(-1,1),
+ na.color = "#00000000",
+ alpha = FALSE,
+ reverse = FALSE
+ ),
+ group = "SSP2 (RCP4.5)"
+ ) %>%
+ addRasterImage(
+ delta[[3]],
+ project = F,
+ colors = colorNumeric(
+ palette = "BrBG",
+ domain = c(-1,1),
+ na.color = "#00000000",
+ alpha = FALSE,
+ reverse = FALSE
+ ),
+ group = "SSP3 (RCP7.0)"
+ ) %>%
+ addLegend(pal = colorNumeric(
+ palette = "BrBG",
+ domain = c(-1,1),
+ na.color = "#00000000",
+ alpha = FALSE,
+ reverse = T
+ ), values = seq(-1, 1, by = 0.1), title = "Delta ROR", opacity = 1, position = "bottomright",
+ labFormat = labelFormat(transform = function(x) sort(x, decreasing = TRUE))) %>%
+ addCircleMarkers(lng = pts$x, lat = pts$y,
+ #clusterOptions = markerClusterOptions(),
+ radius = 5, weight = 2.5,
+ group = "Occurrence") %>%
+ addLayersControl(
+ baseGroups = c("SSP1 (RCP2.6)", "SSP2 (RCP4.5)", "SSP3 (RCP7.0)"),
+ overlayGroups = c("Occurrence"),
+ options = layersControlOptions(collapsed = FALSE)
+ ) %>%
+ setView(-60, 0, zoom=3)
+
+
+
+Current and future distribution (mechanistic model)
+Echinometra lucunter had a similar distribution than L. variegatus, but a smaller extension to the south and unsuitable areas in the northern portion of the Gulf of Mexico. Echinometra lucunter is expected to lose approximately 13, 48, and 65% of its suitable area in the SSP1, SSP2, and SSP3 scenarios.
+
+
+
+
+
+
+
+
+
+Code
+# Load layers and prepare
+# Load threshold data ----
+load("../data/sst_limits/allspecies_oisst_thvalues.RData")
+
+# Load results ----
+sp <- "eclu" # Each species is run separately [try "eclu" and "trve"]
+
+# Load rasters generated before
+curr <- rast(paste0("../data/sst_limits/", sp, "_current_thresh.tif"))
+ssp1 <- rast(paste0("../data/sst_limits/", sp, "_", "ssp126", "_thresh.tif"))
+ssp2 <- rast(paste0("../data/sst_limits/", sp, "_", "ssp245", "_thresh.tif"))
+ssp3 <- rast(paste0("../data/sst_limits/", sp, "_", "ssp370", "_thresh.tif"))
+
+curr <- project(curr, "EPSG:3857")
+ssp1 <- project(ssp1, "EPSG:3857")
+ssp2 <- project(ssp2, "EPSG:3857")
+ssp3 <- project(ssp3, "EPSG:3857")
+
+# Get the percentage of time to use as threshold (mean of min and max point)
+lval <- round(((thresholds[[sp]]$time_inrange_hottest_point +
+ thresholds[[sp]]$time_inrange_coolest_point)/2),
+ 2) # round to 2 digits
+
+
+# Get the polygons of the areas that are suitable
+get.pol <- function(x){
+ # temp <- terra::app(x, function(x){
+ # x[x < lval] <- NA
+ # x[x >= lval] <- 1
+ # x
+ # })
+ # temp <- as.polygons(temp)
+ # temp <- aggregate(buffer(temp,0.0001)) # We use a negligible value here
+ # # to solve problems in the pols
+ # # conversion.
+ # temp <- project(temp, "EPSG:4326")
+ # # temp <- st_as_sf(temp)
+ # # temp <- st_set_crs(temp, crs("EPSG:4326"))
+ # temp
+ x[x < lval] <- NA
+ x[x >= lval] <- 1
+ terra::project(terra::as.polygons(x), "EPSG:4326")
+}
+
+curr.p <- get.pol(curr)
+ssp1.p <- get.pol(ssp1)
+ssp2.p <- get.pol(ssp2)
+ssp3.p <- get.pol(ssp3)
+
+
+
+# Plot maps
+leaflet() %>%
+ #addProviderTiles("OpenStreetMap.Mapnik", group = "OSM") %>%
+ addProviderTiles("Esri.WorldGrayCanvas", group = "ESRI Gray") %>%
+ addRasterImage(
+ curr,
+ project = F,
+ colors = colorNumeric(
+ palette = "Spectral",
+ domain = c(0,1),
+ na.color = "#00000000",
+ alpha = FALSE,
+ reverse = FALSE
+ ),
+ group = "Current"
+ ) %>%
+ addPolygons(data = curr.p,
+ group = "Current Suitable") %>%
+ addRasterImage(
+ ssp1,
+ project = F,
+ colors = colorNumeric(
+ palette = "Spectral",
+ domain = c(0,1),
+ na.color = "#00000000",
+ alpha = FALSE,
+ reverse = FALSE
+ ),
+ group = "SSP1 (RCP2.6)"
+ ) %>%
+ addPolygons(data = ssp1.p,
+ group = "SSP1 Suitable") %>%
+ addRasterImage(
+ ssp2,
+ project = F,
+ colors = colorNumeric(
+ palette = "Spectral",
+ domain = c(0,1),
+ na.color = "#00000000",
+ alpha = FALSE,
+ reverse = FALSE
+ ),
+ group = "SSP2 (RCP4.5)"
+ ) %>%
+ addPolygons(data = ssp2.p,
+ group = "SSP2 Suitable") %>%
+ addRasterImage(
+ ssp3,
+ project = F,
+ colors = colorNumeric(
+ palette = "Spectral",
+ domain = c(0,1),
+ na.color = "#00000000",
+ alpha = FALSE,
+ reverse = FALSE
+ ),
+ group = "SSP3 (RCP7.0)"
+ ) %>%
+ addPolygons(data = ssp3.p,
+ group = "SSP3 Suitable") %>%
+ addLegend(pal = colorNumeric(
+ palette = "Spectral",
+ domain = c(0,1),
+ na.color = "#00000000",
+ alpha = FALSE,
+ reverse = TRUE
+ ), values = seq(0, 1, by = 0.1), title = "% time", opacity = 1, position = "bottomright",
+ labFormat = labelFormat(transform = function(x) sort(x, decreasing = TRUE) * 100)) %>%
+ addCircleMarkers(lng = pts$x, lat = pts$y,
+ #clusterOptions = markerClusterOptions(),
+ radius = 5, weight = 2.5,
+ group = "Occurrence") %>%
+ addLayersControl(
+ baseGroups = c("Current", "SSP1 (RCP2.6)", "SSP2 (RCP4.5)", "SSP3 (RCP7.0)"),
+ overlayGroups = c("Occurrence", "Current Suitable", "SSP1 Suitable", "SSP2 Suitable", "SSP3 Suitable"),
+ options = layersControlOptions(collapsed = FALSE)
+ ) %>%
+ hideGroup(c("Occurrence", "Current Suitable", "SSP1 Suitable", "SSP2 Suitable", "SSP3 Suitable")) %>%
+ setView(-60, 0, zoom=3)
+
+
+
+