forked from EcoForecast/PhenologyForecast
-
Notifications
You must be signed in to change notification settings - Fork 2
/
run.SS.model.R
160 lines (131 loc) · 6.62 KB
/
run.SS.model.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
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
run.SS.model <- function(site.number){
# source functions
source("find.extreme.GCC.NDVI.R")
source("RunJAGS.R")
source("make.SS.plots.R")
require(plyr)
# load site gcc and ndvi data:
gcc.filename <- sprintf("gcc_data_site%i.csv",site.number)
ndvi.filename <- sprintf("ndvi_data_site%i.csv",site.number)
gcc.data <- read.csv(gcc.filename)
ndvi.data <- read.csv(ndvi.filename)
# combine them, to make sure that dates match properly (they should already, but
# just in case...)
all.data <- merge(gcc.data,ndvi.data)
# Throw out years with no data!
no.data.vec <- is.na(all.data$gcc.mean) & is.na(all.data$ndvi) # TRUE/FALSE vector of na locations
years <- as.numeric(strftime(all.data$date, "%Y"))
years.with.data <- unique(years[!no.data.vec])
all.data <- subset(all.data, years %in% years.with.data )
current.year <- as.numeric(strftime(Sys.Date(), "%Y"))
if(!is.null(global_input_parameters$training.end.date)){
current.year = (as.numeric(strftime(global_input_parameters$training.end.date,"%Y"))+1):current.year
}
# Define fall dates:
source("global_input_parameters.R")
first.day.of.season <- global_input_parameters$model.start.DOY
day.of.year <- as.numeric( strftime(all.data$date, format="%j") ) # a vector of the DOY for each date
model = global_input_parameters$model
# We are only going to use the fall data:
fall.days <- (day.of.year >= first.day.of.season) & (day.of.year < 366) # solves the leap year problem
if(global_input_parameters$season == "FALL"){
all.data = all.data[fall.days,]
} else {
if(global_input_parameters$season == "SPRING"){
all.data = all.data[!fall.days,]
}
}
## build data object for JAGS in year loop.
# pull out time, get year for each data point.
time = all.data$date
time_year = as.numeric(format(as.Date(time), "%Y"))
# Define the number of iterations and chains
n.iter = global_input_parameters$number.of.SS.model.iterations;
n.chains = global_input_parameters$number.of.SS.model.chains;
# Stores output from each year
num.years.data <- length(unique(time_year)) # number of years (including the current
# year which won't go into the SS model)
# find max/min of ndvi and gcc over all years of record except current
# outputs (ndvi_max,ndvi_min,gcc_max,gcc_min)
max_min_ndvi_gcc = find.extreme.GCC.NDVI(site.number, min(time_year), max(time_year)-1,
use.interannual.means=TRUE)
ndvi_max = max_min_ndvi_gcc[1]
ndvi_min = max_min_ndvi_gcc[2]
gcc_max = max_min_ndvi_gcc[3]
gcc_min = max_min_ndvi_gcc[4]
# Rescale data to be between 0 and 1 (using max and min NDVI, GCC values from
# all years except current year):
rescaled_NDVI <- (all.data$ndvi-ndvi_min)/(ndvi_max-ndvi_min)
rescaled_NDVI[rescaled_NDVI < -0.3] = NA
rescaled_GCC <- (all.data$gcc.mean-gcc_min)/(gcc_max-gcc_min) # using gcc90
ratio_scale_NDVI = (max(rescaled_NDVI,na.rm=TRUE)-min(rescaled_NDVI,na.rm=TRUE))/(max(all.data$ndvi,na.rm=TRUE) - min(all.data$ndvi,na.rm=TRUE))
ratio_scale_GCC = (max(rescaled_GCC,na.rm=TRUE)-min(rescaled_GCC,na.rm=TRUE))/(max(all.data$gcc.mean,na.rm=TRUE) - min(all.data$gcc.mean,na.rm=TRUE))
SS.years <- setdiff(years.with.data,current.year)
n = max(table(time_year))
y = z = matrix(0.0,length(SS.years),n)
for (i in 1:length(SS.years)) { # loop over all years except current year
II = which(time_year == SS.years[i])
y[i,] = rescaled_NDVI[II] # get ndvi just for ONE year
z[i,] = rescaled_GCC[II] # get gcc just for ONE year
}
# Make list of "data" to be used as input for RunJAGS
x_ic <- global_input_parameters$x_ic
tau_ic <- global_input_parameters$tau_ic
a_ndvi <- (global_input_parameters$a_ndvi)*ratio_scale_NDVI
r_ndvi <- global_input_parameters$r_ndvi
a_gcc <- global_input_parameters$a_gcc
r_gcc <- global_input_parameters$r_gcc
a_add <- (global_input_parameters$a_add)*ratio_scale_GCC
r_add <- global_input_parameters$r_add
data <- list(y = y,z = z,n=n, x_ic=x_ic, tau_ic=tau_ic,
a_ndvi=a_ndvi, r_ndvi=r_ndvi, a_gcc=a_gcc, r_gcc=r_gcc,
a_add=a_add, r_add=r_add,ny = nrow(y))
# run JAGS model
jags.out=RunJAGS(data,n.iter,n.chains)
jags.out.matrix <- as.matrix(jags.out)
is.x = grep(pattern = "x[",x = colnames(jags.out.matrix),fixed = TRUE)
x.out = array(t(jags.out.matrix[,is.x]),c(data$ny,data$n,nrow(jags.out.matrix)))
ci = aaply(x.out,1:2,function(x) quantile(x,c(0.025,0.5,0.975)))
out = list(parms = jags.out.matrix[,-is.x],ci=ci)
if(FALSE){
# output parameters: [r tau_add tau_gcc tau_ndvi x], x will use (366-first.day.of.season columns)
jags.out.all.years.array = array(rep(NA, n.iter*n.chains*(4+366-first.day.of.season)*(num.years.data-1)),
c(n.iter*n.chains,(4+366-first.day.of.season),(num.years.data-1)))
# counts for loop
count = 0
for (YR in SS.years) { # loop over all years except current year
cat(sprintf("Running state-space model for site %i, year %i\n\n",site.number,YR))
count = count + 1;
# gets index of year
II = which(time_year == YR)
rescaled_NDVI_one_year = rescaled_NDVI[II] # get ndvi just for ONE year
rescaled_GCC_one_year = rescaled_GCC[II] # get gcc just for ONE year
# Make list of "data" to be used as input for RunJAGS
x_ic <- global_input_parameters$x_ic
tau_ic <- global_input_parameters$tau_ic
a_ndvi <- (global_input_parameters$a_ndvi)*ratio_scale_NDVI
r_ndvi <- global_input_parameters$r_ndvi
a_gcc <- global_input_parameters$a_gcc
r_gcc <- global_input_parameters$r_gcc
a_add <- (global_input_parameters$a_add)*ratio_scale_GCC
r_add <- global_input_parameters$r_add
data <- list(y = rescaled_NDVI_one_year,z = rescaled_GCC_one_year,
n=length(rescaled_NDVI_one_year),x_ic=x_ic, tau_ic=tau_ic,
a_ndvi=a_ndvi, r_ndvi=r_ndvi, a_gcc=a_gcc, r_gcc=r_gcc,
a_add=a_add, r_add=r_add)
# run JAGS model
jags.out=RunJAGS(data,n.iter,n.chains)
jags.out.matrix <- as.matrix(jags.out)
jags.out.all.years.array[,,count] <- jags.out.matrix
}
}
# Make filename and save jags output
file_name = paste('Jags.SS.out.site',as.character(site.number),model,'RData',sep=".")
# save(jags.out.all.years.array, SS.years, file = file_name)
save(out, SS.years, file = file_name)
# Make plots
# make.SS.plots(jags.out.all.years.array,all.data$date,
# rescaled_NDVI,rescaled_GCC,site.number)
make.SS.plots(out,all.data$date,
rescaled_NDVI,rescaled_GCC,site.number,model)
}