-
Notifications
You must be signed in to change notification settings - Fork 12
/
update.FM.model.R
218 lines (171 loc) · 8.74 KB
/
update.FM.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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
update.FM.model <- function(site_num) {
# The function update.FM.model updates an already existing particle filter
# forecast model. It first checks for new data, then assimilates that data into
# that forecast using a resampling particle filter. Outputs are generated one day
# at a time. Days with no new data are ignored, and the previous forecast values
# for that day are used instead. The function stops when all observed data has
# been assimilated.
# The forecast for each day is plotted and saved to a pdf begining with
# ParticleFilterForecast (with a site number and date appended). The output from
# the current forecast is saved in a file begining with ForecastModel.X.out (with
# a site number and date appended).
source("SSLPM.R")
source("ciEnvelope.R")
source("find.extreme.GCC.NDVI.R")
current.year <- strftime(Sys.Date(),"%Y")
source("global_input_parameters.R")
model.start.DOY <- global_input_parameters$model.start.DOY
##### get the date of the last forecast:
last.date.filename <- paste("last.update.site", as.character(site_num),
"txt",sep=".")
read.in <- source(last.date.filename)
last.forecast.date <- as.Date(read.in$value)
last.date.assimilated <- last.forecast.date
# load the GCC data:
gcc.data <- read.csv( sprintf("gcc_data_site%i.csv",site_num) )
# load the NDVI data:
ndvi.data <- read.csv( sprintf("ndvi_data_site%i.csv",site_num) )
# Merge them:
all.data <- merge(gcc.data,ndvi.data)
##### Rescale the data:
# find max/min of ndvi and gcc over all years of record except current
# outputs (ndvi_max,ndvi_min,gcc_max,gcc_min)
first.year <- as.numeric(strftime(global_input_parameters$data.start.date, "%Y"))
max_min_ndvi_gcc = find.extreme.GCC.NDVI(site_num, first.year,
as.numeric(current.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):
# rescale NDVI (and overwrite all.data$ndvi!)
all.data$ndvi <- (all.data$ndvi-ndvi_min)/(ndvi_max-ndvi_min)
# rescale GCC:
all.data$gcc.90 <- (all.data$gcc.90 - gcc_min)/(gcc_max - gcc_min)
all.data$gcc.mean <- (all.data$gcc.mean - gcc_min)/(gcc_max - gcc_min)
all.data$gcc.min <- (all.data$gcc.min - gcc_min)/(gcc_max - gcc_min)
all.data$gcc.max <- (all.data$gcc.max - gcc_min)/(gcc_max - gcc_min)
# load the forecast model output:
output_file_name = paste("ForecastModel.X.out.site", as.character(site_num),
"RData",sep=".")
load(output_file_name) # loads a num.days x num.ensemble x 2 array called output
# Number of ensemble members:
num.ensemble <- dim(output)[2]
forecast.date <- last.forecast.date + 1
current.date <- Sys.Date()
# Get standard deviations for measurement error from tau_gcc and tau_ndvi from
# our state-space model for now?
file_name = paste('Jags.SS.out.site',as.character(site_num), 'RData',sep=".")
load(file_name)
# get the median precisions from the state space model output, convert to stdevs:
tau.gcc.all <- jags.out.all.years.array[,3,] # num.ensemble members x num.years
gcc.stdev <- 1/sqrt(median(as.vector(tau.gcc.all)))
tau.ndvi.all <- jags.out.all.years.array[,4,] # num.ensemble members x num.years
ndvi.stdev <- 1/sqrt(median(as.vector(tau.ndvi.all)))
#### process error:
# Get process error from the SS model output (tau_add):
tau.add.all <- jags.out.all.years.array[,2,] # num.ensemble members x num.years
process.stdev <- 1/sqrt(median(as.vector(tau.add.all)))
# while loop until you get to the present day:
repeat{
# Keep this break statement floating at the top of the repeat loop:
if(forecast.date > current.date) {break} # This will end the loop
print(paste("Running particle filter for",forecast.date,"at site",site_num))
todays.data <- all.data[as.Date(all.data$date) == forecast.date,]
new.data <- !(is.na(todays.data$gcc.90) & is.na(todays.data$ndvi)) # TRUE/FALSE
# Only need to do anything when there is new data
if(new.data) {
# Let's get today's incoming X and r values:
output.days <- dim(output)[1]
output.index <- output.days - as.numeric(as.Date(paste(current.year,"12-31",sep="-")) - forecast.date,
unit="days")
X <- output[output.index,,1] # vector
r <- output[output.index,,2] # vector
#### Analysis step:
# Calculate the likelihood of our ensemble members given the data:
if(is.na(todays.data$ndvi)){
likelihood.ndvi <- rep(0,3000) # no likelihood if no data...
} else {
log.likelihood.ndvi <- dnorm(X,todays.data$ndvi,ndvi.stdev,log=TRUE)
}
if(is.na(todays.data$gcc.90)){
likelihood.gcc <- rep(0,3000) # no likelihood if no data...
} else {
log.likelihood.gcc <- dnorm(X,todays.data$gcc.90,gcc.stdev,log=TRUE)
}
likelihood <- exp(log.likelihood.gcc + log.likelihood.ndvi)
# if there is an outlier, so bad that it crashed the model, we set
# the likelihoods to all the same (smallish) value
if (sum(likelihood)==0){
likelihood = rep(0.00001,length(likelihood))
}
#### Resampling step:
index = sample.int(num.ensemble, num.ensemble, replace = TRUE, prob = likelihood)
# replace our previous guess with the PF output:
output[output.index,,1] = X[index] # or maybe pmin(1,pmax(0,X[index]
output[output.index,,2] = r[index]
#### Forecast step:
# as long as we're not at the end of the year:
if(forecast.date < as.Date(paste(current.year,"12-31",sep="-"))) {
# Forecast!
for(t in (output.index+1):output.days){
X = output[t-1,,1]
r = output[t-1,,2]
## forward step
output[t,,] = SSLPM(X,r) # num.ensembles x 2
# Add the process error to the state estimate:
output[t,,1] = output[t,,1] + rnorm(num.ensemble,0,process.stdev)
}
}
##### end of forecast loop
# Plot the forecast!
X.mat = output[,,1]
X.ci = apply(X.mat,1,quantile,c(0.025,0.5,0.975))
#### save plot produced to PDF
## name of output file
dir.name <- paste("pdfs/site",as.character(site_num),sep="")
## name of output file
pdf.file.name = paste("ParticleFilterForecast",as.character(site_num),
as.character(forecast.date),"pdf",sep=".")
## saves as PDF
pdf(file=paste(dir.name,pdf.file.name,sep="/"))
#### plot forecast:
# get rid of data from the future!
plottable.data <- subset(all.data,as.Date(all.data$date) <= forecast.date)
# get rid of data from previous years:
plottable.data <- subset(plottable.data,
strftime(as.Date(plottable.data$date),"%Y") == current.year)
# Get rid of early part of year:
plottable.data <- subset(plottable.data,
as.Date(plottable.data$date) >= model.start.DOY)
plot(model.start.DOY:365,X.ci[2,],type='n',
main=paste("Particle Filter Forecast:",forecast.date),
xlab="Day of Year",ylab="Pheno-state",ylim=c(0,1.2))
ciEnvelope(model.start.DOY:365,X.ci[1,],X.ci[3,],col="light grey")
lines(model.start.DOY:365,X.ci[2,],
main=paste("Particle Filter Forecast:",forecast.date),
xlab="Day of Year",ylab="Pheno-state")
non.leap.year.doys <- as.numeric(strftime(plottable.data$date,"%j")) - (as.numeric(current.year)%%4 == 0)
points(non.leap.year.doys, plottable.data$ndvi, pch="+",cex=0.8)
points(non.leap.year.doys, plottable.data$gcc.90, pch="o",cex=0.5)
## ends plot output to PDF
dev.off()
#### append output to pdf files that were created in the forecast model:
# This is important as it is the date to save in the file tracking the last
# date assimilated
last.date.assimilated <- forecast.date
} # end if(new.data)
# Increment the date, and update again!
forecast.date <- forecast.date + 1
}
# Save the most recent output data to file:
# Write the last forecast date to file:
date.string <- as.character(last.date.assimilated)
last.date.filename <- paste("last.update.site", as.character(site_num),
"txt",sep=".")
sink(last.date.filename, append = FALSE)
cat("\"",date.string,"\"",sep="")
sink()
}