-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpredict 2050 lc.R
51 lines (46 loc) · 2.67 KB
/
predict 2050 lc.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
library(raster)
library(rgdal)
library(dismo)
library(sp)
library(rJava)
library(readxl)
library(mapview)
library(ggplot2)
library(jsonlite)
require(sf)
library(devtools)
library(xlsx)
library(dplyr)
# options(java.parameters = "-Xmx8000m")
#latlong <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
#SR.ORG8287 <- CRS('+proj=cea +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs')
#extent.asia.sr <- extent(11000000, 13000000, 900000, 2600000)
# Current
list.raster.full <- list.files("F:/Working/2018/PhD_research/enviromental_variables_new/current", full.names = T, pattern = ".tif")
predictors <- stack(list.raster.full)
list.raster.names <- list.files("F:/Working/2018/PhD_research/enviromental_variables_new/current", pattern = ".tif")
list.raster.names <- gsub('.{4}$','',list.raster.names)
# predictors <- stack(list.raster.full)
names(predictors) <- list.raster.names
# names(predictors) <- c('bio1','bio10','bio11','bio12','bio13','bio14','bio15','bio16','bio17','bio18','bio19','bio2','bio3','bio4','bio5','bio6','bio7','bio8','bio9',
# 'cropland','decidous_forest','evergreen_forest','flooded_vegetation','grassland_scrub','karst','urban')
climate.current <- subset(predictors, c("bio2","bio10","bio11","bio12","bio18","bio19",'karst'))
# karst <- raster('F:/Working/2018/PhD_research/enviromental_variables/current/karst.tif')
# 2050
land.cover.2050 <- list.files("F:/Working/2018/PhD_research/enviromental_variables_new/2050/landcover", full.names = T, pattern = ".tif")
names(land.cover.2050) <- c('cropland','decidous_forest','evergreen_forest','flooded_vegetation','grassland_scrub','urban')
land.cover.2050 <- stack(land.cover.2050)
#import species list
species.list <- read_excel('F:/Working/2018/PhD_research/SDM output R/10_7_23_Transboundary/list of species.xlsx')
species.list <- species.list%>%filter(Num_occurence_all > 19)
nrow(species.list)
#predict
for(sp in species.list$Scientific_name){
load(paste0('F:/Working/2018/PhD_research/SDM output R/10_7_23_Transboundary/response/',sp,'.RData'))
maxent.result <- read.csv(paste0('F:/Working/2018/PhD_research/SDM output R/10_7_23_Transboundary/response/',sp,'/maxentResults.csv'), header=T, sep=',')
predictors.2050 <- stack(climate.current, land.cover.2050)
px.2050 <- predict(predictors.2050, xm, ext=extent(11000000, 13000000, 900000, 2600000))
species.2050.bin <- px.2050 > maxent.result$X10.percentile.training.presence.Cloglog.threshold
writeRaster(species.2050.bin, filename = paste0('F:/Working/2018/PhD_research/SDM output R/10_7_23_Transboundary/output_maps/',sp,'_2050_lc.tif'),format='GTiff', overwrite=T)
print(sp)
}