Revision 3e1b1ed4
Added by Benoit Parmentier almost 12 years ago
climate/research/oregon/interpolation/GAM_fusion_function_multisampling.R | ||
---|---|---|
5 | 5 |
#TODO: |
6 | 6 |
#Add log file and calculate time and sizes for processes-outputs |
7 | 7 |
|
8 |
|
|
9 |
runClimFusion<-function(j){ |
|
8 |
runClim_KGFusion<-function(j){ |
|
10 | 9 |
#Make this a function with multiple argument that can be used by mcmapply?? |
11 | 10 |
#This creates clim fusion layers... |
12 | 11 |
|
... | ... | |
75 | 74 |
data_month$LSTD_bias<-data_month$LST-data_month$TMax |
76 | 75 |
data_month$y_var<-data_month$LSTD_bias #Adding bias as the variable modeled |
77 | 76 |
mod_list<-fit_models(list_formulas,data_month) #only gam at this stage |
78 |
cname<-paste("mod",1:length(mod_list),sep="") |
|
77 |
cname<-paste("mod",1:length(mod_list),sep="") #change to more meaningful name?
|
|
79 | 78 |
names(mod_list)<-cname |
80 | 79 |
#Adding layer LST to the raster stack |
81 |
pos<-match("elev",layerNames(s_raster))
|
|
80 |
pos<-match("elev",names(s_raster))
|
|
82 | 81 |
layerNames(s_raster)[pos]<-"elev_1" |
83 | 82 |
|
84 | 83 |
pos<-match("LST",names(s_raster)) #Find the position of the layer with name "LST", if not present pos=NA |
... | ... | |
91 | 90 |
s_raster<-addLayer(s_raster,LST) #Adding current month |
92 | 91 |
|
93 | 92 |
#Now generate file names for the predictions... |
94 |
list_out_filename<-vector("list",length(in_models)) |
|
95 |
|
|
93 |
list_out_filename<-vector("list",length(mod_list)) |
|
94 |
names(list_out_filename)<-cname |
|
95 |
|
|
96 | 96 |
for (k in 1:length(list_out_filename)){ |
97 |
data_name<-paste("bias_LST_month_",j,"_mod_",k,"_",prop_month, |
|
97 |
#j indicate which month is predicted |
|
98 |
data_name<-paste("bias_LST_month_",j,"_",cname[k],"_",prop_month, |
|
98 | 99 |
"_",run_samp,sep="") |
99 | 100 |
raster_name<-paste("fusion_",data_name,out_prefix,".tif", sep="") |
100 | 101 |
list_out_filename[[k]]<-raster_name |
101 | 102 |
} |
102 | 103 |
|
103 | 104 |
#now predict values for raster image... |
104 |
rast_list<-predict_raster_model(mod_list,s_raster,out_filename) |
|
105 |
rast_list<-rast_list[!sapply(rast_list,is.null)] #remove NULL elements in list |
|
106 |
|
|
107 |
mod_rast<-stack(rast_list) |
|
105 |
rast_bias_list<-predict_raster_model(mod_list,s_raster,list_out_filename) |
|
106 |
names(rast_bias_list)<-cname |
|
107 |
#Some modles will not be predicted...remove them |
|
108 |
rast_bias_list<-rast_bias_list[!sapply(rast_bias_list,is.null)] #remove NULL elements in list |
|
109 |
|
|
110 |
mod_rast<-stack(rast_bias_list) #stack of bias raster images from models |
|
108 | 111 |
rast_clim_list<-vector("list",nlayers(mod_rast)) |
112 |
names(rast_clim_list)<-names(rast_bias_list) |
|
109 | 113 |
for (k in 1:nlayers(mod_rast)){ |
110 | 114 |
clim_fus_rast<-LST-subset(mod_rast,k) |
111 |
data_name<-paste("clim_LST_month_",j,"_mod_",k,"_",prop_month,
|
|
115 |
data_name<-paste("clim_LST_month_",j,"_",names(rast_clim_list)[k],"_",prop_month,
|
|
112 | 116 |
"_",run_samp,sep="") |
113 | 117 |
raster_name<-paste("fusion_",data_name,out_prefix,".tif", sep="") |
114 | 118 |
rast_clim_list[[k]]<-raster_name |
115 | 119 |
writeRaster(clim_fus_rast, filename=raster_name,overwrite=TRUE) #Wri |
116 | 120 |
} |
117 |
clim_obj<-list(rast_list,rast_clim_list,data_month,mod_list,list_formulas) |
|
121 |
|
|
122 |
#Adding Kriging for Climatology options |
|
123 |
|
|
124 |
bias_xy<-coordinates(data_month) |
|
125 |
fitbias<-Krig(bias_xy,data_month$LSTD_bias,theta=1e5) #use TPS or krige |
|
126 |
mod_krtmp1<-fitbias |
|
127 |
model_name<-"mod_kr" |
|
128 |
|
|
129 |
bias_rast<-interpolate(LST,fitbias) #interpolation using function from raster package |
|
130 |
#Saving kriged surface in raster images |
|
131 |
data_name<-paste("bias_LST_month_",j,"_",model_name,"_",prop_month, |
|
132 |
"_",run_samp,sep="") |
|
133 |
raster_name_bias<-paste("fusion_",data_name,out_prefix,".tif", sep="") |
|
134 |
writeRaster(bias_rast, filename=raster_name_bias,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
135 |
|
|
136 |
#now climatology layer |
|
137 |
clim_rast<-LST-bias_rast |
|
138 |
data_name<-paste("clim_LST_month_",j,"_",model_name,"_",prop_month, |
|
139 |
"_",run_samp,sep="") |
|
140 |
raster_name_clim<-paste("fusion_",data_name,out_prefix,".tif", sep="") |
|
141 |
writeRaster(clim_rast, filename=raster_name_clim,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
142 |
|
|
143 |
#Adding to current objects |
|
144 |
mod_list[[model_name]]<-mod_krtmp1 |
|
145 |
rast_bias_list[[model_name]]<-raster_name_bias |
|
146 |
rast_clim_list[[model_name]]<-raster_name_clim |
|
147 |
|
|
148 |
#Prepare object to return |
|
149 |
clim_obj<-list(rast_bias_list,rast_clim_list,data_month,mod_list,list_formulas) |
|
118 | 150 |
names(clim_obj)<-c("bias","clim","data_month","mod","formulas") |
151 |
|
|
152 |
save(clim_obj,file= paste("clim_obj_month_",j,"_",out_prefix,".RData",sep="")) |
|
153 |
|
|
119 | 154 |
return(clim_obj) |
120 | 155 |
} |
121 | 156 |
|
122 | 157 |
## Run function for kriging...? |
123 | 158 |
|
124 |
|
|
125 | 159 |
runGAMFusion <- function(i) { # loop over dates |
160 |
#Change this to allow explicitly arguments... |
|
161 |
#Arguments: |
|
162 |
#1)list of climatology files for all models...(12*nb of models) |
|
163 |
#2)data_s:training |
|
164 |
#3)data_v:testing |
|
165 |
#4)list of dates?? |
|
166 |
#5)stack of covariates: not needed at this this stage |
|
167 |
#6)dst: data at the monthly time scale |
|
126 | 168 |
|
127 | 169 |
#Function used in the script |
128 | 170 |
|
129 |
lookup<-function(r,lat,lon) { |
|
130 |
#This functions extracts values in a projected raster |
|
131 |
#given latitude and longitude vector locations |
|
132 |
#Output: matrix with values |
|
133 |
xy<-project(cbind(lon,lat),proj_str); |
|
134 |
cidx<-cellFromXY(r,xy); |
|
135 |
return(r[cidx]) |
|
136 |
} |
|
137 |
|
|
138 |
|
|
139 | 171 |
date<-strptime(sampling_dat$date[i], "%Y%m%d") # interpolation date being processed |
140 | 172 |
month<-strftime(date, "%m") # current month of the date being processed |
141 | 173 |
LST_month<-paste("mm_",month,sep="") # name of LST month to be matched |
142 | 174 |
proj_str<-proj4string(dst) |
143 |
#Adding layer LST to the raster stack |
|
144 | 175 |
|
145 |
pos<-match("LST",layerNames(s_raster)) #Find the position of the layer with name "LST", if not present pos=NA |
|
146 |
s_raster<-dropLayer(s_raster,pos) # If it exists drop layer |
|
147 |
pos<-match(LST_month,layerNames(s_raster)) #Find column with the current month for instance mm12 |
|
148 |
r1<-raster(s_raster,layer=pos) #Select layer from stack |
|
149 |
layerNames(r1)<-"LST" |
|
150 |
#Screen for extreme values" 10/30 |
|
151 |
#min_val<-(-15+273.16) #if values less than -15C then screen out (note the Kelvin units that will need to be changed later in all datasets) |
|
152 |
#r1[r1 < (min_val)]<-NA |
|
153 |
s_raster<-addLayer(s_raster,r1) #Adding current month |
|
154 |
|
|
155 |
pos<-match("elev",layerNames(s_raster)) |
|
156 |
layerNames(s_raster)[pos]<-"elev_1" |
|
157 | 176 |
###Regression part 1: Creating a validation dataset by creating training and testing datasets |
158 | 177 |
data_day<-ghcn.subsets[[i]] |
159 | 178 |
mod_LST <- ghcn.subsets[[i]][,match(LST_month, names(ghcn.subsets[[i]]))] #Match interpolation date and monthly LST average |
... | ... | |
173 | 192 |
mo<-as.integer(strftime(date_proc, "%m")) # current month of the date being processed |
174 | 193 |
day<-as.integer(strftime(date_proc, "%d")) |
175 | 194 |
year<-as.integer(strftime(date_proc, "%Y")) |
176 |
|
|
177 |
datelabel=format(ISOdate(year,mo,day),"%b %d, %Y") |
|
178 |
|
|
179 |
########### |
|
180 |
# STEP 1 - LST 10 year monthly averages |
|
181 |
########### |
|
182 |
pos<-match("LST",layerNames(s_raster)) #Find the position of the layer with name "LST", |
|
183 |
themolst<-raster(s_raster,layer=pos) |
|
184 |
#themolst<-raster(molst,mo) #current month being processed saved in a raster image |
|
185 |
#min_val<-(-15) #Screening for extreme values |
|
186 |
#themolst[themolst < (min_val)]<-NA |
|
187 |
|
|
188 |
########### |
|
189 |
# STEP 2 - Weather station means across same days: Monthly mean calculation |
|
190 |
########### |
|
191 | 195 |
|
192 | 196 |
modst<-dst[dst$month==mo,] #Subsetting dataset for the relevant month of the date being processed |
193 |
|
|
194 |
########## |
|
195 |
# STEP 3 - get LST at stations |
|
196 |
########## |
|
197 |
|
|
198 |
#sta_lola=modst[,c("lon","lat")] #Extracting locations of stations for the current month.. |
|
199 |
|
|
200 |
#proj_str="+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs"; |
|
197 |
#Change to y_var...could be TMin |
|
198 |
#modst$LSTD_bias <- modst$LST-modst$y_var |
|
199 |
modst$LSTD_bias <- modst$LST-modst$TMax; #That is the difference between the monthly LST mean and monthly station mean |
|
201 | 200 |
|
202 |
#sta_tmax_from_lst=lookup(themolst,sta_lola$lat,sta_lola$lon) #Extracted values of LST for the stations |
|
203 |
sta_tmax_from_lst<-modst$LST |
|
204 |
######### |
|
205 |
# STEP 4 - bias at stations |
|
206 |
######### |
|
207 |
|
|
208 |
sta_bias=sta_tmax_from_lst-modst$TMax; #That is the difference between the monthly LST mean and monthly station mean |
|
209 |
modst$LSTD_bias<-sta_bias #Adding bias to data frame modst containning the monthly average for 10 years |
|
210 |
|
|
211 |
#bias_xy=project(as.matrix(sta_lola),proj_str) |
|
212 |
png(paste("LST_TMax_scatterplot_",sampling_dat$date[i],"_",sampling_dat$prop[i],"_",sampling_dat$run_samp[i], out_prefix,".png", sep="")) |
|
213 |
plot(modst$TMax,sta_tmax_from_lst,xlab="Station mo Tmax",ylab="LST mo Tmax",main=paste("LST vs TMax for",datelabel,sep=" ")) |
|
214 |
abline(0,1) |
|
215 |
nb_point<-paste("n=",length(modst$TMax),sep="") |
|
216 |
mean_bias<-paste("LST bias= ",format(mean(modst$LSTD_bias,na.rm=TRUE),digits=3),sep="") |
|
217 |
#Add the number of data points on the plot |
|
218 |
legend("topleft",legend=c(mean_bias,nb_point),bty="n") |
|
219 |
dev.off() |
|
220 |
|
|
221 |
#added by Benoit |
|
222 |
#x<-ghcn.subsets[[i]] #Holds both training and testing for instance 161 rows for Jan 1 |
|
223 | 201 |
x<-as.data.frame(data_v) |
224 | 202 |
d<-as.data.frame(data_s) |
225 | 203 |
#x[x$value==-999.9]<-NA |
... | ... | |
245 | 223 |
names(x)[pos]<-c("id") |
246 | 224 |
names(modst)[1]<-c("id") #modst contains the average tmax per month for every stations... |
247 | 225 |
|
248 |
dmoday=merge(modst,d,by="id",suffixes=c("",".y2")) #LOOSING DATA HERE!!! from 113 t0 103
|
|
249 |
xmoday=merge(modst,x,by="id",suffixes=c("",".y2")) #LOOSING DATA HERE!!! from 48 t0 43
|
|
226 |
dmoday <-merge(modst,d,by="id",suffixes=c("",".y2")) #LOOSING DATA HERE!!! from 113 t0 103
|
|
227 |
xmoday <-merge(modst,x,by="id",suffixes=c("",".y2")) #LOOSING DATA HERE!!! from 48 t0 43
|
|
250 | 228 |
mod_pat<-glob2rx("*.y2") |
251 | 229 |
var_pat<-grep(mod_pat,names(dmoday),value=FALSE) # using grep with "value" extracts the matching names |
252 | 230 |
dmoday<-dmoday[,-var_pat] |
... | ... | |
255 | 233 |
xmoday<-xmoday[,-var_pat] #Removing duplicate columns |
256 | 234 |
|
257 | 235 |
data_v<-xmoday |
258 |
### |
|
259 | 236 |
|
260 | 237 |
#dmoday contains the daily tmax values for training with TMax being the monthly station tmax mean |
261 | 238 |
#xmoday contains the daily tmax values for validation with TMax being the monthly station tmax mean |
262 | 239 |
|
263 |
png(paste("Daily_tmax_monthly_TMax_scatterplot_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
264 |
"_",sampling_dat$run_samp[i],out_prefix,".png", sep="")) |
|
265 |
plot(dailyTmax~TMax,data=dmoday,xlab="Mo Tmax",ylab=paste("Daily for",datelabel),main="across stations in OR") |
|
266 |
dev.off() |
|
267 |
|
|
268 |
######## |
|
269 |
# STEP 5 - interpolate bias |
|
270 |
######## |
|
271 |
|
|
272 |
#Adding options to use only training stations : 07/11/2012 |
|
273 |
#bias_xy=project(as.matrix(sta_lola),proj_str) |
|
274 |
#bias_xy2=project(as.matrix(c(dmoday$lon,dmoday$lat),proj_str) |
|
275 |
bias_xy<-coordinates(modst) |
|
276 |
if(bias_val==1){ |
|
277 |
sta_bias<-dmoday$LSTD_bias |
|
278 |
bias_xy<-cbind(dmoday$x_OR83M,dmoday$y_OR83M) |
|
279 |
} |
|
280 |
|
|
281 |
fitbias<-Krig(bias_xy,sta_bias,theta=1e5) #use TPS or krige |
|
282 |
#The output is a krig object using fields: modif 10/30 |
|
283 |
#mod9a<-fitbias |
|
284 |
mod_krtmp1<-fitbias |
|
285 |
model_name<-paste("mod_kr","month",sep="_") |
|
286 |
assign(model_name,mod_krtmp1) |
|
287 |
|
|
288 |
|
|
289 | 240 |
########## |
290 | 241 |
# STEP 7 - interpolate delta across space |
291 | 242 |
########## |
292 | 243 |
|
293 |
daily_sta_lola=dmoday[,c("lon","lat")] #could be same as before but why assume merge does this - assume not |
|
294 |
daily_sta_xy=project(as.matrix(daily_sta_lola),proj_str) |
|
295 |
|
|
296 |
daily_delta=dmoday$dailyTmax-dmoday$TMax |
|
297 |
|
|
298 |
#fitdelta<-Tps(daily_sta_xy,daily_delta) #use TPS or krige |
|
299 |
fitdelta<-Krig(daily_sta_xy,daily_delta,theta=1e5) #use TPS or krige |
|
300 |
#Kriging using fields package: modif 10/30 |
|
301 |
#mod9b<-fitdelta |
|
244 |
daily_delta<-dmoday$dailyTmax-dmoday$TMax |
|
245 |
daily_delta_xy<-as.matrix(cbind(dmoday$x,dmoday$y)) |
|
246 |
fitdelta<-Krig(daily_delta_xy,daily_delta,theta=1e5) #use TPS or krige |
|
302 | 247 |
mod_krtmp2<-fitdelta |
303 | 248 |
model_name<-paste("mod_kr","day",sep="_") |
304 |
assign(model_name,mod_krtmp2) |
|
305 |
|
|
306 |
png(paste("Delta_surface_LST_TMax_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
307 |
"_",sampling_dat$run_samp[i],out_prefix,".png", sep="")) |
|
308 |
surface(fitdelta,col=rev(terrain.colors(100)),asp=1,main=paste("Interpolated delta for",datelabel,sep=" ")) |
|
309 |
dev.off() |
|
310 |
# |
|
311 |
|
|
312 |
#### Added by Benoit on 06/19 |
|
313 | 249 |
data_s<-dmoday #put the |
314 | 250 |
data_s$daily_delta<-daily_delta |
315 | 251 |
|
316 |
#data_s$y_var<-daily_delta #y_var is the variable currently being modeled, may be better with BIAS!! |
|
317 |
#data_s$y_var<-data_s$LSTD_bias |
|
318 |
#### Added by Benoit ends |
|
319 |
|
|
320 | 252 |
######### |
321 | 253 |
# STEP 8 - assemble final answer - T=LST+Bias(interpolated)+delta(interpolated) |
322 | 254 |
######### |
323 |
|
|
324 |
bias_rast=interpolate(themolst,fitbias) #interpolation using function from raster package |
|
325 |
#themolst is raster layer, fitbias is "Krig" object from bias surface |
|
326 |
#plot(bias_rast,main="Raster bias") #This not displaying... |
|
327 | 255 |
|
328 |
#Saving kriged surface in raster images |
|
329 |
data_name<-paste("bias_LST_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
330 |
"_",sampling_dat$run_samp[i],sep="") |
|
331 |
raster_name<-paste("fusion_",data_name,out_prefix,".rst", sep="") |
|
332 |
writeRaster(bias_rast, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
256 |
rast_clim_list<-rast_clim_yearlist[[mo]] #select relevant month |
|
257 |
rast_clim_month<-raster(rast_clim_list[[1]]) |
|
333 | 258 |
|
334 |
daily_delta_rast=interpolate(themolst,fitdelta) #Interpolation of the bias surface... |
|
335 |
|
|
336 |
#plot(daily_delta_rast,main="Raster Daily Delta") |
|
259 |
daily_delta_rast<-interpolate(rast_clim_month,fitdelta) #Interpolation of the bias surface... |
|
337 | 260 |
|
338 | 261 |
#Saving kriged surface in raster images |
339 | 262 |
data_name<-paste("daily_delta_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
340 | 263 |
"_",sampling_dat$run_samp[i],sep="") |
341 |
raster_name<-paste("fusion_",data_name,out_prefix,".rst", sep="") |
|
342 |
writeRaster(daily_delta_rast, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
343 |
|
|
344 |
tmax_predicted=themolst+daily_delta_rast-bias_rast #Final surface as a raster layer... |
|
345 |
#tmax_predicted=themolst+daily_delta_rast+bias_rast #Added by Benoit, why is it -bias_rast |
|
346 |
#plot(tmax_predicted,main="Predicted daily") |
|
347 |
|
|
348 |
#Saving kriged surface in raster images |
|
349 |
data_name<-paste("tmax_predicted_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
350 |
"_",sampling_dat$run_samp[i],sep="") |
|
351 |
raster_name<-paste("fusion_",data_name,out_prefix,".rst", sep="") |
|
352 |
writeRaster(tmax_predicted, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
353 |
|
|
354 |
######## |
|
355 |
# check: assessment of results: validation |
|
356 |
######## |
|
357 |
RMSE<-function(res) {return(((mean(res,na.rm=TRUE))^2)^0.5)} |
|
358 |
MAE_fun<-function(res) {return(mean(abs(res),na.rm=TRUE))} |
|
359 |
#ME_fun<-function(x,y){return(mean(abs(y)))} |
|
360 |
#FIT ASSESSMENT |
|
361 |
sta_pred_data_s=lookup(tmax_predicted,data_s$lat,data_s$lon) |
|
362 |
|
|
363 |
rmse_fit=RMSE(sta_pred_data_s-data_s$dailyTmax) |
|
364 |
mae_fit=MAE_fun(sta_pred_data_s-data_s$dailyTmax) |
|
264 |
raster_name_delta<-paste("fusion_",data_name,out_prefix,".tif", sep="") |
|
265 |
writeRaster(daily_delta_rast, filename=raster_name_delta,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
266 |
|
|
267 |
#Now predict daily after having selected the relevant month |
|
268 |
temp_list<-vector("list",length(rast_clim_list)) |
|
269 |
for (j in 1:length(rast_clim_list)){ |
|
270 |
rast_clim_month<-raster(rast_clim_list[[j]]) |
|
271 |
temp_predicted<-rast_clim_month+daily_delta_rast |
|
365 | 272 |
|
366 |
sta_pred=lookup(tmax_predicted,data_v$lat,data_v$lon) |
|
367 |
#sta_pred=lookup(tmax_predicted,daily_sta_lola$lat,daily_sta_lola$lon) |
|
368 |
#rmse=RMSE(sta_pred,dmoday$dailyTmax) |
|
369 |
#pos<-match("value",names(data_v)) #Find column with name "value" |
|
370 |
#names(data_v)[pos]<-c("dailyTmax") |
|
371 |
tmax<-data_v$dailyTmax |
|
372 |
#data_v$dailyTmax<-tmax |
|
373 |
rmse=RMSE(sta_pred-tmax) |
|
374 |
mae<-MAE_fun(sta_pred-tmax) |
|
375 |
r2<-cor(sta_pred,tmax)^2 #R2, coef. of var |
|
376 |
me<-mean(sta_pred-tmax,na.rm=T) |
|
377 |
|
|
378 |
png(paste("Predicted_tmax_versus_observed_scatterplot_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
379 |
"_",sampling_dat$run_samp[i],out_prefix,".png", sep="")) |
|
380 |
plot(sta_pred~tmax,xlab=paste("Actual daily for",datelabel),ylab="Pred daily",main=paste("RMSE=",rmse)) |
|
381 |
abline(0,1) |
|
382 |
dev.off() |
|
383 |
#resid=sta_pred-dmoday$dailyTmax |
|
384 |
resid=sta_pred-tmax |
|
385 |
|
|
386 |
###BEFORE GAM prediction the data object must be transformed to SDF |
|
387 |
|
|
388 |
coords<- data_v[,c('x','y')] |
|
389 |
coordinates(data_v)<-coords |
|
390 |
proj4string(data_v)<-proj_str #Need to assign coordinates... |
|
391 |
coords<- data_s[,c('x','y')] |
|
392 |
coordinates(data_s)<-coords |
|
393 |
proj4string(data_s)<-proj_str #Need to assign coordinates.. |
|
394 |
coords<- modst[,c('x','y')] |
|
395 |
#coordinates(modst)<-coords |
|
396 |
#proj4string(modst)<-proj_str #Need to assign coordinates.. |
|
397 |
|
|
398 |
ns<-nrow(data_s) #This is added to because some loss of data might have happened because of the averaging... |
|
399 |
nv<-nrow(data_v) |
|
400 |
|
|
401 |
###GAM PREDICTION |
|
402 |
|
|
403 |
if (bias_prediction==1){ |
|
404 |
data_s$y_var<-data_s$LSTD_bias #This shoudl be changed for any variable!!! |
|
405 |
data_v$y_var<-data_v$LSTD_bias |
|
406 |
data_month<-modst |
|
407 |
data_month$y_var<-modst$LSTD_bias |
|
408 |
} |
|
409 |
|
|
410 |
if (bias_prediction==0){ |
|
411 |
data_v$y_var<-data_v[[y_var_name]] |
|
412 |
data_s$y_var<-data_s[[y_var_name]] |
|
413 |
} |
|
414 |
|
|
415 |
#Model and response variable can be changed without affecting the script |
|
416 |
|
|
417 |
list_formulas<-vector("list",nmodels) |
|
418 |
|
|
419 |
list_formulas[[1]] <- as.formula("y_var ~ s(elev_1)", env=.GlobalEnv) |
|
420 |
list_formulas[[2]] <- as.formula("y_var ~ s(LST)", env=.GlobalEnv) |
|
421 |
list_formulas[[3]] <- as.formula("y_var ~ s(elev_1,LST)", env=.GlobalEnv) |
|
422 |
list_formulas[[4]] <- as.formula("y_var ~ s(lat) + s(lon)+ s(elev_1)", env=.GlobalEnv) |
|
423 |
list_formulas[[5]] <- as.formula("y_var ~ s(lat,lon,elev_1)", env=.GlobalEnv) |
|
424 |
list_formulas[[6]] <- as.formula("y_var ~ s(lat,lon) + s(elev_1) + s(N_w,E_w) + s(LST)", env=.GlobalEnv) |
|
425 |
list_formulas[[7]] <- as.formula("y_var ~ s(lat,lon) + s(elev_1) + s(N_w,E_w) + s(LST) + s(LC2)", env=.GlobalEnv) |
|
426 |
list_formulas[[8]] <- as.formula("y_var ~ s(lat,lon) + s(elev_1) + s(N_w,E_w) + s(LST) + s(LC6)", env=.GlobalEnv) |
|
427 |
list_formulas[[9]] <- as.formula("y_var ~ s(lat,lon) + s(elev_1) + s(N_w,E_w) + s(LST) + s(DISTOC)", env=.GlobalEnv) |
|
428 |
|
|
429 |
mod_list<-fit_models(list_formulas,data_month) #only gam at this stage |
|
430 |
|
|
431 |
#Added |
|
432 |
#tmax_predicted=themolst+daily_delta_rast-bias_rast #Final surface?? but daily_rst |
|
433 |
|
|
434 |
### Added by benoit |
|
435 |
#Store results using TPS |
|
436 |
j=nmodels+1 |
|
437 |
results_RMSE[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
438 |
results_RMSE[2]<- ns #number of stations used in the training stage |
|
439 |
results_RMSE[3]<- "RMSE" |
|
440 |
|
|
441 |
results_RMSE[j+3]<- rmse #Storing RMSE for the model j |
|
442 |
|
|
443 |
results_RMSE_f[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
444 |
results_RMSE_f[2]<- ns #number of stations used in the training stage |
|
445 |
results_RMSE_f[3]<- "RMSE_f" |
|
446 |
results_RMSE_f[j+3]<- rmse_fit #Storing RMSE for the model j |
|
447 |
|
|
448 |
results_MAE_f[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
449 |
results_MAE_f[2]<- ns #number of stations used in the training stage |
|
450 |
results_MAE_f[3]<- "RMSE_f" |
|
451 |
results_MAE_f[j+3]<- mae_fit #Storing RMSE for the model j |
|
452 |
|
|
453 |
results_MAE[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
454 |
results_MAE[2]<- ns #number of stations used in the training stage |
|
455 |
results_MAE[3]<- "MAE" |
|
456 |
results_MAE[j+3]<- mae #Storing RMSE for the model j |
|
457 |
|
|
458 |
results_ME[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
459 |
results_ME[2]<- ns #number of stations used in the training stage |
|
460 |
results_ME[3]<- "ME" |
|
461 |
results_ME[j+3]<- me #Storing RMSE for the model j |
|
462 |
|
|
463 |
results_R2[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
464 |
results_R2[2]<- ns #number of stations used in the training stage |
|
465 |
results_R2[3]<- "R2" |
|
466 |
results_R2[j+3]<- r2 #Storing RMSE for the model j |
|
467 |
|
|
468 |
#ns<-nrow(data_s) #This is added to because some loss of data might have happened because of the averaging... |
|
469 |
#nv<-nrow(data_v) |
|
470 |
|
|
471 |
pred_mod<-paste("pred_mod",j,sep="") |
|
472 |
#Adding the results back into the original dataframes. |
|
473 |
data_s[[pred_mod]]<-sta_pred_data_s |
|
474 |
data_v[[pred_mod]]<-sta_pred |
|
475 |
|
|
476 |
#Model assessment: RMSE and then krig the residuals....! |
|
477 |
|
|
478 |
res_mod_s<- data_s$dailyTmax - data_s[[pred_mod]] #Residuals from kriging training |
|
479 |
res_mod_v<- data_v$dailyTmax - data_v[[pred_mod]] #Residuals from kriging validation |
|
480 |
|
|
481 |
name2<-paste("res_mod",j,sep="") |
|
482 |
data_v[[name2]]<-as.numeric(res_mod_v) |
|
483 |
data_s[[name2]]<-as.numeric(res_mod_s) |
|
484 |
|
|
485 |
mod_obj<-vector("list",nmodels+2) #This will contain the model objects fitting: 10/30 |
|
486 |
mod_obj[[nmodels+1]]<-mod_kr_month #Storing climatology object |
|
487 |
mod_obj[[nmodels+2]]<-mod_kr_day #Storing delta object |
|
488 |
pred_sgdf<-as(raster_pred,"SpatialGridDataFrame") #Conversion to spatial grid data frame |
|
489 |
|
|
490 |
for (j in 1:nmodels){ |
|
491 |
|
|
492 |
##Model assessment: specific diagnostic/metrics for GAM |
|
493 |
|
|
494 |
name<-paste("mod",j,sep="") #modj is the name of The "j" model (mod1 if j=1) |
|
495 |
#mod<-get(name) #accessing GAM model ojbect "j" |
|
496 |
mod<-mod_list[[j]] |
|
497 |
mod_obj[[j]]<-mod #storing current model object |
|
498 |
#If mod "j" is not a model object |
|
499 |
if (inherits(mod,"try-error")) { |
|
500 |
results_AIC[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
501 |
results_AIC[2]<- ns #number of stations used in the training stage |
|
502 |
results_AIC[3]<- "AIC" |
|
503 |
results_AIC[j+3]<- NA |
|
504 |
|
|
505 |
results_GCV[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
506 |
results_GCV[2]<- ns #number of stations used in the training |
|
507 |
results_GCV[3]<- "GCV" |
|
508 |
results_GCV[j+3]<- NA |
|
509 |
|
|
510 |
results_DEV[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
511 |
results_DEV[2]<- ns #number of stations used in the training stage |
|
512 |
results_DEV[3]<- "DEV" |
|
513 |
results_DEV[j+3]<- NA |
|
514 |
|
|
515 |
results_RMSE_f[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
516 |
results_RMSE_f[2]<- ns #number of stations used in the training stage |
|
517 |
results_RMSE_f[3]<- "RSME_f" |
|
518 |
results_RMSE_f[j+3]<- NA |
|
519 |
|
|
520 |
results_MAE_f[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
521 |
results_MAE_f[2]<- ns #number of stations used in the training stage |
|
522 |
results_MAE_f[3]<- "MAE_f" |
|
523 |
results_MAE_f[j+3]<-NA |
|
524 |
|
|
525 |
results_RMSE[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
526 |
results_RMSE[2]<- ns #number of stations used in the training stage |
|
527 |
results_RMSE[3]<- "RMSE" |
|
528 |
results_RMSE[j+3]<- NA #Storing RMSE for the model j |
|
529 |
results_MAE[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
530 |
results_MAE[2]<- ns #number of stations used in the training stage |
|
531 |
results_MAE[3]<- "MAE" |
|
532 |
results_MAE[j+3]<- NA #Storing MAE for the model j |
|
533 |
results_ME[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
534 |
results_ME[2]<- ns #number of stations used in the training stage |
|
535 |
results_ME[3]<- "ME" |
|
536 |
results_ME[j+3]<- NA #Storing ME for the model j |
|
537 |
results_R2[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
538 |
results_R2[2]<- ns #number of stations used in the training stage |
|
539 |
results_R2[3]<- "R2" |
|
540 |
|
|
541 |
|
|
542 |
results_R2[j+3]<- NA #Storing R2 for the model j |
|
543 |
|
|
544 |
} |
|
545 |
|
|
546 |
#If mod is a modelobject |
|
547 |
|
|
548 |
#If mod "j" is not a model object |
|
549 |
if (inherits(mod,"gam")) { |
|
550 |
|
|
551 |
results_AIC[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
552 |
results_AIC[2]<- ns #number of stations used in the training stage |
|
553 |
results_AIC[3]<- "AIC" |
|
554 |
results_AIC[j+3]<- AIC (mod) |
|
555 |
|
|
556 |
results_GCV[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
557 |
results_GCV[2]<- ns #number of stations used in the training |
|
558 |
results_GCV[3]<- "GCV" |
|
559 |
results_GCV[j+3]<- mod$gcv.ubre |
|
560 |
|
|
561 |
results_DEV[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
562 |
results_DEV[2]<- ns #number of stations used in the training stage |
|
563 |
results_DEV[3]<- "DEV" |
|
564 |
results_DEV[j+3]<- mod$deviance |
|
565 |
|
|
566 |
y_var_fit= mod$fit |
|
567 |
|
|
568 |
results_RMSE_f[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
569 |
results_RMSE_f[2]<- ns #number of stations used in the training stage |
|
570 |
results_RMSE_f[3]<- "RSME_f" |
|
571 |
#results_RMSE_f[j+3]<- sqrt(sum((y_var_fit-data_s$y_var)^2)/ns) |
|
572 |
results_RMSE_f[j+3]<-sqrt(mean(mod$residuals^2,na.rm=TRUE)) |
|
573 |
|
|
574 |
results_MAE_f[1]<- sampling_dat$date[i] #storing the interpolation sampling_dat$date in the first column |
|
575 |
results_MAE_f[2]<- ns #number of stations used in the training stage |
|
576 |
results_MAE_f[3]<- "MAE_f" |
|
577 |
#results_MAE_f[j+3]<-sum(abs(y_var_fit-data_s$y_var))/ns |
|
578 |
results_MAE_f[j+3]<-mean(abs(mod$residuals),na.rm=TRUE) |
|
579 |
|
|
580 |
##Model assessment: general diagnostic/metrics |
|
581 |
##validation: using the testing data |
|
582 |
if (predval==1) { |
|
583 |
|
|
584 |
##Model assessment: specific diagnostic/metrics for GAM |
|
585 |
|
|
586 |
name<-paste("mod",j,sep="") #modj is the name of The "j" model (mod1 if j=1) |
|
587 |
mod<-get(name) #accessing GAM model ojbect "j" |
|
588 |
|
|
589 |
s_sgdf<-as(s_raster,"SpatialGridDataFrame") #Conversion to spatial grid data frame |
|
590 |
|
|
591 |
rpred<- predict(mod, newdata=s_sgdf, se.fit = TRUE) #Using the coeff to predict new values. |
|
592 |
y_pred<-rpred$fit #rpred is a list with fit being and array |
|
593 |
raster_pred<-r1 |
|
594 |
layerNames(raster_pred)<-"y_pred" |
|
595 |
values(raster_pred)<-as.numeric(y_pred) |
|
596 |
|
|
597 |
if (bias_prediction==1){ |
|
598 |
data_name<-paste("predicted_mod",j,"_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
599 |
"_",sampling_dat$run_samp[i],sep="") |
|
600 |
raster_name<-paste("GAM_bias_",data_name,out_prefix,".rst", sep="") |
|
601 |
writeRaster(raster_pred, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
602 |
bias_rast<-raster_pred |
|
603 |
|
|
604 |
raster_pred=themolst+daily_delta_rast-bias_rast #Final surface as a raster layer...wiht daily surface calculated earlier... |
|
605 |
layerNames(raster_pred)<-"y_pred" |
|
606 |
#=themolst+daily_delta_rast-bias_rast #Final surface as a raster layer... |
|
607 |
|
|
608 |
data_name<-paste("predicted_mod",j,"_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
609 |
"_",sampling_dat$run_samp[i],sep="") |
|
610 |
raster_name<-paste("GAM_bias_tmax_",data_name,out_prefix,".rst", sep="") |
|
611 |
writeRaster(raster_pred, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
612 |
|
|
613 |
} |
|
614 |
|
|
615 |
#rpred_val_s <- overlay(raster_pred,data_s) #This overlays the kriged surface tmax and the location of weather stations |
|
616 |
|
|
617 |
rpred_val_s <- overlay(pred_sgdf,data_s) #This overlays the interpolated surface tmax and the location of weather stations |
|
618 |
rpred_val_v <- overlay(pred_sgdf,data_v) #This overlays the interpolated surface tmax and the location of weather stations |
|
619 |
|
|
620 |
pred_mod<-paste("pred_mod",j,sep="") |
|
621 |
#Adding the results back into the original dataframes. |
|
622 |
data_s[[pred_mod]]<-rpred_val_s$y_pred |
|
623 |
|
|
624 |
data_v[[pred_mod]]<-rpred_val_v$y_pred |
|
625 |
|
|
626 |
#Model assessment: RMSE and then krig the residuals....! |
|
627 |
|
|
628 |
res_mod_s<-data_s[[y_var_name]] - data_s[[pred_mod]] #residuals from modeling training |
|
629 |
res_mod_v<-data_v[[y_var_name]] - data_v[[pred_mod]] #residuals from modeling validation |
|
630 |
|
|
631 |
} |
|
632 |
|
|
633 |
if (predval==0) { |
|
634 |
|
|
635 |
y_mod<- predict(mod, newdata=data_v, se.fit = TRUE) #Using the coeff to predict new values. |
|
636 |
|
|
637 |
pred_mod<-paste("pred_mod",j,sep="") |
|
638 |
#Adding the results back into the original dataframes. |
|
639 |
data_s[[pred_mod]]<-as.numeric(mod$fit) |
|
640 |
data_v[[pred_mod]]<-as.numeric(y_mod$fit) |
|
641 |
|
|
642 |
#Model assessment: RMSE and then krig the residuals....! |
|
643 |
|
|
644 |
#res_mod_s<- data_s$y_var - data_s[[pred_mod]] #Residuals from modeling training |
|
645 |
#res_mod_v<- data_v$y_var - data_v[[pred_mod]] #Residuals from modeling validation |
|
646 |
res_mod_s<-data_s[[y_var_name]] - data_s[[pred_mod]] |
|
647 |
res_mod_v<-data_v[[y_var_name]] - data_v[[pred_mod]] |
|
648 |
|
|
649 |
} |
|
650 |
|
|
651 |
####ADDED ON JULY 20th |
|
652 |
res_mod<-res_mod_v |
|
653 |
|
|
654 |
#RMSE_mod <- sqrt(sum(res_mod^2)/nv) #RMSE FOR REGRESSION STEP 1: GAM |
|
655 |
RMSE_mod<- sqrt(mean(res_mod^2,na.rm=TRUE)) |
|
656 |
#MAE_mod<- sum(abs(res_mod),na.rm=TRUE)/(nv-sum(is.na(res_mod))) #MAE from kriged surface validation |
|
657 |
MAE_mod<- mean(abs(res_mod), na.rm=TRUE) |
|
658 |
#ME_mod<- sum(res_mod,na.rm=TRUE)/(nv-sum(is.na(res_mod))) #ME, Mean Error or bias FOR REGRESSION STEP 1: GAM |
|
659 |
ME_mod<- mean(res_mod,na.rm=TRUE) #ME, Mean Error or bias FOR REGRESSION STEP 1: GAM |
|
660 |
#R2_mod<- cor(data_v$y_var,data_v[[pred_mod]])^2 #R2, coef. of var FOR REGRESSION STEP 1: GAM |
|
661 |
R2_mod<- cor(data_v$y_var,data_v[[pred_mod]], use="complete")^2 |
|
662 |
results_RMSE[1]<- sampling_dat$date[i] #storing the interpolation sampling_dat$date in the first column |
|
663 |
results_RMSE[2]<- ns #number of stations used in the training stage |
|
664 |
results_RMSE[3]<- "RMSE" |
|
665 |
results_RMSE[j+3]<- RMSE_mod #Storing RMSE for the model j |
|
666 |
results_MAE[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
667 |
results_MAE[2]<- ns #number of stations used in the training stage |
|
668 |
results_MAE[3]<- "MAE" |
|
669 |
results_MAE[j+3]<- MAE_mod #Storing MAE for the model j |
|
670 |
results_ME[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
671 |
results_ME[2]<- ns #number of stations used in the training stage |
|
672 |
results_ME[3]<- "ME" |
|
673 |
results_ME[j+3]<- ME_mod #Storing ME for the model j |
|
674 |
results_R2[1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
|
675 |
results_R2[2]<- ns #number of stations used in the training stage |
|
676 |
results_R2[3]<- "R2" |
|
677 |
results_R2[j+3]<- R2_mod #Storing R2 for the model j |
|
678 |
|
|
679 |
#Saving residuals and prediction in the dataframes: tmax predicted from GAM |
|
680 |
|
|
681 |
name2<-paste("res_mod",j,sep="") |
|
682 |
data_v[[name2]]<-as.numeric(res_mod_v) |
|
683 |
data_s[[name2]]<-as.numeric(res_mod_s) |
|
684 |
#end of loop calculating RMSE |
|
685 |
} |
|
273 |
data_name<-paste(y_var_name,"_predicted_",names(rast_clim_list)[j],"_", |
|
274 |
sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
275 |
"_",sampling_dat$run_samp[i],sep="") |
|
276 |
raster_name<-paste("fusion_",data_name,out_prefix,".tif", sep="") |
|
277 |
writeRaster(temp_predicted, filename=raster_name,overwrite=TRUE) |
|
278 |
temp_list[[j]]<-raster_name |
|
686 | 279 |
} |
687 | 280 |
|
688 |
#if (i==length(dates)){ |
|
689 |
|
|
690 |
#Specific diagnostic measures related to the testing datasets |
|
691 |
|
|
692 |
results_table_RMSE<-as.data.frame(results_RMSE) |
|
693 |
results_table_MAE<-as.data.frame(results_MAE) |
|
694 |
results_table_ME<-as.data.frame(results_ME) |
|
695 |
results_table_R2<-as.data.frame(results_R2) |
|
696 |
results_table_RMSE_f<-as.data.frame(results_RMSE_f) |
|
697 |
results_table_MAE_f<-as.data.frame(results_MAE_f) |
|
698 |
|
|
699 |
results_table_AIC<-as.data.frame(results_AIC) |
|
700 |
results_table_GCV<-as.data.frame(results_GCV) |
|
701 |
results_table_DEV<-as.data.frame(results_DEV) |
|
702 |
|
|
703 |
tb_metrics1<-rbind(results_table_RMSE,results_table_MAE, results_table_ME, results_table_R2,results_table_RMSE_f,results_table_MAE_f) # |
|
704 |
tb_metrics2<-rbind(results_table_AIC,results_table_GCV, results_table_DEV) |
|
705 |
|
|
706 |
#Preparing labels 10/30 |
|
707 |
mod_labels<-rep("mod",nmodels+1) |
|
708 |
index<-as.character(1:(nmodels+1)) |
|
709 |
mod_labels<-paste(mod_labels,index,sep="") |
|
710 |
cname<-c("dates","ns","metric", mod_labels) |
|
711 |
#cname<-c("dates","ns","metric","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7","mod8","mod9") |
|
712 |
colnames(tb_metrics1)<-cname |
|
713 |
#cname<-c("dates","ns","metric","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7","mod8") |
|
714 |
#colnames(tb_metrics2)<-cname |
|
715 |
colnames(tb_metrics2)<-cname[1:(nmodels+3)] |
|
716 |
|
|
717 |
#write.table(tb_diagnostic1, file= paste(path,"/","results_fusion_Assessment_measure1",out_prefix,".txt",sep=""), sep=",") |
|
718 |
|
|
719 |
#} |
|
720 |
print(paste(sampling_dat$date[i],"processed")) |
|
721 |
# end of the for loop1 |
|
722 |
#mod_obj<-list(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,mod9a,mod9b) |
|
723 |
#names(mod_obj)<-c("mod1","mod2","mod3","mod4","mod5","mod6","mod7","mod8","mod9a","mod9b") |
|
724 |
mod_labels_kr<-c("mod_kr_month", "mod_kr_day") |
|
725 |
names(mod_obj)<-c(mod_labels[1:nmodels],mod_labels_kr) |
|
726 |
results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2,mod_obj,data_month,list_formulas) |
|
727 |
names(results_list)<-c("data_s","data_v","tb_metrics1","tb_metrics2","mod_obj","data_month","formulas") |
|
728 |
|
|
729 |
#results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2) |
|
730 |
#results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2,mod_obj,sampling_dat[i,],data_month) |
|
731 |
#names(results_list)<-c("data_s","data_v","tb_metrics1","tb_metrics2","mod_obj","sampling_dat","data_month") |
|
732 |
save(results_list,file= paste(path,"/","results_list_metrics_objects_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
281 |
mod_krtmp2<-fitdelta |
|
282 |
model_name<-paste("mod_kr","day",sep="_") |
|
283 |
names(tmp_list)<-names(rast_clim_list) |
|
284 |
coordinates(data_s)<-cbind(data_s$x,data_s$y) |
|
285 |
proj4string(data_s)<-proj_str |
|
286 |
coordinates(data_v)<-cbind(data_v$x,data_v$y) |
|
287 |
proj4string(data_v)<-proj_str |
|
288 |
|
|
289 |
delta_obj<-list(temp_list,rast_clim_list,raster_name_delta,data_s, |
|
290 |
data_v,sampling_dat[i,],mod_krtmp2) |
|
291 |
|
|
292 |
obj_names<-c(y_var_name,"clim","delta","data_s","data_v", |
|
293 |
"sampling_dat",model_name) |
|
294 |
names(delta_obj)<-obj_names |
|
295 |
save(delta_obj,file= paste("delta_obj_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
733 | 296 |
"_",sampling_dat$run_samp[i],out_prefix,".RData",sep="")) |
734 |
return(results_list) |
|
735 |
#return(tb_diagnostic1) |
|
736 |
} |
|
297 |
return(delta_obj) |
|
298 |
|
|
299 |
} |
|
300 |
|
Also available in: Unified diff
GAM Fusion function major modification, separation of monhtly and daily steps, Venezuela prediction