Revision 55056785
Added by Benoit Parmentier almost 12 years ago
climate/research/oregon/interpolation/GAM_fusion_function_multisampling.R | ||
---|---|---|
8 | 8 |
|
9 | 9 |
#Adding layer LST to the raster stack |
10 | 10 |
|
11 |
pos<-match("LST",layerNames(s_raster)) #Find column with name "LST" |
|
12 |
s_raster<-dropLayer(s_raster,pos) # If it exists drop layer |
|
11 | 13 |
pos<-match(LST_month,layerNames(s_raster)) #Find column with the current month for instance mm12 |
12 | 14 |
r1<-raster(s_raster,layer=pos) #Select layer from stack |
13 | 15 |
layerNames(r1)<-"LST" |
14 |
s_raster<-addLayer(s_raster,r1) #Adding current month |
|
16 |
s_raster<-addLayer(s_raster,r1) #Adding current month as "LST"
|
|
15 | 17 |
|
16 | 18 |
###Regression part 1: Creating a validation dataset by creating training and testing datasets |
17 | 19 |
|
18 | 20 |
mod_LST <-ghcn.subsets[[i]][,match(LST_month, names(ghcn.subsets[[i]]))] #Match interpolation date and monthly LST average |
19 | 21 |
ghcn.subsets[[i]] = transform(ghcn.subsets[[i]],LST = mod_LST) #Add the variable LST to the subset dataset |
22 |
dst$LST<-dst[[LST_month]] |
|
20 | 23 |
#n<-nrow(ghcn.subsets[[i]]) |
21 | 24 |
#ns<-n-round(n*prop) #Create a sample from the data frame with 70% of the rows |
22 | 25 |
#nv<-n-ns #create a sample for validation with prop of the rows |
... | ... | |
48 | 51 |
# STEP 2 - Weather station means across same days: Monthly mean calculation |
49 | 52 |
########### |
50 | 53 |
|
51 |
modst=dst[dst$month==mo,] #Subsetting dataset for the relevant month of the date being processed
|
|
54 |
modst<-dst[dst$month==mo,] #Subsetting dataset for the relevant month of the date being processed
|
|
52 | 55 |
|
53 | 56 |
########## |
54 |
# STEP 3 - get LST at stations |
|
57 |
# STEP 3 - get LST at stations
|
|
55 | 58 |
########## |
56 | 59 |
|
57 | 60 |
sta_lola=modst[,c("lon","lat")] #Extracting locations of stations for the current month.. |
... | ... | |
94 | 97 |
names(d)[pos]<-c("id") |
95 | 98 |
names(x)[pos]<-c("id") |
96 | 99 |
names(modst)[1]<-c("id") #modst contains the average tmax per month for every stations... |
97 |
dmoday=merge(modst,d,by="id") #LOOSING DATA HERE!!! from 113 t0 103 |
|
98 |
xmoday=merge(modst,x,by="id") #LOOSING DATA HERE!!! from 48 t0 43 |
|
99 |
names(dmoday)[4]<-c("lat") |
|
100 |
names(dmoday)[5]<-c("lon") #dmoday contains all the the information: BIAS, monn |
|
101 |
names(xmoday)[4]<-c("lat") |
|
102 |
names(xmoday)[5]<-c("lon") #dmoday contains all the the information: BIAS, monn |
|
100 |
|
|
101 |
dmoday=merge(modst,d,by="id",suffixes=c("",".y2")) #LOOSING DATA HERE!!! from 113 t0 103 |
|
102 |
xmoday=merge(modst,x,by="id",suffixes=c("",".y2")) #LOOSING DATA HERE!!! from 48 t0 43 |
|
103 |
mod_pat<-glob2rx("*.y2") |
|
104 |
var_pat<-grep(mod_pat,names(dmoday),value=FALSE) # using grep with "value" extracts the matching names |
|
105 |
dmoday<-dmoday[,-var_pat] |
|
106 |
mod_pat<-glob2rx("*.y2") |
|
107 |
var_pat<-grep(mod_pat,names(xmoday),value=FALSE) # using grep with "value" extracts the matching names |
|
108 |
xmoday<-xmoday[,-var_pat] #Removing duplicate columns |
|
103 | 109 |
|
104 | 110 |
data_v<-xmoday |
105 | 111 |
### |
... | ... | |
107 | 113 |
#dmoday contains the daily tmax values for training with TMax being the monthly station tmax mean |
108 | 114 |
#xmoday contains the daily tmax values for validation with TMax being the monthly station tmax mean |
109 | 115 |
|
110 |
# windows() |
|
111 |
#png(paste("LST_TMax_scatterplot_",sampling_dat$date[i],out_prefix,".png", sep="")) |
|
112 | 116 |
png(paste("Daily_tmax_monthly_TMax_scatterplot_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
113 | 117 |
"_",sampling_dat$run_samp[i],out_prefix,".png", sep="")) |
114 | 118 |
plot(dailyTmax~TMax,data=dmoday,xlab="Mo Tmax",ylab=paste("Daily for",datelabel),main="across stations in OR") |
115 |
#savePlot(paste("Daily_tmax_monthly_TMax_scatterplot_",sampling_dat$date[i],out_prefix,".png", sep=""), type="png") |
|
116 |
#png(paste("LST_TMax_scatterplot_",sampling_dat$date[i],out_prefix,".png", sep="")) |
|
117 | 119 |
dev.off() |
118 | 120 |
|
119 | 121 |
######## |
120 | 122 |
# STEP 5 - interpolate bias |
121 | 123 |
######## |
122 | 124 |
|
123 |
# ?? include covariates like elev, distance to coast, cloud frequency, tree height |
|
124 |
#library(fields) |
|
125 |
#windows() |
|
126 |
quilt.plot(sta_lola,sta_bias,main="Bias at stations",asp=1) |
|
127 |
US(add=T,col="magenta",lwd=2) |
|
125 |
# ?? include covariates like elev, distance to coast, cloud frequency, tree heig |
|
128 | 126 |
#fitbias<-Tps(bias_xy,sta_bias) #use TPS or krige |
129 | 127 |
|
130 |
#Adding options to use only training stations: 07/11/2012 |
|
128 |
#Adding options to use only training stations : 07/11/2012
|
|
131 | 129 |
bias_xy=project(as.matrix(sta_lola),proj_str) |
132 | 130 |
#bias_xy2=project(as.matrix(c(dmoday$lon,dmoday$lat),proj_str) |
133 | 131 |
if(bias_val==1){ |
... | ... | |
138 | 136 |
fitbias<-Krig(bias_xy,sta_bias,theta=1e5) #use TPS or krige |
139 | 137 |
#The output is a krig object using fields |
140 | 138 |
mod9a<-fitbias |
141 |
# Creating plot of bias surface and saving it |
|
142 |
#X11() |
|
143 |
png(paste("Bias_surface_LST_TMax_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
144 |
"_",sampling_dat$run_samp[i],out_prefix,".png", sep="")) #Create file to write a plot |
|
145 |
datelabel2=format(ISOdate(year,mo,day),"%B ") #added by Benoit, label |
|
146 |
surface(fitbias,col=rev(terrain.colors(100)),asp=1,main=paste("Interpolated bias for",datelabel2,sep=" ")) #Plot to file |
|
147 |
#savePlot(paste("Bias_surface_LST_TMax_",sampling_dat$date[i],out_prefix,".png", sep=""), type="png") |
|
148 |
dev.off() #Release the hold to the file |
|
149 |
|
|
150 |
#US(add=T,col="magenta",lwd=2) |
|
151 |
|
|
139 |
|
|
152 | 140 |
########## |
153 | 141 |
# STEP 7 - interpolate delta across space |
154 | 142 |
########## |
... | ... | |
156 | 144 |
daily_sta_lola=dmoday[,c("lon","lat")] #could be same as before but why assume merge does this - assume not |
157 | 145 |
daily_sta_xy=project(as.matrix(daily_sta_lola),proj_str) |
158 | 146 |
daily_delta=dmoday$dailyTmax-dmoday$TMax |
159 |
#windows() |
|
160 |
quilt.plot(daily_sta_lola,daily_delta,asp=1,main="Station delta for Jan 15") |
|
161 |
US(add=T,col="magenta",lwd=2) |
|
147 |
|
|
162 | 148 |
#fitdelta<-Tps(daily_sta_xy,daily_delta) #use TPS or krige |
163 | 149 |
fitdelta<-Krig(daily_sta_xy,daily_delta,theta=1e5) #use TPS or krige |
164 | 150 |
#Kriging using fields package |
165 | 151 |
mod9b<-fitdelta |
166 |
# Creating plot of bias surface and saving it |
|
167 |
#X11() |
|
152 |
|
|
168 | 153 |
png(paste("Delta_surface_LST_TMax_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
169 | 154 |
"_",sampling_dat$run_samp[i],out_prefix,".png", sep="")) |
170 | 155 |
surface(fitdelta,col=rev(terrain.colors(100)),asp=1,main=paste("Interpolated delta for",datelabel,sep=" ")) |
171 |
#savePlot(paste("Delta_surface_LST_TMax_",sampling_dat$date[i],out_prefix,".png", sep=""), type="png") |
|
172 | 156 |
dev.off() |
173 |
#US(add=T,col="magenta",lwd=2) |
|
174 | 157 |
# |
175 | 158 |
|
176 | 159 |
#### Added by Benoit on 06/19 |
... | ... | |
237 | 220 |
mae<-MAE_fun(sta_pred,tmax) |
238 | 221 |
r2<-cor(sta_pred,tmax)^2 #R2, coef. of var |
239 | 222 |
me<-mean(sta_pred-tmax) |
240 |
|
|
241 |
#plot(sta_pred~dmoday$dailyTmax,xlab=paste("Actual daily for",datelabel),ylab="Pred daily",main=paste("RMSE=",rmse)) |
|
242 |
|
|
223 |
|
|
243 | 224 |
png(paste("Predicted_tmax_versus_observed_scatterplot_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
244 | 225 |
"_",sampling_dat$run_samp[i],out_prefix,".png", sep="")) |
245 | 226 |
plot(sta_pred~tmax,xlab=paste("Actual daily for",datelabel),ylab="Pred daily",main=paste("RMSE=",rmse)) |
246 | 227 |
abline(0,1) |
247 |
#savePlot(paste("Predicted_tmax_versus_observed_scatterplot_",sampling_dat$date[i],out_prefix,".png", sep=""), type="png") |
|
248 | 228 |
dev.off() |
249 | 229 |
#resid=sta_pred-dmoday$dailyTmax |
250 | 230 |
resid=sta_pred-tmax |
251 |
#quilt.plot(daily_sta_lola,resid) |
|
252 |
|
|
253 |
### END OF BRIAN's code |
|
254 |
|
|
255 |
### Added by benoit |
|
256 | 231 |
|
257 | 232 |
###BEFORE GAM prediction the data object must be transformed to SDF |
258 | 233 |
|
... | ... | |
262 | 237 |
coords<- data_s[,c('x_OR83M','y_OR83M')] |
263 | 238 |
coordinates(data_s)<-coords |
264 | 239 |
proj4string(data_s)<-CRS #Need to assign coordinates.. |
240 |
coords<- modst[,c('x_OR83M','y_OR83M')] |
|
241 |
coordinates(modst)<-coords |
|
242 |
proj4string(modst)<-CRS #Need to assign coordinates.. |
|
265 | 243 |
|
266 | 244 |
ns<-nrow(data_s) #This is added to because some loss of data might have happened because of the averaging... |
267 | 245 |
nv<-nrow(data_v) |
268 | 246 |
|
269 | 247 |
###GAM PREDICTION |
270 | 248 |
|
271 |
#data_s$y_var<-data_s$dailyTmax #This shoudl be changed for any variable!!! |
|
272 |
#data_v$y_var<-data_v$dailyTmax |
|
273 |
data_v$y_var<-data_v[[y_var_name]] |
|
274 |
data_s$y_var<-data_s[[y_var_name]] |
|
249 |
if (bias_prediction==1){ |
|
250 |
data_s$y_var<-data_s$LSTD_bias #This shoudl be changed for any variable!!! |
|
251 |
data_v$y_var<-data_v$LSTD_bias |
|
252 |
data_month<-modst |
|
253 |
data_month$y_var<-modst$LSTD_bias |
|
254 |
} |
|
255 |
|
|
256 |
if (bias_prediction==0){ |
|
257 |
data_v$y_var<-data_v[[y_var_name]] |
|
258 |
data_s$y_var<-data_s[[y_var_name]] |
|
259 |
} |
|
275 | 260 |
|
276 | 261 |
#Model and response variable can be changed without affecting the script |
277 | 262 |
|
... | ... | |
282 | 267 |
formula5 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)", env=.GlobalEnv) |
283 | 268 |
formula6 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC1)", env=.GlobalEnv) |
284 | 269 |
formula7 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC3)", env=.GlobalEnv) |
285 |
formula8 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST) + s(LST,LC1)", env=.GlobalEnv) |
|
286 |
|
|
287 |
#formula1 <- as.formula("y_var ~ s(lat,lon,ELEV_SRTM)", env=.GlobalEnv) |
|
288 |
#formula2 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST) + s(CANHEIGHT)", env=.GlobalEnv) |
|
289 |
#formula3 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST,CANHEIGHT)", env=.GlobalEnv) |
|
290 |
#formula4 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST,LC1)", env=.GlobalEnv) |
|
291 |
#formula5 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST,LC3)", env=.GlobalEnv) |
|
292 |
#formula6 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC1)", env=.GlobalEnv) |
|
293 |
#formula7 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC3)", env=.GlobalEnv) |
|
294 |
#formula8 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST,LC1) + s(LST,LC3)", env=.GlobalEnv) |
|
295 |
|
|
296 |
mod1<- try(gam(formula1, data=data_s)) |
|
297 |
mod2<- try(gam(formula2, data=data_s)) #modified nesting....from 3 to 2 |
|
298 |
mod3<- try(gam(formula3, data=data_s)) |
|
299 |
mod4<- try(gam(formula4, data=data_s)) |
|
300 |
mod5<- try(gam(formula5, data=data_s)) |
|
301 |
mod6<- try(gam(formula6, data=data_s)) |
|
302 |
mod7<- try(gam(formula7, data=data_s)) |
|
303 |
mod8<- try(gam(formula8, data=data_s)) |
|
304 |
|
|
305 |
# mod1<- try(gam(formula1, data=data_s)) |
|
306 |
# mod2<- try(gam(formula2, data=data_s)) #modified nesting....from 3 to 2 |
|
307 |
# mod3<- try(gam(y_var~ s(lat) + s (lon) + s (ELEV_SRTM) + s (Northness)+ s (Eastness) + s(DISTOC), data=data_s)) |
|
308 |
# mod4<- try(gam(y_var~ s(lat) + s (lon) + s(ELEV_SRTM) + s(Northness) + s (Eastness) + s(DISTOC) + s(LST), data=data_s)) |
|
309 |
# mod5<- try(gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST), data=data_s)) |
|
310 |
# mod6<- try(gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC1), data=data_s)) |
|
311 |
# mod7<- try(gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC3), data=data_s)) |
|
312 |
# mod8<- try(gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST) + s(LC1,LC3), data=data_s)) |
|
313 |
# |
|
270 |
formula8 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST) + s(LC1,LC3)", env=.GlobalEnv) |
|
271 |
|
|
272 |
if (bias_prediction==1){ |
|
273 |
mod1<- try(gam(formula1, data=data_month)) |
|
274 |
mod2<- try(gam(formula2, data=data_month)) #modified nesting....from 3 to 2 |
|
275 |
mod3<- try(gam(formula3, data=data_month)) |
|
276 |
mod4<- try(gam(formula4, data=data_month)) |
|
277 |
mod5<- try(gam(formula5, data=data_month)) |
|
278 |
mod6<- try(gam(formula6, data=data_month)) |
|
279 |
mod7<- try(gam(formula7, data=data_month)) |
|
280 |
mod8<- try(gam(formula8, data=data_month)) |
|
281 |
|
|
282 |
} else if (bias_prediction==0){ |
|
283 |
|
|
284 |
mod1<- try(gam(formula1, data=data_s)) |
|
285 |
mod2<- try(gam(formula2, data=data_s)) #modified nesting....from 3 to 2 |
|
286 |
mod3<- try(gam(formula3, data=data_s)) |
|
287 |
mod4<- try(gam(formula4, data=data_s)) |
|
288 |
mod5<- try(gam(formula5, data=data_s)) |
|
289 |
mod6<- try(gam(formula6, data=data_s)) |
|
290 |
mod7<- try(gam(formula7, data=data_s)) |
|
291 |
mod8<- try(gam(formula8, data=data_s)) |
|
292 |
} |
|
293 |
|
|
314 | 294 |
#Added |
315 | 295 |
#tmax_predicted=themolst+daily_delta_rast-bias_rast #Final surface?? but daily_rst |
316 | 296 |
|
... | ... | |
351 | 331 |
#ns<-nrow(data_s) #This is added to because some loss of data might have happened because of the averaging... |
352 | 332 |
#nv<-nrow(data_v) |
353 | 333 |
|
334 |
pred_mod<-paste("pred_mod",j,sep="") |
|
335 |
#Adding the results back into the original dataframes. |
|
336 |
data_s[[pred_mod]]<-sta_pred_data_s |
|
337 |
data_v[[pred_mod]]<-sta_pred |
|
338 |
|
|
339 |
#Model assessment: RMSE and then krig the residuals....! |
|
340 |
|
|
341 |
res_mod_s<- data_s$dailyTmax - data_s[[pred_mod]] #Residuals from kriging training |
|
342 |
res_mod_v<- data_v$dailyTmax - data_v[[pred_mod]] #Residuals from kriging validation |
|
343 |
|
|
344 |
name2<-paste("res_mod",j,sep="") |
|
345 |
data_v[[name2]]<-as.numeric(res_mod_v) |
|
346 |
data_s[[name2]]<-as.numeric(res_mod_s) |
|
347 |
|
|
354 | 348 |
|
355 | 349 |
for (j in 1:nmodels){ |
356 | 350 |
|
... | ... | |
453 | 447 |
s_sgdf<-as(s_raster,"SpatialGridDataFrame") #Conversion to spatial grid data frame |
454 | 448 |
|
455 | 449 |
rpred<- predict(mod, newdata=s_sgdf, se.fit = TRUE) #Using the coeff to predict new values. |
456 |
y_pred<-rpred$fit |
|
450 |
y_pred<-rpred$fit #rpred is a list with fit being and array
|
|
457 | 451 |
raster_pred<-r1 |
458 | 452 |
layerNames(raster_pred)<-"y_pred" |
459 | 453 |
values(raster_pred)<-as.numeric(y_pred) |
460 |
data_name<-paste("predicted_mod",j,"_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
461 |
"_",sampling_dat$run_samp[i],sep="") |
|
462 |
raster_name<-paste("GAM_",data_name,out_prefix,".rst", sep="") |
|
463 |
writeRaster(raster_pred, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
464 |
#writeRaster(r2, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
454 |
|
|
455 |
if (bias_prediction==1){ |
|
456 |
data_name<-paste("predicted_mod",j,"_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
457 |
"_",sampling_dat$run_samp[i],sep="") |
|
458 |
raster_name<-paste("GAM_bias_",data_name,out_prefix,".rst", sep="") |
|
459 |
writeRaster(raster_pred, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
460 |
bias_rast<-raster_pred |
|
461 |
|
|
462 |
raster_pred=themolst+daily_delta_rast-bias_rast #Final surface as a raster layer... |
|
463 |
layerNames(raster_pred)<-"y_pred" |
|
464 |
#=themolst+daily_delta_rast-bias_rast #Final surface as a raster layer... |
|
465 |
|
|
466 |
data_name<-paste("predicted_mod",j,"_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
467 |
"_",sampling_dat$run_samp[i],sep="") |
|
468 |
raster_name<-paste("GAM_bias_tmax_",data_name,out_prefix,".rst", sep="") |
|
469 |
writeRaster(raster_pred, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
470 |
|
|
471 |
} |
|
472 |
|
|
473 |
if (bias_prediction==0){ |
|
474 |
data_name<-paste("predicted_mod",j,"_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
475 |
"_",sampling_dat$run_samp[i],sep="") |
|
476 |
raster_name<-paste("GAM_",data_name,out_prefix,".rst", sep="") |
|
477 |
writeRaster(raster_pred, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
478 |
#writeRaster(r2, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
|
479 |
|
|
480 |
} |
|
481 |
|
|
465 | 482 |
|
466 | 483 |
pred_sgdf<-as(raster_pred,"SpatialGridDataFrame") #Conversion to spatial grid data frame |
467 | 484 |
#rpred_val_s <- overlay(raster_pred,data_s) #This overlays the kriged surface tmax and the location of weather stations |
468 | 485 |
|
469 |
rpred_val_s <- overlay(pred_sgdf,data_s) #This overlays the kriged surface tmax and the location of weather stations
|
|
470 |
rpred_val_v <- overlay(pred_sgdf,data_v) #This overlays the kriged surface tmax and the location of weather stations
|
|
486 |
rpred_val_s <- overlay(pred_sgdf,data_s) #This overlays the interpolated surface tmax and the location of weather stations
|
|
487 |
rpred_val_v <- overlay(pred_sgdf,data_v) #This overlays the interpolated surface tmax and the location of weather stations
|
|
471 | 488 |
|
472 | 489 |
pred_mod<-paste("pred_mod",j,sep="") |
473 | 490 |
#Adding the results back into the original dataframes. |
... | ... | |
477 | 494 |
|
478 | 495 |
#Model assessment: RMSE and then krig the residuals....! |
479 | 496 |
|
480 |
res_mod_s<- data_s$y_var - data_s[[pred_mod]] #Residuals from kriging training
|
|
481 |
res_mod_v<- data_v$y_var - data_v[[pred_mod]] #Residuals from kriging validation
|
|
497 |
res_mod_s<-data_s[[y_var_name]] - data_s[[pred_mod]] #residuals from modeling training
|
|
498 |
res_mod_v<-data_v[[y_var_name]] - data_v[[pred_mod]] #residuals from modeling validation
|
|
482 | 499 |
|
483 | 500 |
} |
484 | 501 |
|
... | ... | |
493 | 510 |
|
494 | 511 |
#Model assessment: RMSE and then krig the residuals....! |
495 | 512 |
|
496 |
res_mod_s<- data_s$y_var - data_s[[pred_mod]] #Residuals from kriging training |
|
497 |
res_mod_v<- data_v$y_var - data_v[[pred_mod]] #Residuals from kriging validation |
|
513 |
#res_mod_s<- data_s$y_var - data_s[[pred_mod]] #Residuals from modeling training |
|
514 |
#res_mod_v<- data_v$y_var - data_v[[pred_mod]] #Residuals from modeling validation |
|
515 |
res_mod_s<-data_s[[y_var_name]] - data_s[[pred_mod]] |
|
516 |
res_mod_v<-data_v[[y_var_name]] - data_v[[pred_mod]] |
|
517 |
|
|
498 | 518 |
} |
499 | 519 |
|
500 | 520 |
####ADDED ON JULY 20th |
... | ... | |
568 | 588 |
mod_obj<-list(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,mod9a,mod9b) |
569 | 589 |
names(mod_obj)<-c("mod1","mod2","mod3","mod4","mod5","mod6","mod7","mod8","mod9a","mod9b") |
570 | 590 |
#results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2) |
571 |
results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2,mod_obj) |
|
572 |
names(results_list)<-c("data_s","data_v","tb_metrics1","tb_metrics2","mod_obj") |
|
591 |
results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2,mod_obj,sampling_dat[i,],data_month)
|
|
592 |
names(results_list)<-c("data_s","data_v","tb_metrics1","tb_metrics2","mod_obj","sampling_dat","data_month")
|
|
573 | 593 |
save(results_list,file= paste(path,"/","results_list_metrics_objects_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
574 | 594 |
"_",sampling_dat$run_samp[i],out_prefix,".RData",sep="")) |
575 | 595 |
return(results_list) |
Also available in: Unified diff
GAM fusion function, added GAM models for bias surface and extraction of monthly mean tmax OR