Better cartographic image of climate zones at risk of invasion

[SanderDevisscher]

as proposed by @timadriaens

Having discussed with @SanderDevisscher and looked at the results of the bulbul climate matching, we notice that the maps produced are not very clear nor useful to provide an idea to a risk assessor about areas at risk. The reason is that the %overlap (perc_climate) is calculated on the entire global dataset and many records do not "hit" a climate classification that occurs in the risk assessment area (Europe). Consequently, when making the legend, the 'suitability' is extremely low overall.

To better visualize the areas at risk in Europe, we propose to calculate the % overlap using only the records that intersect in a climate classification occurring in Europe. Maybe we could also think about an automated way of reclassification (highly suitable, moderately suitable, low suitability).

=> I'll add this "subset" as a new type of map output based on the code below, the cm output and the single species maps
https://github.com/inbo/riparias-prep/blob/5c3039199c6104a30adf0a382651e16fe2c202b3/scripts/bul_bul_cm.Rmd#L89-L156

this output has to be NULL when maps == FALSE and/or when missing(region)

n_totaal, and consequently the perc_climate, needs to be recalculated to only take observations in climate zones present in the region of intrest into account.

the new outputs shall be:

• region_cm: the recalculated cm
• region_single_species_maps: The visualized region_cm with bbox around the region_shape

[SanderDevisscher]

