Revision c851123d
Added by Benoit Parmentier almost 12 years ago
| climate/research/oregon/interpolation/GAM_CAI_function_multisampling.R | ||
|---|---|---|
| 1 | 1 |
runGAMCAI <- function(i) { # loop over dates
|
| 2 | 2 |
|
| 3 |
#With upates from 10/26/2012 |
|
| 3 | 4 |
#date<-strptime(dates[i], "%Y%m%d") # interpolation date being processed |
| 4 | 5 |
date<-strptime(sampling_dat$date[i], "%Y%m%d") # interpolation date being processed, converting the string using specific format |
| 5 | 6 |
month<-strftime(date, "%m") # current month of the date being processed |
| ... | ... | |
| 12 | 13 |
pos<-match(LST_month,layerNames(s_raster)) #Find column with the current month for instance mm12 |
| 13 | 14 |
r1<-raster(s_raster,layer=pos) #Select layer from stack |
| 14 | 15 |
layerNames(r1)<-"LST" |
| 16 |
#Screen for extreme values" 10/30 |
|
| 17 |
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) |
|
| 18 |
r1[r1 < (min_val)]<-NA |
|
| 15 | 19 |
s_raster<-addLayer(s_raster,r1) #Adding current month |
| 16 | 20 |
|
| 17 | 21 |
###Regression part 1: Creating a validation dataset by creating training and testing datasets |
| ... | ... | |
| 45 | 49 |
########### |
| 46 | 50 |
|
| 47 | 51 |
themolst<-raster(molst,mo) #current month being processed saved in a raster image |
| 52 |
min_val<-(-15) #Screening for extreme values |
|
| 53 |
themolst[themolst < (min_val)]<-NA |
|
| 54 |
|
|
| 48 | 55 |
plot(themolst) |
| 49 | 56 |
|
| 50 | 57 |
########### |
| ... | ... | |
| 154 | 161 |
#fitbias<-Krig(bias_xy,sta_bias,theta=1e5) #use TPS or krige |
| 155 | 162 |
fitclim<-Krig(clim_xy,sta_clim,theta=1e5) |
| 156 | 163 |
|
| 157 |
#The output is a krig object using fields |
|
| 164 |
#The output is a krig object using fields: modif 10/30
|
|
| 158 | 165 |
#mod9a<-fitbias |
| 159 |
mod9a<-fitclim |
|
| 166 |
mod_krtmp1<-fitclim |
|
| 167 |
model_name<-paste("mod_kr","month",sep="_")
|
|
| 168 |
assign(model_name,mod_krtmp1) |
|
| 160 | 169 |
|
| 161 | 170 |
# Creating plot of bias surface and saving it |
| 162 | 171 |
#X11() |
| ... | ... | |
| 184 | 193 |
fitdelta<-Krig(daily_sta_xy,daily_delta,theta=1e5) #use TPS or krige |
| 185 | 194 |
fitdeltaclim<-Krig(daily_sta_xy,daily_deltaclim,theta=1e5) #use TPS or krige |
| 186 | 195 |
|
| 187 |
#Kriging using fields package |
|
| 196 |
#Kriging using fields package: modif 10/30
|
|
| 188 | 197 |
#mod9b<-fitdelta |
| 189 |
mod9b<-fitdeltaclim |
|
| 198 |
mod_krtmp2<-fitdeltaclim |
|
| 199 |
model_name<-paste("mod_kr","day",sep="_")
|
|
| 200 |
assign(model_name,mod_krtmp2) |
|
| 201 |
|
|
| 190 | 202 |
# Creating plot of bias surface and saving it |
| 191 | 203 |
#X11() |
| 192 | 204 |
png(paste("Deltaclim_surface_TMax_",sampling_dat$date[i],"_",sampling_dat$prop[i],"_",sampling_dat$run_samp[i],
|
| ... | ... | |
| 310 | 322 |
|
| 311 | 323 |
#Model and response variable can be changed without affecting the script |
| 312 | 324 |
|
| 313 |
formula1 <- as.formula("y_var ~ s(lat) + s(lon) + s(ELEV_SRTM)", env=.GlobalEnv)
|
|
| 314 |
formula2 <- as.formula("y_var~ s(lat,lon)+ s(ELEV_SRTM)", env=.GlobalEnv)
|
|
| 315 |
formula3 <- as.formula("y_var~ s(lat) + s (lon) + s (ELEV_SRTM) + s (Northness)+ s (Eastness) + s(DISTOC)", env=.GlobalEnv)
|
|
| 316 |
formula4 <- as.formula("y_var~ s(lat) + s (lon) + s(ELEV_SRTM) + s(Northness) + s (Eastness) + s(DISTOC) + s(LST)", env=.GlobalEnv)
|
|
| 317 |
formula5 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)", env=.GlobalEnv)
|
|
| 318 |
formula6 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC1)", env=.GlobalEnv)
|
|
| 319 |
formula7 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC3)", env=.GlobalEnv)
|
|
| 320 |
formula8 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST) + s(LC1,LC3)", env=.GlobalEnv)
|
|
| 325 |
#formula1 <- as.formula("y_var ~ s(lat) + s(lon) + s(ELEV_SRTM)", env=.GlobalEnv)
|
|
| 326 |
#formula2 <- as.formula("y_var~ s(lat,lon)+ s(ELEV_SRTM)", env=.GlobalEnv)
|
|
| 327 |
#formula3 <- as.formula("y_var~ s(lat) + s (lon) + s (ELEV_SRTM) + s (Northness)+ s (Eastness) + s(DISTOC)", env=.GlobalEnv)
|
|
| 328 |
#formula4 <- as.formula("y_var~ s(lat) + s (lon) + s(ELEV_SRTM) + s(Northness) + s (Eastness) + s(DISTOC) + s(LST)", env=.GlobalEnv)
|
|
| 329 |
#formula5 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)", env=.GlobalEnv)
|
|
| 330 |
#formula6 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC1)", env=.GlobalEnv)
|
|
| 331 |
#formula7 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC3)", env=.GlobalEnv)
|
|
| 332 |
#formula8 <- as.formula("y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST) + s(LC1,LC3)", env=.GlobalEnv)
|
|
| 333 |
|
|
| 334 |
list_formulas<-vector("list",nmodels)
|
|
| 335 |
|
|
| 336 |
#This can be entered as textfile or option later...ok for running now on 10/30/2012 |
|
| 337 |
list_formulas[[1]] <- as.formula("y_var~ s(ELEV_SRTM)", env=.GlobalEnv)
|
|
| 338 |
list_formulas[[2]] <- as.formula("y_var~ s(LST)", env=.GlobalEnv)
|
|
| 339 |
list_formulas[[3]] <- as.formula("y_var~ s(LST) + s(ELEV_SRTM)", env=.GlobalEnv)
|
|
| 340 |
list_formulas[[4]] <- as.formula("y_var~ s(LST,ELEV_SRTM)", env=.GlobalEnv)
|
|
| 341 |
list_formulas[[5]] <- as.formula("y_var~ s(lat,lon,ELEV_SRTM)", env=.GlobalEnv)
|
|
| 342 |
|
|
| 321 | 343 |
|
| 322 | 344 |
#mod1<- try(gam(formula1, data=data_s)) |
| 323 | 345 |
#mod2<- try(gam(formula2, data=data_s)) #modified nesting....from 3 to 2 |
| ... | ... | |
| 329 | 351 |
#mod8<- try(gam(formula8, data=data_s)) |
| 330 | 352 |
|
| 331 | 353 |
if (climgam==1){ #This will automatically use monthly station data in the second step
|
| 332 |
mod1<- try(gam(formula1, data=data_month)) |
|
| 333 |
mod2<- try(gam(formula2, data=data_month)) #modified nesting....from 3 to 2 |
|
| 334 |
mod3<- try(gam(formula3, data=data_month)) |
|
| 335 |
mod4<- try(gam(formula4, data=data_month)) |
|
| 336 |
mod5<- try(gam(formula5, data=data_month)) |
|
| 337 |
mod6<- try(gam(formula6, data=data_month)) |
|
| 338 |
mod7<- try(gam(formula7, data=data_month)) |
|
| 339 |
mod8<- try(gam(formula8, data=data_month)) |
|
| 354 |
#mod1<- try(gam(formula1, data=data_month)) |
|
| 355 |
#mod2<- try(gam(formula2, data=data_month)) #modified nesting....from 3 to 2 |
|
| 356 |
#mod3<- try(gam(formula3, data=data_month)) |
|
| 357 |
#mod4<- try(gam(formula4, data=data_month)) |
|
| 358 |
#mod5<- try(gam(formula5, data=data_month)) |
|
| 359 |
#mod6<- try(gam(formula6, data=data_month)) Change this in a loop around model !!! mod-j .... |
|
| 360 |
#mod7<- try(gam(formula7, data=data_month)) |
|
| 361 |
#mod8<- try(gam(formula8, data=data_month)) |
|
| 362 |
|
|
| 363 |
for (j in 1:nmodels){
|
|
| 364 |
formula<-list_formulas[[j]] |
|
| 365 |
mod<- try(gam(formula, data=data_month)) |
|
| 366 |
model_name<-paste("mod",j,sep="")
|
|
| 367 |
assign(model_name,mod) |
|
| 368 |
} |
|
| 340 | 369 |
|
| 341 | 370 |
} else if (climgam==0){ #This will use daily delta in the second step
|
| 342 | 371 |
|
| 343 |
mod1<- try(gam(formula1, data=data_s)) |
|
| 344 |
mod2<- try(gam(formula2, data=data_s)) #modified nesting....from 3 to 2 |
|
| 345 |
mod3<- try(gam(formula3, data=data_s)) |
|
| 346 |
mod4<- try(gam(formula4, data=data_s)) |
|
| 347 |
mod5<- try(gam(formula5, data=data_s)) |
|
| 348 |
mod6<- try(gam(formula6, data=data_s)) |
|
| 349 |
mod7<- try(gam(formula7, data=data_s)) |
|
| 350 |
mod8<- try(gam(formula8, data=data_s)) |
|
| 372 |
#mod1<- try(gam(formula1, data=data_s)) |
|
| 373 |
#mod2<- try(gam(formula2, data=data_s)) #modified nesting....from 3 to 2 |
|
| 374 |
#mod3<- try(gam(formula3, data=data_s)) |
|
| 375 |
#mod4<- try(gam(formula4, data=data_s)) |
|
| 376 |
#mod5<- try(gam(formula5, data=data_s)) |
|
| 377 |
#mod6<- try(gam(formula6, data=data_s)) |
|
| 378 |
#mod7<- try(gam(formula7, data=data_s)) |
|
| 379 |
#mod8<- try(gam(formula8, data=data_s)) |
|
| 380 |
for (j in 1:nmodels){
|
|
| 381 |
formula<-list_formulas[[j]] |
|
| 382 |
mod<- try(gam(formula, data=data_s)) |
|
| 383 |
model_name<-paste("mod",j,sep="")
|
|
| 384 |
assign(model_name,mod) |
|
| 385 |
} |
|
| 386 |
|
|
| 351 | 387 |
} |
| 352 | 388 |
|
| 353 | 389 |
### Added by benoit |
| ... | ... | |
| 402 | 438 |
#nv<-nrow(data_v) |
| 403 | 439 |
#browser() |
| 404 | 440 |
|
| 441 |
mod_obj<-vector("list",nmodels+2) #This will contain the model objects fitting: 10/30
|
|
| 442 |
mod_obj[[nmodels+1]]<-mod_kr_month #Storing climatology object |
|
| 443 |
mod_obj[[nmodels+2]]<-mod_kr_day #Storing delta object |
|
| 444 |
|
|
| 405 | 445 |
for (j in 1:nmodels){
|
| 406 | 446 |
|
| 407 | 447 |
##Model assessment: specific diagnostic/metrics for GAM |
| 408 | 448 |
|
| 409 | 449 |
name<-paste("mod",j,sep="") #modj is the name of The "j" model (mod1 if j=1)
|
| 410 | 450 |
mod<-get(name) #accessing GAM model ojbect "j" |
| 411 |
|
|
| 451 |
mod_obj[[j]]<-mod #storing current model object |
|
| 452 |
|
|
| 412 | 453 |
#If mod "j" is not a model object |
| 413 | 454 |
if (inherits(mod,"try-error")) {
|
| 414 | 455 |
results_m1[1,1]<- sampling_dat$date[i] #storing the interpolation dates in the first column |
| ... | ... | |
| 627 | 668 |
tb_metrics1<-rbind(results_table_RMSE,results_table_MAE, results_table_ME, |
| 628 | 669 |
results_table_R2,results_table_RMSE_f,results_table_MAE_f,results_table_R2_f) # |
| 629 | 670 |
tb_metrics2<-rbind(results_table_m1,results_table_m2, results_table_m3) |
| 630 |
cname<-c("dates","ns","metric","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7","mod8","mod9")
|
|
| 631 |
colnames(tb_metrics1)<-cname |
|
| 632 |
cname<-c("dates","ns","metric","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7","mod8")
|
|
| 633 |
colnames(tb_metrics2)<-cname |
|
| 634 |
#colnames(results_table_RMSE)<-cname |
|
| 635 |
#colnames(results_table_RMSE_f)<-cname |
|
| 636 |
#tb_diagnostic1<-results_table_RMSE #measures of validation |
|
| 637 |
#tb_diagnostic2<-results_table_RMSE_f #measures of fit |
|
| 638 |
|
|
| 639 |
#write.table(tb_diagnostic1, file= paste(path,"/","results_fusion_Assessment_measure1",out_prefix,".txt",sep=""), sep=",") |
|
| 640 | 671 |
|
| 641 |
#} |
|
| 672 |
#Preparing labels |
|
| 673 |
mod_labels<-rep("mod",nmodels+1)
|
|
| 674 |
index<-as.character(1:(nmodels+1)) |
|
| 675 |
mod_labels<-paste(mod_labels,index,sep="") |
|
| 676 |
cname<-c("dates","ns","metric", mod_labels)
|
|
| 677 |
#cname<-c("dates","ns","metric","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7","mod8","mod9")
|
|
| 678 |
colnames(tb_metrics1)<-cname |
|
| 679 |
#cname<-c("dates","ns","metric","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7","mod8")
|
|
| 680 |
colnames(tb_metrics2)<-cname[1:(nmodels+3)] |
|
| 681 |
|
|
| 642 | 682 |
print(paste(sampling_dat$date[i],"processed")) |
| 643 | 683 |
# Kriging object may need to be modified...because it contains the full image of prediction!! |
| 644 | 684 |
##loop through model objects data frame and set field to zero... |
| 645 | 685 |
|
| 646 |
mod_obj<-list(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,mod9a,mod9b) |
|
| 647 |
names(mod_obj)<-c("mod1","mod2","mod3","mod4","mod5","mod6","mod7","mod8","mod9a","mod9b") #generate names automatically??
|
|
| 648 |
#results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2)
|
|
| 649 |
#results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2,mod_obj)
|
|
| 650 |
results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2,mod_obj,data_month) |
|
| 651 |
names(results_list)<-c("data_s","data_v","tb_metrics1","tb_metrics2","mod_obj","data_month")
|
|
| 686 |
#mod_obj<-list(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,mod9a,mod9b)
|
|
| 687 |
#names(mod_obj)<-c("mod1","mod2","mod3","mod4","mod5","mod6","mod7","mod8","mod9a","mod9b") #generate names automatically??
|
|
| 688 |
mod_labels_kr<-c("mod_kr_month", "mod_kr_day")
|
|
| 689 |
names(mod_obj)<-c(mod_labels[1:nmodels],mod_labels_kr)
|
|
| 690 |
results_list<-list(data_s,data_v,tb_metrics1,tb_metrics2,mod_obj,data_month,list_formulas)
|
|
| 691 |
names(results_list)<-c("data_s","data_v","tb_metrics1","tb_metrics2","mod_obj","data_month","formulas")
|
|
| 652 | 692 |
save(results_list,file= paste(path,"/","results_list_metrics_objects_",sampling_dat$date[i],"_",sampling_dat$prop[i],"_",sampling_dat$run_samp[i], |
| 653 | 693 |
out_prefix,".RData",sep="")) |
| 654 | 694 |
return(results_list) |
Also available in: Unified diff
GAM CAI function correction constant sampling and GAM climatology models, tmax OR