-
Notifications
You must be signed in to change notification settings - Fork 0
/
3_prepStatesGrid.r
110 lines (86 loc) · 3.42 KB
/
3_prepStatesGrid.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
#clean
rm(list=ls())
# load lib
library(doParallel)
library(raster)
# load params
pars <- read.table("./pars/GenSA_rf_all_2_5y_rep1.txt",row.names=1)
# list landscape
clim_files <- list.files("data/futClimGrid/stmClim",full.names=TRUE)
# source res analytics function
source('./lib/resAnalytics.r')
source('./lib/convertStates.r')
# create output dir
system('mkdir -p ./data/futStatesGrid/probs')
system('mkdir -p ./data/futStatesGrid/stm100')
system('mkdir -p ./data/futStatesGrid/stm1000')
# open cluster
cl <- makeCluster(35)
registerDoParallel(cl)
###### PROBS
# foreach(i=1:length(clim_files),.packages=c('raster','rgdal'))%dopar%{
#
# ### READ
# climGrid <- read.csv(clim_files[i],header=TRUE,stringsAsFactors=FALSE)
# names(climGrid)[4:5] <- c("tp","pp")
#
# # set NA in order to solve
# climGrid$tp[which(climGrid$tp == -9999)] <- NA
# climGrid$pp[which(climGrid$pp == -9999)] <- NA
#
# #### SOLVE
# probsGrid <- solve_stm(climGrid_scale,pars)
# probsGrid <- data.frame(x=climGrid$x,y=climGrid$y,probsGrid,stringsAsFactors=FALSE)
#
# # get metadata
# filename <- strsplit(clim_files[i],"[/.]")[[1]][4]
#
# # write
# write.csv(probsGrid,file=paste0("./data/futStatesGrid/probs/",filename,".csv"),row.names=FALSE)
#
# }
###### STATES
probs_files <- list.files("./data/futStatesGrid/probs",full.names=TRUE,pattern="1985-2000")
#load raster ref
ref_rs100 <- readRDS("data/futClimGrid/rasterClim/res100/rcp85-ACCESS1_0-1985-2000.rda")
ref_rs1000 <- readRDS("data/futClimGrid/rasterClim/res1000/rcp85-ACCESS1_0-1985-2000.rda")
foreach(i=1:length(probs_files),.packages=c('raster','rgdal'))%dopar%{
probsGrid <- read.csv(probs_files[i])
states <- character(0)
for(r in 1:nrow(probsGrid)){
if(any(is.na(probsGrid[r,3:5]))){
states <- append(states,NA)
} else {
states <- append(states,names(which.max(probsGrid[r,3:5])))
}
}
stmGrid <- data.frame(coordinates(ref_rs1000),state=stateToId(states),stringsAsFactors=FALSE)
# turn into raster using clim rs
rs_stmGrid <- stmGrid
coordinates(rs_stmGrid) <- ~ x + y
gridded(rs_stmGrid) <- TRUE
rs_stmGrid <- raster(rs_stmGrid)
projection(rs_stmGrid) <- projection(ref_rs1000)
########## CREATE RES 1000
df_stmGrid_1000 <- as.data.frame(rs_stmGrid,xy=TRUE)
# turn coord into facteur
df_stmGrid_1000$x <- as.numeric(as.factor(df_stmGrid_1000$x)) - 1
df_stmGrid_1000$y <- as.numeric(as.factor(df_stmGrid_1000$y)) - 1
df_stmGrid_1000$state <- idToState(df_stmGrid_1000$state)
df_stmGrid_1000[is.na(df_stmGrid_1000$state),"state"] <- 0
########## CREATE RES 100
# reample using ref_rs100
rs_stmGrid_100 <- resample(rs_stmGrid,ref_rs100,method="ngb")
df_stmGrid_100 <- as.data.frame(rs_stmGrid_100,xy=TRUE)
# get metadata
filename <- strsplit(probs_files[i],"[/.]")[[1]][6]
filename <- paste(strsplit(filename,"[-]")[[1]][c(1,2,4)],collapse="-")
# turn coord into facteur
df_stmGrid_100$x <- as.numeric(as.factor(df_stmGrid_100$x)) - 1
df_stmGrid_100$y <- as.numeric(as.factor(df_stmGrid_100$y)) - 1
df_stmGrid_100$state <- idToState(df_stmGrid_100$state)
df_stmGrid_100[is.na(df_stmGrid_100$state),"state"] <- 0
# save for stm input
write.table(df_stmGrid_100,file=paste0("./data/futStatesGrid/stm100/",filename,".stm"),quote=FALSE,row.names=FALSE,col.names=FALSE,sep=",")
write.table(df_stmGrid_1000,file=paste0("./data/futStatesGrid/stm1000/",filename,".stm"),quote=FALSE,row.names=FALSE,col.names=FALSE,sep=",")
}