D_01_Make_change_maps.R

D - 01: Make change maps

Libraries

Code

suppressPackageStartupMessages({

library(sf)
library(dplyr)
library(here)
library(skimr)
library(ggplot2)
library(purrr)
  
})

Warning: Paket 'sf' wurde unter R Version 4.4.3 erstellt

Warning: Paket 'ggplot2' wurde unter R Version 4.4.3 erstellt

Warning: Paket 'purrr' wurde unter R Version 4.4.3 erstellt

Part A: Create spatial change data

Load data

Code

grids <-
  readRDS(here("Data/input/grid.rds")) %>%
  filter(scalingID == 1) %>%
  select(datasetID, siteID, verbatimFootprintSRS, croppedArea, geometry)

dat <-
  readRDS(here("Data/output/1_data_filtered.rds")) %>%
  filter(scalingID == 1) %>%
  filter(cell_sampling_repeats == 2 & sp_sampling_repeats == 2)

Merge grids & data

Code

data_sf <-
  full_join(grids, dat)

Extract all sites

Code

all_cells <-
  data_sf %>%
  st_drop_geometry() %>%
  distinct(
    datasetID,
    siteID,
    cell_sampling_repeats
  ) %>%
  mutate(cell_sampling_repeats = as.factor(case_when(
    is.na(cell_sampling_repeats) ~ 0,
    .default = cell_sampling_repeats
  )))

Split data by groups into a list

Code

list_dat <-
  dat %>%
  group_by(datasetID) %>%
  group_split()

Function to get spatial change in occupied sites for 2 atlas periods

Code

get_occupancy_change <- function(data_i, all_cells) {

  # dataset identificator
  atlas_id <- unique(data_i$datasetID)

  #get species list for this atlas
  sp_list <- na.omit(unique(data_i$verbatimIdentification))

  #filter all cells to atlas
  atlas_cells <- all_cells %>% filter(datasetID == atlas_id)
  
  #Process each species in parallel using map
  map_dfr(sp_list, function(this_species) {
    
    #subset data to species
    data_species <- filter(data_i, verbatimIdentification == this_species)
    
    #get occpied sites from period 1 and 2
    sites1 <- unique(data_species %>% filter(samplingPeriodID == 1) %>% pull(siteID))
    sites2 <- unique(data_species %>% filter(samplingPeriodID == 2) %>% pull(siteID))

    
    #get occupancy status for all cells
    stable <- base::intersect(sites1, sites2)
    colonized <- setdiff(sites2, sites1)
    extirpated <- setdiff(sites1, sites2)

    
    #get not sampled cells
    not_sampled <- atlas_cells %>%
      filter(cell_sampling_repeats == 0) %>%
      pull(siteID) %>%
      unique()
    
    #get cells not occupied
    not_occupied <- setdiff(atlas_cells$siteID, c(stable, colonized, extirpated, not_sampled))

    
    #bind to dataframe
    atlas_cells %>%
      mutate(
        verbatimIdentification = this_species,
        change = factor(case_when(
          siteID %in% stable ~ "stable occupancy",
          siteID %in% colonized ~ "newly colonized",
          siteID %in% extirpated ~ "newly extirpated",
          siteID %in% not_sampled ~ "not sampled",
          siteID %in% not_occupied ~ "not occupied",
          TRUE ~ NA_character_
        ))
      )
  })
}

Apply function to each dataset

Code

list_change_species <-
  map(list_dat, ~ get_occupancy_change(.x, all_cells))

make sf data with change info

Code

change_sf <-
  list_change_species %>%
  bind_rows() %>%
  full_join(grids) %>%
  st_as_sf()

check NAs

Code

change_sf %>%
  is.na() %>%
  colSums()

split change data by dataset into lists

Code

list_change_sf <-
  change_sf %>%
  group_by(datasetID) %>%
  group_split()

Get country boundaries for plotting:

Code

library(rnaturalearth)
sf_use_s2(TRUE)

Europe:

Code

my_crs_eu <-
  grids %>%
  st_drop_geometry() %>%
  filter(datasetID == 26) %>%
  distinct(verbatimFootprintSRS) %>%
  pull()

crop europe

Code

buffered_bbox <-
  st_buffer(
    grids %>%
      filter(datasetID == 26) %>%
      st_transform(crs = my_crs_eu) %>%
      select(siteID) %>%
      unique(),
    dist = 15
  )

