Lytechinus variegatus

Lytechinus variegatus occurs from North Carolina up to the south of Brazil. It occurs on rocky reefs, sand flats, and seagrass meadows, where it can have great importance in grazing seagrass.

Alvaro E. Migotto. Sea urchin. Cifonauta image database. Available at: http://cifonauta.cebimar.usp.br/media/10635/ Accessed: 2024-04-22. CC BY-NC-SA 3.0 DEED

Current and future distribution (SDMs)

Areas that are suitable for L. variegatus span from the south of Brazil up to North Carolina, which are comparable to the ones for the other species. However, L. variegatus showed areas with a high Relative Occurrence Rate (ROR) in the western Gulf of Mexico that are not shared by neither E. lucunter or T. ventricosus.

Code
suppressPackageStartupMessages(library(terra))
suppressPackageStartupMessages(library(sf))
library(leaflet)
library(leaflet.providers)
library(leafem)

sp <- "lyva"

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$lyva_mean_m4_cont_current)[1,])) / (terra::minmax(proj_lays$lyva_mean_m4_cont_current)[2,] - terra::minmax(proj_lays$lyva_mean_m4_cont_current)[1,])

# Get areas of extrapolation
extrap_lays <- proj_lays[[2:4]]
extrap_lays[extrap_lays <= terra::minmax(proj_lays$lyva_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)

Major changes (including loss in suitable areas) are expected in the distribution of L. variegatus.

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)

Lytechinus variegatus had the widest suitable area according to the mechanistic model. The species would lose approximately 15, 49, and 64% 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 <- "lyva" # 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)