In HMB3/nenswniche: Rapidly estimate species ranges and intersect with other layers

\

Aim :: The Code below takes the species distribution models (SDMs) for each invertebrate taxa, and the SDMs for their host plant taxa, and intersects the SDMs with a Vegetation layer (), then intersects the combined habitat layer with a fire severity layer.

\

To install, run :

## Set environments 
## java memory limit and temporary raster dir
rm(list = ls())
options(java.parameters = "-Xmx64000m")


## Function to load or install packages
ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg))
    install.packages(new.pkg, dependencies = TRUE, repos="https://cran.csiro.au/")
  sapply(pkg, require,      character.only = TRUE)
}


## Load packages 
library(nenswniche)
data('sdmgen_packages')
ipak(sdmgen_packages)

## Set temporary raster dir for the terra package
terraOptions(memfrac = 0.5, 
             tempdir = 'G:/North_east_NSW_fire_recovery/TEMP')

\

1). Read in Vegetation and Fire Rasters to combine with ivertebrate occurrence data

\

First, get all the rasters to compare with SDM modelling :

• SDM output (Binary habitat suitability, 0-1, 100m) (Source)
• NSW Vegetation (Veg maps, categorical, 100m) (Source )
• Fire layers (FESM, 0-5, 100m) (https://datasets.seed.nsw.gov.au/dataset/fire-extent-and-severity-mapping-fesm)

\

## To Do :
## 1). Check the output with Shawn
## 2). Check all the rasters align in QGIS
## 2). Get examples working for the 18 bugs with host species
source('./R/SDM_GEN_MAXENT_FUNCTIONS.R')
source('./R/SDM_GEN_MAPPING_FUNCTIONS.R')


## Set the projections
sp_epsg54009 <- "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0"
sp_epsg3577  <- "+proj=aea +lat_0=0 +lon_0=132 +lat_1=-18 +lat_2=-36 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"


## SVTM Rasters : https://www.environment.nsw.gov.au/vegetation/state-vegetation-type-map.htm
## The State Vegetation Type Map (SVTM)  of Plant Community Types across NSW, originally @ 5m resolution, re-sampled to 100m
GBAM        <- raster('./data/Remote_sensing/aligned_rasters/GBAM_100m.tif')


## SDM output, re-sampled to 100m
study_sdm_binary <- stack(
  list.files('./output/invert_maxent/Habitat_suitability/SDM_thresholds',
             'current_suit_not_novel_above', full.names = TRUE))


## FESM   : https://datasets.seed.nsw.gov.au/dataset/fire-extent-and-severity-mapping-fesm
## VALUES : 1-4, burn intensity from 2019-2020 fires, originally @ 10m resolution, re-sampled to 100m
FESM_2_100m     <- raster('./data/Remote_sensing/aligned_rasters/FESM_2_100m_align.tif')
FESM_100m       <- raster('./data/Remote_sensing/aligned_rasters/FESM_100m_align.tif')
FESM_100m_align <- readRDS('./data/Remote_sensing/aligned_rasters/FESM_100m_align.rds')


## Read in the SDM data, to intersect with the Veg layers
SVTM_Veg_Class_GDA        <- readRDS('./data/Remote_sensing/aligned_rasters/SVTM_Veg_Class_GDA.rds')
SDM.SPAT.OCC.BG.GDA       <- readRDS('./output/results/SDM_SPAT_OCC_BG_GDA_ALL_TARGET_INVERT_TAXA.rds')
SDM.PLANT.SPAT.OCC.BG.GDA <- readRDS('./output/results/SDM_SPAT_OCC_BG_TARGET_HOST_PLANTS.rds')


## Check projections and resolutions
projection(FESM_100m);projection(study_sdm_binary[[1]]);projection(GBAM);projection(SDM.SPAT.OCC.BG.GDA)
raster::xres(FESM_100m);raster::xres(study_sdm_binary[[1]]);raster::xres(GBAM)

\

2). Select pixels from the Vegetation Habitat Raster that intersect the records for each Invertebrate taxa, as habitat surrogates

\

Let's intersect the SDMs for each invertebrate taxa with the Vegetation Habitat Raster

\

## These still need more exception handling, some taxa fail...why do they fail?
## Are they failing outside the loop, but not inside?


## Select the Vegetation pixels that intersect with the records of each invertebrate species
taxa_records_habitat_intersect(analysis_df    = SDM.SPAT.OCC.BG.GDA,
                               taxa_list      = target.insect.spp,
                               taxa_level     = 'species',
                               habitat_poly   = SVTM_Veg_Class_GDA,
                               output_path    = './output/invert_maxent/Habitat_suitability/SVTM_intersect/',
                               buffer         = 5000)