europe <-
  ne_countries(
    continent = "europe",
    scale = 10
  ) %>%
  st_transform(crs = my_crs_eu) %>%
  st_crop(buffered_bbox) %>%
  select(name)

Czechia:

Code

my_crs_cz <-
  grids %>%
  st_drop_geometry() %>%
  filter(datasetID == 5) %>%
  distinct(verbatimFootprintSRS) %>%
  pull()

czechia <-
  ne_countries(country = "czechia", scale = 10) %>%
  st_transform(crs = my_crs_cz) %>%
  select(name)

Japan:

Code

my_crs_jp <-
  grids %>%
  st_drop_geometry() %>%
  filter(datasetID == 13) %>%
  distinct(verbatimFootprintSRS) %>%
  pull()

japan <-
  ne_countries(country = "japan", scale = 10) %>%
  st_transform(crs = my_crs_jp) %>%
  select(name)

New York state:

Code

my_crs_ny <-
  grids %>%
  st_drop_geometry() %>%
  filter(datasetID == 6) %>%
  distinct(verbatimFootprintSRS) %>%
  pull()

new_york <-
  ne_states(country = "United States of America") %>%
  filter(name == "New York") %>%
  st_transform(crs = my_crs_ny) %>%
  select(name)

Part B: Plotting maps

Transform change data to local CRS

Code

list_country_borders <-
  list(czechia, new_york, japan, europe)

list_crs <-
  list(
    my_crs_cz,
    my_crs_ny,
    my_crs_jp,
    my_crs_eu
  )

list_change_sf_transformed <-
  map2(list_change_sf, list_crs, ~ st_transform(.x, .y))

Colors:

Code

colors <-
  setNames(
    c(
      "#574f7d", "#998ec3", "#f1a340",
      "#F7F7F7", "lightgrey"
    ),
    c(
      "stable occupancy", "newly colonized", "newly extirpated",
      "not occupied", "not sampled"
    )
  )

Theme:

Code

my_theme <-
  ggthemes::theme_map() + # Use theme_map() as a base
  theme(
    legend.position = "right", # Move legend outside
    legend.box = "vertical", # Arrange legend items vertically
    legend.margin = margin(10, 10, 10, 10), # Add space around legend
    plot.margin = margin(10, 50, 10, 10) # Expand right margin to fit legend
  )

Map:

Code

ggplot() +
  geom_sf(
    data = list_change_sf[[1]] %>% st_transform(crs = list_crs[[1]]) %>%
      filter(verbatimIdentification == "Cinclus cinclus"),
    aes(fill = change),
    color = "lightgrey"
  ) +
  geom_sf(
    data = list_country_borders[[1]],
    fill = NA,
    col = "black"
  ) +
  scale_fill_manual(values = colors) +
  my_theme

Map all species via loop

Code

out_paths <- list(
  here("Figures/D_maps/czechia//"),
  here("Figures/D_maps/new_york//"),
  here("Figures/D_maps/japan//"),
  here("Figures/D_maps/europe//")
)

for (atlas_i in seq_along(list_change_sf)) {
  this_atlas <- list_change_sf[[atlas_i]]
  atlas_id <- unique(this_atlas$datasetID)
  sp_list <- unique(this_atlas$verbatimIdentification)
  local_crs <- list_crs[[atlas_i]]

  for (sp_i in seq_along(sp_list)) {
    this_species <- sp_list[sp_i]
    data_species <- this_atlas %>%
      filter(verbatimIdentification == this_species)

    species_expr <- bquote(italic(.(this_species)))
    atlas_expr <- paste0("dataset = ", atlas_id)
    title_expr <- bquote(.(species_expr) * " " * .(atlas_expr))

    map <-
      ggplot() +
      my_theme +
      geom_sf(
        data = data_species %>% st_transform(crs = local_crs),
        aes(fill = change),
        col = "lightgrey"
      ) +
      geom_sf(
        data = list_country_borders[[atlas_i]] %>% st_transform(crs = local_crs),
        fill = NA,
        col = "black"
      ) +
      scale_fill_manual(values = colors)

    # export to ppt
    output_file <- paste0(
      out_paths[[atlas_i]],
      gsub(" ", "_", this_species),
      "_change_map.pptx"
    )

    if (!file.exists(output_file)) { # do not overwrite existing files #
      export::graph2ppt(map,
        width = 9,
        height = 9,
        file = output_file
      )
    }
  }
}



