Revision da42ba27
Added by Benoit Parmentier over 11 years ago
| climate/research/oregon/interpolation/GAM_CAI_analysis_raster_prediction_multisampling.R | ||
|---|---|---|
| 6 | 6 |
#Method is assedsed using constant sampling with variation of validation sample with different # |
| 7 | 7 |
#hold out proportions. # |
| 8 | 8 |
#AUTHOR: Benoit Parmentier # |
| 9 |
#DATE: 10/30/2012 #
|
|
| 9 |
#DATE: 12/27/2012 #
|
|
| 10 | 10 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#491-- # |
| 11 | 11 |
################################################################################################### |
| 12 | 12 |
|
| ... | ... | |
| 43 | 43 |
#GHCN Database for 1980-2010 for study area (OR) |
| 44 | 44 |
data3<-read.table(paste(path,"/","ghcn_data_TMAXy1980_2010_OR_0602012.txt",sep=""),sep=",", header=TRUE) |
| 45 | 45 |
|
| 46 |
nmodels<-5 #number of models running
|
|
| 47 |
y_var_name<-"dailyTmax" |
|
| 48 |
climgam=1 #if 1, then GAM is run on the climatology rather than the daily deviation surface...
|
|
| 49 |
predval<-1 |
|
| 50 |
prop<-0.3 #Proportion of testing retained for validation
|
|
| 46 |
nmodels<-9 #number of models running
|
|
| 47 |
y_var_name<-"dailyTmax" #climate variable interpolated
|
|
| 48 |
climgam<-1 #if 1, then GAM is run on the climatology rather than the daily deviation surface...
|
|
| 49 |
predval<-1 #if 1, produce raster prediction
|
|
| 50 |
prop<-0.3 #Proportion of testing retained for validation |
|
| 51 | 51 |
|
| 52 | 52 |
seed_number<- 100 #Seed number for random sampling, if seed_number<0, no seed number is used.. |
| 53 |
#out_prefix<-"_365d_GAM_CAI2_const_10222012_" #User defined output prefix
|
|
| 54 |
#out_prefix<-"_365d_GAM_CAI2_const_all_lstd_10272012" #User defined output prefix
|
|
| 55 |
out_prefix<-"_365d_GAM_CAI3_all_10302012" #User defined output prefix
|
|
| 53 |
#out_prefix<-"_365d_GAM_CAI2_const_10222012_" #User defined output prefix |
|
| 54 |
#out_prefix<-"_365d_GAM_CAI2_const_all_lstd_10272012" #User defined output prefix |
|
| 55 |
out_prefix<-"_365d_GAM_CAI4_all_12272012" #User defined output prefix
|
|
| 56 | 56 |
|
| 57 | 57 |
bias_val<-0 #if value 1 then daily training data is used in the bias surface rather than the all monthly stations (added on 07/11/2012) |
| 58 | 58 |
bias_prediction<-1 #if value 1 then use GAM for the BIAS prediction otherwise GAM direct reprediction for y_var (daily tmax) |
| ... | ... | |
| 60 | 60 |
prop_min<-0.3 #if prop_min=prop_max and step=0 then predicitons are done for the number of dates... |
| 61 | 61 |
prop_max<-0.3 |
| 62 | 62 |
step<-0 |
| 63 |
constant<-0 #if value 1 then use the same samples as date one for the all set of dates
|
|
| 63 |
constant<-0 #if value 1 then use the same sample used in the first date for interpolation over the set of dates
|
|
| 64 | 64 |
#projection used in the interpolation of the study area |
| 65 | 65 |
CRS_interp<-"+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"; |
| 66 | 66 |
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| 67 | 67 |
|
| 68 |
#This can be entered as textfile or option later...ok for running now on 12/07/2012 |
|
| 69 |
list_formulas<-vector("list",nmodels)
|
|
| 70 |
|
|
| 71 |
list_formulas[[1]] <- as.formula("y_var~ s(ELEV_SRTM)", env=.GlobalEnv)
|
|
| 72 |
list_formulas[[2]] <- as.formula("y_var~ s(LST)", env=.GlobalEnv)
|
|
| 73 |
list_formulas[[3]] <- as.formula("y_var~ s(ELEV_SRTM,LST)", env=.GlobalEnv)
|
|
| 74 |
list_formulas[[4]] <- as.formula("y_var~ s(lat)+s(lon)+s(ELEV_SRTM)", env=.GlobalEnv)
|
|
| 75 |
list_formulas[[5]] <- as.formula("y_var~ s(lat,lon,ELEV_SRTM)", env=.GlobalEnv)
|
|
| 76 |
list_formulas[[6]] <- as.formula("y_var~ s(lat,lon)+s(ELEV_SRTM)+s(Northness_w,Eastness_w)+s(LST)", env=.GlobalEnv)
|
|
| 77 |
list_formulas[[7]] <- as.formula("y_var~ s(lat,lon)+s(ELEV_SRTM)+s(Northness_w,Eastness_w)+s(LST)+s(LC1)", env=.GlobalEnv)
|
|
| 78 |
list_formulas[[8]] <- as.formula("y_var~ s(lat,lon)+s(ELEV_SRTM)+s(Northness_w,Eastness_w)+s(LST)+s(LC3)", env=.GlobalEnv)
|
|
| 79 |
list_formulas[[9]] <- as.formula("y_var~ s(x_OR83M,y_OR83M)", env=.GlobalEnv)
|
|
| 80 |
|
|
| 68 | 81 |
#source("GAM_CAI_function_multisampling_10252012.R")
|
| 69 |
source("GAM_CAI_function_multisampling_10302012.R")
|
|
| 82 |
source("GAM_CAI_function_multisampling_12072012.R")
|
|
| 70 | 83 |
|
| 71 | 84 |
############ START OF THE SCRIPT ################## |
| 72 | 85 |
|
| ... | ... | |
| 315 | 328 |
#Start loop here... |
| 316 | 329 |
|
| 317 | 330 |
#gam_CAI_mod<-mclapply(1:length(dates), runGAMCAI,mc.preschedule=FALSE,mc.cores = 8) #This is the end bracket from mclapply(...) statement |
| 318 |
gam_CAI_mod<-mclapply(1:length(ghcn.subsets), runGAMCAI,mc.preschedule=FALSE,mc.cores = 8) #This is the end bracket from mclapply(...) statement |
|
| 319 |
#gam_CAI_mod<-mclapply(1:1, runGAMCAI,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement |
|
| 331 |
gam_CAI_mod<-mclapply(1:length(ghcn.subsets), runGAMCAI,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement |
|
| 332 |
#gam_CAI_mod<-mclapply(1:2, runGAMCAI,mc.preschedule=FALSE,mc.cores = 2) #This is the end bracket from mclapply(...) statement# |
|
| 333 |
#gam_CAI_mod<-mclapply(1:2, runGAMCAI,mc.preschedule=FALSE,mc.cores = 2) #This is the end bracket from mclapply(...) statement |
|
| 320 | 334 |
|
| 321 | 335 |
tb<-gam_CAI_mod[[1]][[3]][0,] #empty data frame with metric table structure that can be used in rbinding... |
| 322 | 336 |
tb_tmp<-gam_CAI_mod #copy |
| ... | ... | |
| 341 | 355 |
assign(metric_name,tb_metric) |
| 342 | 356 |
tb_metric_list[[i]]<-tb_metric |
| 343 | 357 |
} |
| 358 |
mod_labels<-rep("mod",nmodels+1)
|
|
| 359 |
index<-as.character(1:(nmodels+1)) |
|
| 360 |
mod_labels<-paste(mod_labels,index,sep="") |
|
| 344 | 361 |
|
| 345 | 362 |
tb_diagnostic<-do.call(rbind,tb_metric_list) #produce a data.frame from the list ... |
| 346 | 363 |
tb_diagnostic[["prop"]]<-as.factor(tb_diagnostic[["prop"]]) |
| 347 | 364 |
|
| 348 | 365 |
t<-melt(tb_diagnostic, |
| 349 |
measure=c("mod1","mod2","mod3","mod4", "mod5", "mod6", "mod7", "mod8","mod9"),
|
|
| 366 |
measure=mod_labels,
|
|
| 350 | 367 |
id=c("dates","metric","prop"),
|
| 351 | 368 |
na.rm=F) |
| 352 | 369 |
avg_tb<-cast(t,metric+prop~variable,mean) |
| ... | ... | |
| 367 | 384 |
save(sampling_obj, file= paste(path,"/","results2_CAI_sampling_obj",out_prefix,".RData",sep="")) |
| 368 | 385 |
save(gam_CAI_mod,file= paste(path,"/","results2_CAI_Assessment_measure_all",out_prefix,".RData",sep="")) |
| 369 | 386 |
|
| 387 |
#new combined object used since november 2012 |
|
| 388 |
gam_CAI_mod_obj<-list(gam_CAI_mod=gam_CAI_mod,sampling_obj=sampling_obj) |
|
| 389 |
save(gam_CAI_mod_obj,file= paste(path,"/","results_mod_obj_",out_prefix,".RData",sep="")) |
|
| 390 |
|
|
| 370 | 391 |
#### END OF SCRIPT |
Also available in: Unified diff
GAM CAI raster predictions, modifcations and model running for IBS conference 2013