-
Notifications
You must be signed in to change notification settings - Fork 1
/
03_change_suitability.R
145 lines (121 loc) · 4.98 KB
/
03_change_suitability.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
####################################################
###### Calculate change in suitability
# .........................................
# .........................................
# Packages ####
library("tidyverse")
library("magrittr")
library("svglite")
library("gridExtra")
library("sp")
library("dismo")
library("raster")
library("rgdal")
library("rgeos")
# .........................................
# .........................................
# Read species names
setwd("C:/Users/Pablo/OneDrive - Universidad de Córdoba/1_proyectos/2019_CocoAgroForecast/WP5 - suitability")
BD_calib <- read.csv("./BD/BD_calibrate.csv")
sp <- unique(BD_calib$species)
spnames <- sort(sp) # hacer acrónimos
# .........................................
# .........................................
# Prepare data and files to process maps ####
# names of RCP scenarios
RCP <- c("ssp126","ssp585")
scenario <- c("2021-2040","2041-2060")
#define extention of study area
ext <- raster::extent(-17,15,2,13)
# define projection
proj <- "+proj=longlat +datum=WGS84"
output <- "processing/species_sets/"
dir.create(output, showWarnings = FALSE, recursive = TRUE)
### read shapefile of inland water
# lakes <- "data/shapefiles/water_areas/mesoamerica_water_areas_dcw.shp"
# lakes %<>%
# readOGR(.) %>%
# subset(.$HYC_DESCRI == "Perennial/Permanent") %>%
# raster::crop(. , ext)
#
#
# # Read country borders
# border <- "data/shapefiles/country_borders/Mesoamerica.shp"
# border %<>%
# readOGR(.) %>%
# raster::crop(. , ext)
# .........................................
# .........................................
# Run over spnames ####
# create NULL dataframe to keep data on area changes for all species
changes_all <- NULL
for (i in seq_along(spnames)) {
cat(i, "out of", length(spnames), "\n")
presence <- paste0("processing/enm/", spnames[i], "/ensembles/presence/")
#create directory to save rasters
output2 <- paste0(output, spnames[i])
dir.create(output2, showWarnings = FALSE, recursive = TRUE)
# load rasters of current niches
# raster of current presence-absence
presence_current <- raster::stack(paste0(presence,spnames[i], "__bio_current.grd"))
# crop raster using Mesoamerica extention
presence_current <- raster::crop(presence_current, ext)
# # remove Caribbean islands
# presence_current <- raster::mask(presence_current, border, inverse=FALSE)
# # remove presence in lakes
# presence_current <- raster::mask(presence_current, lakes, inverse=TRUE)
presence_current <- raster::stack(presence_current)
# run over RCP models
for (s in seq_along(scenario)){
for(j in seq_along(RCP)) {
# read presence-absence rasters under modelled RCP scenario##
presence_rcp <- raster::stack(list.files(paste0("processing/enm/", spnames[i] , "/ensembles/presence"),
pattern = paste0(RCP[j],"_",scenario[s],".gri$"),
full.names = TRUE))
# mean values of 1-0 raster (presence and absence)
# this is also the measure agreement raster
presence_rcp_mean <- raster::calc(presence_rcp, fun = mean)
# Define likelihood RCP mask
# more than 66% of the raster must agree on the presence
# of the species in each grid cell
# further information about likelihood at Mastrandrea el al (2010)
# raster for future presence-absence
presence_rcp <- presence_rcp_mean
presence_rcp[presence_rcp[] < 0.659] <- 0
presence_rcp[presence_rcp[] >= 0.659] <- 1
#raster of threshold for future presence (defined by likelihood)
thresh_presence_rcp <- presence_rcp
thresh_presence_rcp[thresh_presence_rcp[] == 0] <- NA
#identify change in suitability in RCP scenario
#future minus current raster
presence_rcp[presence_rcp[]==1] <- 2
#change in suitability codes
# 1 = always suitable
# 0 = never suitable
# -1 = no longer suitable
# 2 = new habitat
change_suit_rcp <- raster::overlay(presence_rcp,
presence_current,
fun = function(x,y) {(x-y)})
# write tif files
# current presence-absence layer
writeRaster(change_suit_rcp,
filename = paste0(output2,"/",spnames[i],"_",RCP[j],"_",scenario[s],"_change.tif"),
format = "GTiff",
overwrite = TRUE)
writeRaster(presence_rcp,
filename = paste0(output2,"/",spnames[i],"_",RCP[j],"_",scenario[s],"_presence.tif"),
format = "GTiff",
overwrite = TRUE)
writeRaster(presence_rcp_mean,
filename = paste0(output2,"/",spnames[i],"_",RCP[j],"_",scenario[s],"_presence_agreement.tif"),
format = "GTiff",
overwrite = TRUE)
}
}
# current presence-absence
writeRaster(presence_current,
filename = paste0(output2,"/",spnames[i],"_presence.tif"),
format="GTiff",
overwrite=TRUE)
}