Revision d93bea39
Added by Benoit Parmentier almost 12 years ago
| climate/research/oregon/interpolation/GAM_fusion_function_multisampling.R | ||
|---|---|---|
| 1 | 1 |
runGAMFusion <- function(i) { # loop over dates
|
| 2 | 2 |
|
| 3 |
#date<-strptime(dates[i], "%Y%m%d") # interpolation date being processed |
|
| 4 | 3 |
date<-strptime(sampling_dat$date[i], "%Y%m%d") # interpolation date being processed |
| 5 |
|
|
| 6 | 4 |
month<-strftime(date, "%m") # current month of the date being processed |
| 7 | 5 |
LST_month<-paste("mm_",month,sep="") # name of LST month to be matched
|
| 8 | 6 |
|
| 9 | 7 |
#Adding layer LST to the raster stack |
| 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 |
|
| 8 |
|
|
| 9 |
pos<-match("LST",layerNames(s_raster)) #Find the position of the layer with name "LST", if not present pos=NA
|
|
| 10 |
s_raster<-dropLayer(s_raster,pos) # If it exists drop layer
|
|
| 13 | 11 |
pos<-match(LST_month,layerNames(s_raster)) #Find column with the current month for instance mm12 |
| 14 | 12 |
r1<-raster(s_raster,layer=pos) #Select layer from stack |
| 15 | 13 |
layerNames(r1)<-"LST" |
| 16 |
s_raster<-addLayer(s_raster,r1) #Adding current month as "LST" |
|
| 14 |
#Screen for extreme values" 10/30 |
|
| 15 |
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) |
|
| 16 |
r1[r1 < (min_val)]<-NA |
|
| 17 |
s_raster<-addLayer(s_raster,r1) #Adding current month |
|
| 17 | 18 |
|
| 18 | 19 |
###Regression part 1: Creating a validation dataset by creating training and testing datasets |
| 19 | 20 |
|
| ... | ... | |
| 45 | 46 |
########### |
| 46 | 47 |
|
| 47 | 48 |
themolst<-raster(molst,mo) #current month being processed saved in a raster image |
| 49 |
min_val<-(-15) #Screening for extreme values |
|
| 50 |
themolst[themolst < (min_val)]<-NA |
|
| 51 |
|
|
| 48 | 52 |
plot(themolst) |
| 49 | 53 |
|
| 50 | 54 |
########### |
| ... | ... | |
| 134 | 138 |
} |
| 135 | 139 |
|
| 136 | 140 |
fitbias<-Krig(bias_xy,sta_bias,theta=1e5) #use TPS or krige |
| 137 |
#The output is a krig object using fields |
|
| 138 |
mod9a<-fitbias |
|
| 141 |
#The output is a krig object using fields: modif 10/30 |
|
| 142 |
#mod9a<-fitbias |
|
| 143 |
mod_krtmp1<-fitbias |
|
| 144 |
model_name<-paste("mod_kr","month",sep="_")
|
|
| 145 |
assign(model_name,mod_krtmp1) |
|
| 146 |
|
|
| 139 | 147 |
|
| 140 | 148 |
########## |
| 141 | 149 |
# STEP 7 - interpolate delta across space |
| ... | ... | |
| 147 | 155 |
|
| 148 | 156 |
#fitdelta<-Tps(daily_sta_xy,daily_delta) #use TPS or krige |
| 149 | 157 |
fitdelta<-Krig(daily_sta_xy,daily_delta,theta=1e5) #use TPS or krige |
| 150 |
#Kriging using fields package |
|
| 151 |
mod9b<-fitdelta |
|
| 152 |
|
|
| 158 |
#Kriging using fields package: modif 10/30 |
|
| 159 |
#mod9b<-fitdelta |
|
| 160 |
mod_krtmp2<-fitdelta |
|
| 161 |
model_name<-paste("mod_kr","day",sep="_")
|
|
| 162 |
assign(model_name,mod_krtmp2) |
|
| 163 |
|
|
| 153 | 164 |
png(paste("Delta_surface_LST_TMax_",sampling_dat$date[i],"_",sampling_dat$prop[i],
|
| 154 | 165 |
"_",sampling_dat$run_samp[i],out_prefix,".png", sep="")) |
| 155 | 166 |
surface(fitdelta,col=rev(terrain.colors(100)),asp=1,main=paste("Interpolated delta for",datelabel,sep=" "))
|
| ... | ... | |
| 260 | 271 |
|
| 261 | 272 |
#Model and response variable can be changed without affecting the script |
| 262 | 273 |
|
| 263 |
formula1 <- as.formula("y_var ~ s(lat) + s(lon) + s(ELEV_SRTM)", env=.GlobalEnv)
|
|
| 264 |
formula2 <- as.formula("y_var~ s(lat,lon)+ s(ELEV_SRTM)", env=.GlobalEnv)
|
|
| 265 |
formula3 <- as.formula("y_var~ s(lat) + s (lon) + s (ELEV_SRTM) + s (Northness)+ s (Eastness) + s(DISTOC)", env=.GlobalEnv)
|
|
| 266 |
formula4 <- as.formula("y_var~ s(lat) + s (lon) + s(ELEV_SRTM) + s(Northness) + s (Eastness) + s(DISTOC) + s(LST)", env=.GlobalEnv)
|
|
| 267 |
formula5 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)", env=.GlobalEnv)
|
|
| 268 |
formula6 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC1)", env=.GlobalEnv)
|
|
| 269 |
formula7 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC3)", env=.GlobalEnv)
|
|
| 270 |
formula8 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST) + s(LC1,LC3)", env=.GlobalEnv)
|
|
| 274 |
list_formulas<-vector("list",nmodels)
|
|
| 275 |
|
|
| 276 |
list_formulas[[1]] <- as.formula("y_var ~ s(ELEV_SRTM)", env=.GlobalEnv)
|
|
| 277 |
list_formulas[[2]] <- as.formula("y_var ~ s(LST)", env=.GlobalEnv)
|
|
| 278 |
list_formulas[[3]] <- as.formula("y_var ~ s(ELEV_SRTM,LST)", env=.GlobalEnv)
|
|
| 279 |
list_formulas[[4]] <- as.formula("y_var ~ s(lat) + s(lon)+ s(ELEV_SRTM)", env=.GlobalEnv)
|
|
| 280 |
list_formulas[[5]] <- as.formula("y_var ~ s(lat,lon,ELEV_SRTM)", env=.GlobalEnv)
|
|
| 281 |
list_formulas[[6]] <- as.formula("y_var ~ s(lat,lon) + s(ELEV_SRTM) + s(Northness_w,Eastness_w) + s(LST)", env=.GlobalEnv)
|
|
| 282 |
list_formulas[[7]] <- as.formula("y_var ~ s(lat,lon) + s(ELEV_SRTM) + s(Northness_w,Eastness_w) + s(LST) + s(LC1)", env=.GlobalEnv)
|
|
| 283 |
list_formulas[[8]] <- as.formula("y_var ~ s(lat,lon) + s(ELEV_SRTM) + s(Northness_w,Eastness_w) + s(LST) + s(LC3)", env=.GlobalEnv)
|
|
| 284 |
list_formulas[[9]] <- as.formula("y_var ~ s(lat,lon) + s(ELEV_SRTM) + s(Northness_w,Eastness_w) + s(LST) + s(DISTOC)", env=.GlobalEnv)
|
|
| 285 |
|
|
| 286 |
#list_formulas[[1]] <- as.formula("y_var ~ s(ELEV_SRTM)", env=.GlobalEnv)
|
|
| 287 |
#list_formulas[[2]] <- as.formula("y_var ~ s(lat,lon)", env=.GlobalEnv)
|
|
| 288 |
#list_formulas[[3]] <- as.formula("y_var~ s(lat,lon,ELEV_SRTM)", env=.GlobalEnv)
|
|
| 289 |
#list_formulas[[4]] <- as.formula("y_var~ s(lat) + s (lon) + s (ELEV_SRTM) + s(DISTOC)", env=.GlobalEnv)
|
|
| 290 |
#list_formulas[[5]] <- as.formula("y_var~ s(lat,lon,ELEV_SRTM) + s(Northness) + s (Eastness) + s(DISTOC)", env=.GlobalEnv)
|
|
| 291 |
#list_formulas[[6]] <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC)", env=.GlobalEnv)
|
|
| 271 | 292 |
|
| 272 | 293 |
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)) |
|
| 294 |
#mod1<- try(gam(formula1, data=data_month))
|
|
| 295 |
#mod2<- try(gam(formula2, data=data_month)) #modified nesting....from 3 to 2
|
|
| 296 |
#mod3<- try(gam(formula3, data=data_month))
|
|
| 297 |
#mod4<- try(gam(formula4, data=data_month))
|
|
| 298 |
#mod5<- try(gam(formula5, data=data_month))
|
|
| 299 |
#mod6<- try(gam(formula6, data=data_month))
|
|
| 300 |
#mod7<- try(gam(formula7, data=data_month))
|
|
| 301 |
#mod8<- try(gam(formula8, data=data_month))
|
|
| 281 | 302 |
|
| 282 |
} else if (bias_prediction==0){
|
|
| 303 |
for (j in 1:nmodels){
|
|
| 304 |
formula<-list_formulas[[j]] |
|
| 305 |
mod<- try(gam(formula, data=data_month)) |
|
| 306 |
model_name<-paste("mod",j,sep="")
|
|
| 307 |
assign(model_name,mod) |
|
| 308 |
} |
|
| 309 |
|
|
| 310 |
} else if (bias_prediction==0){ #Use daily data for direct prediction using GAM
|
|
| 311 |
|
|
| 312 |
#mod1<- try(gam(formula1, data=data_s)) |
|
| 313 |
#mod2<- try(gam(formula2, data=data_s)) #modified nesting....from 3 to 2 |
|
| 314 |
#mod3<- try(gam(formula3, data=data_s)) |
|
| 315 |
#mod4<- try(gam(formula4, data=data_s)) |
|
| 316 |
#mod5<- try(gam(formula5, data=data_s)) |
|
| 317 |
#mod6<- try(gam(formula6, data=data_s)) |
|
| 318 |
#mod7<- try(gam(formula7, data=data_s)) |
|
| 319 |
#mod8<- try(gam(formula8, data=data_s)) |
|
| 320 |
|
|
| 321 |
for (j in 1:nmodels){
|
|
| 322 |
formula<-list_formulas[[j]] |
|
| 323 |
mod<- try(gam(formula, data=data_s)) |
|
| 324 |
model_name<-paste("mod",j,sep="")
|
|
| 325 |
assign(model_name,mod) |
|
| 326 |
} |
|
| 283 | 327 |
|
| 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 | 328 |
} |
| 293 | 329 |
|
| 294 | 330 |
#Added |
| ... | ... | |
| 345 | 381 |
data_v[[name2]]<-as.numeric(res_mod_v) |
| 346 | 382 |
data_s[[name2]]<-as.numeric(res_mod_s) |
| 347 | 383 |
|
| 384 |
mod_obj<-vector("list",nmodels+2) #This will contain the model objects fitting: 10/30
|
|
| 385 |
mod_obj[[nmodels+1]]<-mod_kr_month #Storing climatology object |
|
| 386 |
mod_obj[[nmodels+2]]<-mod_kr_day #Storing delta object |
|
| 348 | 387 |
|
| 349 | 388 |
for (j in 1:nmodels){
|
| 350 | 389 |
|
| ... | ... | |
| 352 | 391 |
|
| 353 | 392 |
name<-paste("mod",j,sep="") #modj is the name of The "j" model (mod1 if j=1)
|
| 354 | 393 |
mod<-get(name) #accessing GAM model ojbect "j" |
| 394 |
mod_obj[[j]]<-mod #storing current model object |
|
| 355 | 395 |
|
| 356 | 396 |
#If mod "j" is not a model object |
| 357 | 397 |
if (inherits(mod,"try-error")) {
|
| ... | ... | |
| 459 | 499 |
writeRaster(raster_pred, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
| 460 | 500 |
bias_rast<-raster_pred |
| 461 | 501 |
|
| 462 |
raster_pred=themolst+daily_delta_rast-bias_rast #Final surface as a raster layer... |
|
| 502 |
raster_pred=themolst+daily_delta_rast-bias_rast #Final surface as a raster layer...wiht daily surface calculated earlier...
|
|
| 463 | 503 |
layerNames(raster_pred)<-"y_pred" |
| 464 | 504 |
#=themolst+daily_delta_rast-bias_rast #Final surface as a raster layer... |
| 465 | 505 |
|
| ... | ... | |
| 571 | 611 |
|
| 572 | 612 |
tb_metrics1<-rbind(results_table_RMSE,results_table_MAE, results_table_ME, results_table_R2,results_table_RMSE_f,results_table_MAE_f) # |
| 573 | 613 |
tb_metrics2<-rbind(results_table_AIC,results_table_GCV, results_table_DEV) |
| 574 |
cname<-c("dates","ns","metric","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7","mod8","mod9")
|
|
| 614 |
|
|
| 615 |
#Preparing labels 10/30 |
|
| 616 |
mod_labels<-rep("mod",nmodels+1)
|
|
| 617 |
index<-as.character(1:(nmodels+1)) |
|
| 618 |
mod_labels<-paste(mod_labels,index,sep="") |
|
| 619 |
cname<-c("dates","ns","metric", mod_labels)
|
|
| 620 |
#cname<-c("dates","ns","metric","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7","mod8","mod9")
|
|
| 575 | 621 |
colnames(tb_metrics1)<-cname |
| 576 |
cname<-c("dates","ns","metric","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7","mod8")
|
|
| 577 |
colnames(tb_metrics2)<-cname |
|
| 578 |
#colnames(results_table_RMSE)<-cname |
|
| 579 |
#colnames(results_table_RMSE_f)<-cname |
|
| 580 |
#tb_diagnostic1<-results_table_RMSE #measures of validation |
|
| 581 |
#tb_diagnostic2<-results_table_RMSE_f #measures of fit |
|
| 622 |
#cname<-c("dates","ns","metric","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7","mod8")
|
|
| 623 |
#colnames(tb_metrics2)<-cname |
|
| 624 |
colnames(tb_metrics2)<-cname[1:(nmodels+3)] |
|
| 582 | 625 |
|
| 583 | 626 |
#write.table(tb_diagnostic1, file= paste(path,"/","results_fusion_Assessment_measure1",out_prefix,".txt",sep=""), sep=",") |
| 584 | 627 |
|
| 585 | 628 |
#} |
| 586 | 629 |
print(paste(sampling_dat$date[i],"processed")) |
| 587 | 630 |
# end of the for loop1 |
| 588 |
mod_obj<-list(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,mod9a,mod9b) |
|
| 589 |
names(mod_obj)<-c("mod1","mod2","mod3","mod4","mod5","mod6","mod7","mod8","mod9a","mod9b")
|
|
| 631 |
#mod_obj<-list(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,mod9a,mod9b) |
|
| 632 |
#names(mod_obj)<-c("mod1","mod2","mod3","mod4","mod5","mod6","mod7","mod8","mod9a","mod9b")
|
|
| 633 |
mod_labels_kr<-c("mod_kr_month", "mod_kr_day")
|
|
| 634 |
names(mod_obj)<-c(mod_labels[1:nmodels],mod_labels_kr) |
|
| 635 |
results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2,mod_obj,data_month,list_formulas) |
|
| 636 |
names(results_list)<-c("data_s","data_v","tb_metrics1","tb_metrics2","mod_obj","data_month","formulas")
|
|
| 637 |
|
|
| 590 | 638 |
#results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2) |
| 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")
|
|
| 639 |
#results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2,mod_obj,sampling_dat[i,],data_month)
|
|
| 640 |
#names(results_list)<-c("data_s","data_v","tb_metrics1","tb_metrics2","mod_obj","sampling_dat","data_month")
|
|
| 593 | 641 |
save(results_list,file= paste(path,"/","results_list_metrics_objects_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
| 594 | 642 |
"_",sampling_dat$run_samp[i],out_prefix,".RData",sep="")) |
| 595 | 643 |
return(results_list) |
Also available in: Unified diff
GAM fusion function, IBS 2013 models and additional cleaning of code