## Select the Vegetation pixels that intersect with the records of each invertebrate genus 
taxa_records_habitat_intersect(analysis_df    = SDM.SPAT.OCC.BG.GDA,
                               taxa_list      = target.insect.genera,
                               taxa_level     = 'genus',
                               habitat_poly   = SVTM_Veg_Class_GDA,
                               output_path    = './output/invert_maxent/Habitat_suitability/SVTM_intersect/',
                               buffer         = 5000)


## Select the Vegetation pixels that intersect with the records of each invertebrate family
taxa_records_habitat_intersect(analysis_df    = SDM.SPAT.OCC.BG.GDA,
                               taxa_list      = target.insect.families,
                               taxa_level     = 'family',
                               habitat_poly   = SVTM_Veg_Class_GDA,
                               output_path    = './output/invert_maxent/Habitat_suitability/SVTM_intersect/',
                               buffer         = 5000)

\

3). Calculate Estimated fire habitat after the 2019-2020 Fires

For each Invertebrate species, calculate the % of suitable habitat that was burnt by the 2019-2020 fires. We can do this by combining pixels in the rasters like this:

\

• [Invert_SDM + Host_plant_SDM + Inv Veg pixels] * Fire_layer

\

This will give us the % of suitable habitat in each burn intensity category(0-5).

\

## Calculate Insect habitat - fails after this species?
## Code is stalling before or after :: Naranjakotta - it should be the taxa either side of that...
calculate_taxa_habitat(taxa_list          = rev(MAXENT.RESULTS.HOSTS$searchTaxon),
                       targ_maxent_table  = MAXENT.RESULTS.HOSTS,
                       host_maxent_table  = PLANT.MAXENT.RESULTS,
                       target_path        = './output/invert_maxent/back_sel_models/',
                       intersect_path     = 'G:/North_east_NSW_fire_recovery/output/invert_maxent/Habitat_suitability/SVTM_intersect',
                       raster_pattern     = '_SVTM_intersection_5000m.tif',
                       fire_raster        = FESM_100m_align,
                       cell_size          = 100,
                       output_path        = './output/invert_maxent/Habitat_suitability/FESM_SDM_intersect/',
                       country_shp        = 'AUS',
                       country_prj        = CRS("+init=EPSG:3577"),
                       write_rasters      = TRUE)

\

For each taxa, we create a table of the area in square kilometers of suitable habitat that intersects with each burn intensity category from the FESM fire intensity layer. Let's combine all those tables together, to create a master table of estimated burnt area.

## Calculate Insect habitat
INVERT.FESM.list <- list.files('./output/invert_maxent/Habitat_suitability/FESM_SDM_intersect/', 
                               pattern     = '_intersect_Fire.csv', 
                               full.names  = TRUE, 
                               recursive   = TRUE) 


## Now combine the SUA tables for each species into one table 
INVERT.FESM.TABLE <- INVERT.FESM.list %>%

  ## pipe the list into lapply
  lapply(function(x) {

    ## create the character string
    f <- paste0(x)

    ## load each .RData file
    d <- read.csv(f)
    d

  }) %>%

  ## finally, bind all the rows together
  bind_rows


## Subset to just the analysis species - some species did not process properly?
INVERT.FESM.TABLE <-  INVERT.FESM.TABLE[INVERT.FESM.TABLE$Habitat_taxa %in% 
                                          sort(unique(MAXENT.RESULTS.HOSTS$searchTaxon)) , ] %>%  
  .[complete.cases(.), ]


## How many taxa have been processed?
length(unique(INVERT.FESM.TABLE$Habitat_taxa))
View(INVERT.FESM.TABLE)


## Save the FESM intersect results to file
write_csv(INVERT.FESM.TABLE, 
          './output/invert_maxent/Habitat_suitability/FESM_SDM_intersect/INVERT_TAXA_SDM_VEG_intersect_Fire.csv')

\

Figure 1. Example of a combined map of change in climatic suitability from current conditions to 2070. Species occurrence points are plotted in red on the left panel. The cells in the right and bottom panels are coded as either lost (orange cells - present now but not in 2070 according to 4 or more GCMs), gained (green cells - absent now, but present in 2070), stable (blue cells - present now and in 2070), or never suitable (white cells - never present).

\

HMB3/nenswniche documentation built on Jan. 31, 2023, 11:46 p.m.