ggpubr::get_legend(map + theme(legend.position = "bottom")) %>% plot()
ggsave(filename = "Change_map_legend.svg", path = here::here("Figures/A_data/"), width = 15)

--- title: "D_01_Make_change_maps.R" format: html: self-contained: true embed-resources: true toc: true # optional: adds a table of contents theme: cosmo # optional: Bootstrap theme code-fold: show # optional: collapsible code blocks code-tools: true # optional: adds copy/paste buttons toc-depth: 4 editor: markdown: wrap: 72 --- # D - 01: Make change maps ```{r} #| echo: false #| message: false #| warning: false #| results: "hide" #| output: false rm(list = ls()) gc() ``` ## Libraries ```{r} suppressPackageStartupMessages({ library(sf) library(dplyr) library(here) library(skimr) library(ggplot2) library(purrr) }) ``` # Part A: Create spatial change data ## Load data ```{r} #| eval: false grids <- readRDS(here("Data/input/grid.rds")) %>% filter(scalingID == 1) %>% select(datasetID, siteID, verbatimFootprintSRS, croppedArea, geometry) dat <- readRDS(here("Data/output/1_data_filtered.rds")) %>% filter(scalingID == 1) %>% filter(cell_sampling_repeats == 2 & sp_sampling_repeats == 2) ``` Merge grids & data ```{r} #| eval: false data_sf <- full_join(grids, dat) ``` Extract all sites ```{r} #| eval: false all_cells <- data_sf %>% st_drop_geometry() %>% distinct( datasetID, siteID, cell_sampling_repeats ) %>% mutate(cell_sampling_repeats = as.factor(case_when( is.na(cell_sampling_repeats) ~ 0, .default = cell_sampling_repeats ))) ``` Split data by groups into a list ```{r} #| eval: false list_dat <- dat %>% group_by(datasetID) %>% group_split() ``` Function to get spatial change in occupied sites for 2 atlas periods ```{r} #| eval: false get_occupancy_change <- function(data_i, all_cells) { # dataset identificator atlas_id <- unique(data_i$datasetID) #get species list for this atlas sp_list <- na.omit(unique(data_i$verbatimIdentification)) #filter all cells to atlas atlas_cells <- all_cells %>% filter(datasetID == atlas_id) #Process each species in parallel using map map_dfr(sp_list, function(this_species) { #subset data to species data_species <- filter(data_i, verbatimIdentification == this_species) #get occpied sites from period 1 and 2 sites1 <- unique(data_species %>% filter(samplingPeriodID == 1) %>% pull(siteID)) sites2 <- unique(data_species %>% filter(samplingPeriodID == 2) %>% pull(siteID)) #get occupancy status for all cells stable <- base::intersect(sites1, sites2) colonized <- setdiff(sites2, sites1) extirpated <- setdiff(sites1, sites2) #get not sampled cells not_sampled <- atlas_cells %>% filter(cell_sampling_repeats == 0) %>% pull(siteID) %>% unique() #get cells not occupied not_occupied <- setdiff(atlas_cells$siteID, c(stable, colonized, extirpated, not_sampled)) #bind to dataframe atlas_cells %>% mutate( verbatimIdentification = this_species, change = factor(case_when( siteID %in% stable ~ "stable occupancy", siteID %in% colonized ~ "newly colonized", siteID %in% extirpated ~ "newly extirpated", siteID %in% not_sampled ~ "not sampled", siteID %in% not_occupied ~ "not occupied", TRUE ~ NA_character_ )) ) }) } ``` Apply function to each dataset ```{r} #| eval: false list_change_species <- map(list_dat, ~ get_occupancy_change(.x, all_cells)) ``` make sf data with change info ```{r} #| eval: false change_sf <- list_change_species %>% bind_rows() %>% full_join(grids) %>% st_as_sf() ``` check NAs ```{r} #| eval: false change_sf %>% is.na() %>% colSums() ``` split change data by dataset into lists ```{r} #| eval: false list_change_sf <- change_sf %>% group_by(datasetID) %>% group_split() ``` Get country boundaries for plotting: ```{r} #| eval: false library(rnaturalearth) sf_use_s2(TRUE) ``` Europe: ```{r} #| eval: false my_crs_eu <- grids %>% st_drop_geometry() %>% filter(datasetID == 26) %>% distinct(verbatimFootprintSRS) %>% pull() ``` crop europe ```{r} #| eval: false buffered_bbox <- st_buffer( grids %>% filter(datasetID == 26) %>% st_transform(crs = my_crs_eu) %>% select(siteID) %>% unique(), dist = 15 ) europe <- ne_countries( continent = "europe", scale = 10 ) %>% st_transform(crs = my_crs_eu) %>% st_crop(buffered_bbox) %>% select(name) ``` Czechia: ```{r} #| eval: false my_crs_cz <- grids %>% st_drop_geometry() %>% filter(datasetID == 5) %>% distinct(verbatimFootprintSRS) %>% pull() czechia <- ne_countries(country = "czechia", scale = 10) %>% st_transform(crs = my_crs_cz) %>% select(name) ``` Japan: ```{r} #| eval: false my_crs_jp <- grids %>% st_drop_geometry() %>% filter(datasetID == 13) %>% distinct(verbatimFootprintSRS) %>% pull() japan <- ne_countries(country = "japan", scale = 10) %>% st_transform(crs = my_crs_jp) %>% select(name) ``` New York state: ```{r} #| eval: false my_crs_ny <- grids %>% st_drop_geometry() %>% filter(datasetID == 6) %>% distinct(verbatimFootprintSRS) %>% pull() new_york <- ne_states(country = "United States of America") %>% filter(name == "New York") %>% st_transform(crs = my_crs_ny) %>% select(name) ``` Part B: Plotting maps Transform change data to local CRS ```{r} #| eval: false list_country_borders <- list(czechia, new_york, japan, europe) list_crs <- list( my_crs_cz, my_crs_ny, my_crs_jp, my_crs_eu ) list_change_sf_transformed <- map2(list_change_sf, list_crs, ~ st_transform(.x, .y)) ``` Colors: ```{r} #| eval: false colors <- setNames( c( "#574f7d", "#998ec3", "#f1a340", "#F7F7F7", "lightgrey" ), c( "stable occupancy", "newly colonized", "newly extirpated", "not occupied", "not sampled" ) ) ``` Theme: ```{r} #| eval: false my_theme <- ggthemes::theme_map() + # Use theme_map() as a base theme( legend.position = "right", # Move legend outside legend.box = "vertical", # Arrange legend items vertically legend.margin = margin(10, 10, 10, 10), # Add space around legend plot.margin = margin(10, 50, 10, 10) # Expand right margin to fit legend ) ``` Map: ```{r} #| eval: false ggplot() + geom_sf( data = list_change_sf[[1]] %>% st_transform(crs = list_crs[[1]]) %>% filter(verbatimIdentification == "Cinclus cinclus"), aes(fill = change), color = "lightgrey" ) + geom_sf( data = list_country_borders[[1]], fill = NA, col = "black" ) + scale_fill_manual(values = colors) + my_theme ``` Map all species via loop ```{r} #| eval: false out_paths <- list( here("Figures/D_maps/czechia//"), here("Figures/D_maps/new_york//"), here("Figures/D_maps/japan//"), here("Figures/D_maps/europe//") ) for (atlas_i in seq_along(list_change_sf)) { this_atlas <- list_change_sf[[atlas_i]] atlas_id <- unique(this_atlas$datasetID) sp_list <- unique(this_atlas$verbatimIdentification) local_crs <- list_crs[[atlas_i]] for (sp_i in seq_along(sp_list)) { this_species <- sp_list[sp_i] data_species <- this_atlas %>% filter(verbatimIdentification == this_species) species_expr <- bquote(italic(.(this_species))) atlas_expr <- paste0("dataset = ", atlas_id) title_expr <- bquote(.(species_expr) * " " * .(atlas_expr)) map <- ggplot() + my_theme + geom_sf( data = data_species %>% st_transform(crs = local_crs), aes(fill = change), col = "lightgrey" ) + geom_sf( data = list_country_borders[[atlas_i]] %>% st_transform(crs = local_crs), fill = NA, col = "black" ) + scale_fill_manual(values = colors) # export to ppt output_file <- paste0( out_paths[[atlas_i]], gsub(" ", "_", this_species), "_change_map.pptx" ) if (!file.exists(output_file)) { # do not overwrite existing files # export::graph2ppt(map, width = 9, height = 9, file = output_file ) } } } ggpubr::get_legend(map + theme(legend.position = "bottom")) %>% plot() ggsave(filename = "Change_map_legend.svg", path = here::here("Figures/A_data/"), width = 15) ```