Not commited since not working as wished. Will look at it later.

 ## Region ##
  if (missing(region)) {
    # no region is provided => worldwide
    message("region not provided, using worldmap")
    region_shape <- getMap(resolution = "low")
  }else{
    if (is.character(region)) {
@ -346,7 +347,8 @@ climate_match <- function(region,
           .data$perc_climate)
  
  # Determine future scenarios ####
  scenarios <- c("2001-2025-A1FI",
  scenarios <- c("1980-2016",
                 "2001-2025-A1FI",
                 "2026-2050-A1FI",
                 "2051-2075-A1FI",
                 "2071-2100_Beck",
@ -360,7 +362,11 @@ climate_match <- function(region,
  
  # Calculate KG codes 
  for(s in scenarios){
    shape <- future[[s]]
    if(s == "1980-2016"){
      shape <- observed[[s]]
    }else{
      shape <- future[[s]]
    }
    
    if(c("gridcode") %in% colnames(shape@data)){
      shape@data <- shape@data %>% 
@ -380,7 +386,8 @@ climate_match <- function(region,
  }
  
  output_1 <- output %>% 
    dplyr::filter(grepl(pattern = "Beck", .data$scenario)) %>% 
    dplyr::filter(grepl(pattern = "Beck", .data$scenario) 
                  | scenario == "1980-2016") %>% 
    left_join(KG_Beck, by = c("KG_GridCode" = "GRIDCODE"))
  
  output_2 <- output %>% 
@ -411,6 +418,17 @@ climate_match <- function(region,
    }
  }
  
  # Calculate region threshold parameters ####
  region_cm <- cm %>% 
    group_by(taxonKey, 
             acceptedScientificName, 
             scenario) %>% 
    mutate(n_region = sum(n_climate, na.rm = TRUE),
           region_prec_climate = n_climate/n_region) %>% 
    ungroup() %>% 
    select(-n_totaal, -perc_climate)
  
  
  # Thresholds ####
  if(missing(n_limit)){
    warning("no n_totaal threshold was provided. defaults to 0!")
@ -623,22 +641,23 @@ climate_match <- function(region,
    single_species_map <- leaflet(sea) %>% 
      addPolygons(data = sea,
                  fillColor = "#e0e0e0",
                  weight = 0.5) %>% 
      addLegend(colors = "black",
                labels = "observations",
                position = "bottomleft")
                  weight = 0.5)
    
    # Create single species maps
    
    ## Create scenarios list (future + current)
    scenarios_2 <- c("1980-2016", scenarios)
    
    ## Prepare output list
    single_species_maps <- list_along(taxonkey)
    names(single_species_maps) <- taxonkey
    
    ## Merge unfiltered climate match results with climate shapes
    for (i in 1:length(taxonkey)) {
      
      t <- taxonkey[i]
      
      # Subset data per species
      temp_data <- data_overlay_unfiltered %>% 
        dplyr::filter(taxonKey == t) %>% 
        select(-Description)
@ -652,6 +671,7 @@ climate_match <- function(region,
        
        for(s in scenarios_2){
          
          # get scenario shape
          if(s == "1980-2016"){
            scenario_shape <- observed[[s]]
          }else{
@ -669,6 +689,7 @@ climate_match <- function(region,
              left_join(legends$KG_A1FI, by = "GRIDCODE")
          }
          
          # Combine data with shape
          temp_climate <- sp::merge(scenario_shape, temp_data, 
                                    by = "Classification",
                                    all.y = TRUE,
@ -713,35 +734,156 @@ climate_match <- function(region,
                    weight = 0.5,
                    group = ~scenario,
                    popup = ~popup) %>% 
        addCircleMarkers(data = data_sp_species_obs,
                         color = "black",
                         radius = 1) %>% 
        addLegend(pal = pal_current,
                  values = seq(from = 0, 
                               to = 1, 
                               by = 0.1),
                  position = "bottomleft",
                  title = "<strong>Climate match</strong>") %>% 
        addCircleMarkers(data = data_sp_species_obs,
                         color = "black",
                         radius = 1) %>% 
        addLegend(colors = "black",
                  labels = "observations",
                  position = "bottomleft") %>% 
        addLayersControl(baseGroups = ~temp_shape@data$scenario)
      
      
      single_species_maps[[i]] <- scenario_map
    }  
    } 
    ## region species maps ####
    
    # Create basemap
    region_species_map <- leaflet(sea) %>% 
      addPolygons(data = sea,
                  fillColor = "#e0e0e0",
                  weight = 0.5)
    
    ## Prepare output list
    region_species_maps <- list_along(taxonkey)
    names(region_species_maps) <- taxonkey
    
    ## Merge unfiltered climate match results with climate shapes
    for (i in 1:length(taxonkey)) {
      
      t <- taxonkey[i]
      
      # Extract cm data from region_cm
      temp_data_1 <- region_cm %>% 
        dplyr::filter(taxonKey == t) %>% 
        select(-Description)
      
      species <- unique(temp_data$acceptedScientificName)
      
      if(is_empty(species)){
        next
      }else{
        temp_shape <- data.frame()
        
        # Get needed scenario shapes
        for(s in scenarios_2){
          
          # subset data per scenario
          scenario_data <- temp_data %>% 
            filter(scenario == s)
          
          if(s == "1980-2016"){
            scenario_shape <- observed[[s]]
          }else{
            scenario_shape <- future[[s]]
          }
          if (grepl("Beck", s) | s == "1980-2016") {
            scenario_shape@data <- scenario_shape@data %>% 
              mutate(GRIDCODE = as.double(.data$gridcode),
                     ID = .data$Id) %>% 
              select(-gridcode, -Id) %>% 
              left_join(legends$KG_Beck, by = "GRIDCODE")
          }else{
            scenario_shape@data <- scenario_shape@data %>% 
              mutate(GRIDCODE = as.double(.data$GRIDCODE)) %>% 
              left_join(legends$KG_A1FI, by = "GRIDCODE")
          }
          
          temp_climate <- sp::merge(scenario_shape, scenario_data, 
                                    by = "Classification",
                                    all.y = TRUE,
                                    duplicateGeoms = TRUE)
          
          temp_climate@data <- temp_climate@data %>% 
            mutate(taxonKey = t,
                   acceptedScientificName = species,
                   scenario = s)
          
          if(class(temp_shape) == "data.frame"){
            temp_shape <- temp_climate
          }else{
            temp_shape <- rbind.SpatialPolygonsDataFrame(temp_shape, 
                                                         temp_climate)
          }
        }
        
        
        temp_shape@data <- temp_shape@data %>% 
          mutate(popup = paste0("<strong>Classification: </strong>", .data$Description, " (", .data$Classification, ")", 
                                "</br><strong>ScientificName: </strong>", 
                                .data$acceptedScientificName,
                                "</br><strong>%obs in climate: </strong>", 
                                round(.data$region_prec_climate*100, 2), "%",
                                "</br><strong>scenario: </strong>",
                                .data$scenario))
        
        temp_shape <- subset(temp_shape, !is.na(temp_shape$Classification))
        
        # Subset observations
        data_sp_species_obs <- subset(data_sp, 
                                      data_sp$taxonKey == t)
        
        # Add layer to map
        scenario_map <- region_species_map %>% 
          addPolygons(data = temp_shape,
                      color = "#bababa",
                      fillColor = ~pal_current(.data$region_perc_climate),
                      fillOpacity = 0.8,
                      stroke = TRUE,
                      weight = 0.5,
                      group = ~scenario,
                      popup = ~popup) %>% 
          addLegend(pal = pal_current,
                    values = seq(from = 0, 
                                 to = 1, 
                                 by = 0.1),
                    position = "bottomright",
                    title = "<strong>Climate match</strong>") %>% 
          addCircleMarkers(data = data_sp_species_obs,
                           color = "black",
                           radius = 1) %>% 
          addLegend(colors = "black",
                    labels = "observations",
                    position = "bottomright") %>% 
          addLayersControl(baseGroups = ~temp_shape@data$scenario)
        
        region_species_maps[[i]] <- scenario_map
      }
    }
  }else{
    message("maps are disabled")
    current_climate_map <- NULL
    future_scenario_maps <- NULL
    single_species_maps <- NULL
    region_species_maps <- NULL
  }
  
  
  
  # Return ####
  return(list(unfiltered = data_overlay_unfiltered, 
              cm = cm,
              region_cm = region_cm,
              filtered = data_overlay_scenario_filtered,
              future = future_climate,
              spatial = data_sp,
              current_map = current_climate_map,
              future_maps = future_scenario_maps,
              single_species_maps = single_species_maps))
              single_species_maps = single_species_maps,
              region_species_maps = region_species_maps))
}

[damianooldoni]

@SanderDevisscher: has this already been tackled?

• If not, can you transfer this issue to the new package/repo? To do so, click the "Transfer issue" here on the right side of the page.
• If you don't see the "Transfer issue" button, let me know.
• If the repo has still not be done, could you please add the blocked label to this issue?
• If you don't have access to the labels, please let me know.

What a beautiful "if cascade" 😄

[SanderDevisscher]

@damianooldoni somehow I can't transfer the issue to inbo/ClimateCastR

[damianooldoni]

@SanderDevisscher: you are right. In GitHub documentation, I read:

Note: You can only transfer issues between repositories owned by the same user or organization account.

So, we cannot transfer this issue to ClimateCastR as its GitHub organization is inbo, while this repo is owned by trias-project GitHub organization.