Revision 65d96c1a
Added by Benoit Parmentier over 11 years ago
| climate/research/oregon/interpolation/GAM_fusion_analysis_raster_prediction_multisampling.R | ||
|---|---|---|
| 1 |
######################### Raster prediction GAM FUSION ####################################
|
|
| 1 |
######################### Raster prediction #################################### |
|
| 2 | 2 |
############################ Interpolation of temperature for given processing region ########################################## |
| 3 | 3 |
#This script interpolates temperature values using MODIS LST, covariates and GHCND station data. |
| 4 | 4 |
#It requires the text file of stations and a shape file of the study area. |
| ... | ... | |
| 8 | 8 |
#2)Constant sampling: use the same sample over the runs |
| 9 | 9 |
#3)over dates: run over for example 365 dates without mulitsampling |
| 10 | 10 |
#4)use seed number: use seed if random samples must be repeatable |
| 11 |
#5)GAM fusion: possibilty of running GAM+FUSION or GAM+CAI and other options added
|
|
| 11 |
#5)possibilty of running GAM+FUSION or GAM+CAI and other options added |
|
| 12 | 12 |
#The interpolation is done first at the monthly time scale then delta surfaces are added. |
| 13 | 13 |
#AUTHOR: Benoit Parmentier |
| 14 |
#DATE: 04/05/2013
|
|
| 14 |
#DATE: 04/22/2013
|
|
| 15 | 15 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#568-- |
| 16 | 16 |
# |
| 17 | 17 |
# TO DO: |
| ... | ... | |
| 22 | 22 |
|
| 23 | 23 |
################################################################################################### |
| 24 | 24 |
|
| 25 |
raster_prediction_gam_fusion<-function(list_param_raster_prediction){
|
|
| 25 |
raster_prediction_fun <-function(list_param_raster_prediction){
|
|
| 26 | 26 |
|
| 27 | 27 |
##Function to predict temperature interpolation with 21 input parameters |
| 28 | 28 |
#9 parameters used in the data preparation stage and input in the current script |
| ... | ... | |
| 191 | 191 |
t1<-proc.time() |
| 192 | 192 |
j=12 |
| 193 | 193 |
#browser() |
| 194 |
list_param_runClim_KGFusion<-list(j,s_raster,covar_names,lst_avg,list_models,dst,var,y_var_name, out_prefix) |
|
| 195 |
names(list_param_runClim_KGFusion)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","var","y_var_name","out_prefix")
|
|
| 196 |
#source(file.path(script_path,"GAM_fusion_function_multisampling_03122013.R")) |
|
| 197 |
gamclim_fus_mod<-mclapply(1:12, list_param=list_param_runClim_KGFusion, runClim_KGFusion,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement |
|
| 198 |
#test<-runClim_KGFusion(1,list_param=list_param_runClim_KGFusion) |
|
| 199 |
#gamclim_fus_mod<-mclapply(1:6, list_param=list_param_runClim_KGFusion, runClim_KGFusion,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement |
|
| 194 |
if (interpolation_method=="gam_fusion"){
|
|
| 195 |
list_param_runClim_KGFusion<-list(j,s_raster,covar_names,lst_avg,list_models,dst,var,y_var_name, out_prefix) |
|
| 196 |
names(list_param_runClim_KGFusion)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","var","y_var_name","out_prefix")
|
|
| 197 |
#source(file.path(script_path,"GAM_fusion_function_multisampling_03122013.R")) |
|
| 198 |
clim_method_mod_obj<-mclapply(1:12, list_param=list_param_runClim_KGFusion, runClim_KGFusion,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement |
|
| 199 |
#test<-runClim_KGFusion(1,list_param=list_param_runClim_KGFusion) |
|
| 200 |
#gamclim_fus_mod<-mclapply(1:6, list_param=list_param_runClim_KGFusion, runClim_KGFusion,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement |
|
| 201 |
|
|
| 202 |
} |
|
| 200 | 203 |
|
| 204 |
if (interpolation_method=="gam_CAI"){
|
|
| 205 |
list_param_runClim_KGCAI<-list(j,s_raster,covar_names,lst_avg,list_models,dst,var,y_var_name, out_prefix) |
|
| 206 |
names(list_param_runClim_KGCAI)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","var","y_var_name","out_prefix")
|
|
| 207 |
clim_method_mod_obj<-mclapply(1:12, list_param=list_param_runClim_KGCAI, runClim_KGCAI,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement |
|
| 208 |
#test<-runClim_KGCAI(1,list_param=list_param_runClim_KGCAI) |
|
| 209 |
#gamclim_fus_mod<-mclapply(1:6, list_param=list_param_runClim_KGFusion, runClim_KGFusion,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement |
|
| 210 |
} |
|
| 211 |
|
|
| 201 | 212 |
#gamclim_fus_mod<-runClim_KGFusion(1,list_param=list_param_runClim_KGFusion) #This is the end bracket from mclapply(...) statement |
| 202 |
save(gamclim_fus_mod,file= paste("gamclim_fus_mod_",y_var_name,out_prefix,".RData",sep=""))
|
|
| 213 |
save(clim_method_mod_obj,file= paste(interpolation_method,"_mod_",y_var_name,out_prefix,".RData",sep=""))
|
|
| 203 | 214 |
t2<-proc.time()-t1 |
| 204 | 215 |
writeLines(as.character(t2),con=log_file,sep="\n") |
| 205 | 216 |
|
| 206 | 217 |
#now get list of raster clim layers |
| 207 | 218 |
|
| 208 |
list_tmp<-vector("list",length(gamclim_fus_mod))
|
|
| 209 |
for (i in 1:length(gamclim_fus_mod)){
|
|
| 210 |
tmp<-gamclim_fus_mod[[i]]$clim
|
|
| 219 |
list_tmp<-vector("list",length(clim_method_mod_obj))
|
|
| 220 |
for (i in 1:length(clim_method_mod_obj)){
|
|
| 221 |
tmp<-clim_method_mod_obj[[i]]$clim
|
|
| 211 | 222 |
list_tmp[[i]]<-tmp |
| 212 | 223 |
} |
| 213 | 224 |
|
| 214 | 225 |
################## PREDICT AT DAILY TIME SCALE ################# |
| 226 |
#Predict at daily time scale from single time scale or multiple time scale methods: 2 methods availabe now |
|
| 215 | 227 |
|
| 216 | 228 |
#put together list of clim models per month... |
| 217 | 229 |
#rast_clim_yearlist<-list_tmp |
| 218 | 230 |
clim_yearlist<-list_tmp |
| 219 | 231 |
#Second predict at the daily time scale: delta |
| 220 | 232 |
|
| 221 |
#gam_fus_mod<-mclapply(1:1, runGAMFusion,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement
|
|
| 233 |
#method_mod_obj<-mclapply(1:1, runGAMFusion,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement
|
|
| 222 | 234 |
writeLines("Predictions at the daily scale:",con=log_file,sep="\n")
|
| 223 | 235 |
t1<-proc.time() |
| 224 | 236 |
|
| 225 |
#input a list:note that ghcn.subsets is not sampling_obj$data_day_ghcn |
|
| 226 |
list_param_runGAMFusion<-list(i,clim_yearlist,sampling_obj,dst,var,y_var_name, out_prefix) |
|
| 227 |
names(list_param_runGAMFusion)<-c("list_index","clim_yearlist","sampling_obj","dst","var","y_var_name","out_prefix")
|
|
| 228 |
#test<-mclapply(1:18, runGAMFusion,list_param=list_param_runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) |
|
| 229 |
#test<-runGAMFusion(1,list_param=list_param_runGAMFusion) |
|
| 230 |
|
|
| 231 |
#MAKE IT GENERAL: for any method: replace "gam_fus_mod" by "method_mod_obj"? |
|
| 232 |
|
|
| 233 |
gam_fus_mod<-mclapply(1:length(sampling_obj$ghcn_data_day),list_param=list_param_runGAMFusion,runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement |
|
| 234 |
|
|
| 235 |
#gam_fus_mod<-mclapply(1:length(sampling_obj$ghcn_data_day),runGAMFusion,list_param_runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement |
|
| 236 |
#gam_fus_mod<-mclapply(1:length(ghcn.subsets), runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement |
|
| 237 |
|
|
| 238 |
save(gam_fus_mod,file= paste("gam_fus_mod_",y_var_name,out_prefix,".RData",sep=""))
|
|
| 237 |
if (interpolation_method=="gam_CAI" | interpolation_method=="gam_fusion"){
|
|
| 238 |
#input a list:note that ghcn.subsets is not sampling_obj$data_day_ghcn |
|
| 239 |
list_param_run_prediction_daily_deviation <-list(i,clim_yearlist,sampling_obj,dst,var,y_var_name, interpolation_method,out_prefix) |
|
| 240 |
names(list_param_run_prediction_daily_deviation)<-c("list_index","clim_yearlist","sampling_obj","dst","var","y_var_name","interpolation_method","out_prefix")
|
|
| 241 |
#test<-mclapply(1:18, runGAMFusion,list_param=list_param_runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) |
|
| 242 |
#test<-runGAMFusion(1,list_param=list_param_runGAMFusion) |
|
| 243 |
|
|
| 244 |
method_mod_obj<-mclapply(1:length(sampling_obj$ghcn_data_day),list_param=list_param_run_prediction_daily_deviation,run_prediction_daily_deviation,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement |
|
| 245 |
#method_mod_obj<-mclapply(1:1,list_param=list_param_run_prediction_daily_deviation,run_prediction_daily_deviation,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement |
|
| 246 |
|
|
| 247 |
#method_mod_obj<-mclapply(1:length(sampling_obj$ghcn_data_day),runGAMFusion,list_param_runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement |
|
| 248 |
#method_mod_obj<-mclapply(1:length(ghcn.subsets), runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement |
|
| 249 |
|
|
| 250 |
} |
|
| 251 |
|
|
| 252 |
save(method_mod_obj,file= paste("method_mod_obj_",interpolation_method,"_",y_var_name,out_prefix,".RData",sep=""))
|
|
| 239 | 253 |
t2<-proc.time()-t1 |
| 240 | 254 |
writeLines(as.character(t2),con=log_file,sep="\n") |
| 241 | 255 |
#browser() |
| ... | ... | |
| 243 | 257 |
############### NOW RUN VALIDATION ######################### |
| 244 | 258 |
#SIMPLIFY THIS PART: one call |
| 245 | 259 |
|
| 246 |
list_tmp<-vector("list",length(gam_fus_mod))
|
|
| 247 |
for (i in 1:length(gam_fus_mod)){
|
|
| 248 |
tmp<-gam_fus_mod[[i]][[y_var_name]] #y_var_name is the variable predicted (dailyTmax or dailyTmin)
|
|
| 260 |
list_tmp<-vector("list",length(method_mod_obj))
|
|
| 261 |
for (i in 1:length(method_mod_obj)){
|
|
| 262 |
tmp<-method_mod_obj[[i]][[y_var_name]] #y_var_name is the variable predicted (dailyTmax or dailyTmin)
|
|
| 249 | 263 |
list_tmp[[i]]<-tmp |
| 250 | 264 |
} |
| 251 | 265 |
rast_day_yearlist<-list_tmp #list of predicted images |
| ... | ... | |
| 253 | 267 |
writeLines("Validation step:",con=log_file,sep="\n")
|
| 254 | 268 |
t1<-proc.time() |
| 255 | 269 |
#calculate_accuary_metrics<-function(i) |
| 256 |
list_param_validation<-list(i,rast_day_yearlist,gam_fus_mod,y_var_name, out_prefix)
|
|
| 270 |
list_param_validation<-list(i,rast_day_yearlist,method_mod_obj,y_var_name, out_prefix)
|
|
| 257 | 271 |
names(list_param_validation)<-c("list_index","rast_day_year_list","method_mod_obj","y_var_name","out_prefix") #same names for any method
|
| 258 | 272 |
|
| 259 |
#gam_fus_validation_mod<-mclapply(1:length(gam_fus_mod), calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement
|
|
| 260 |
gam_fus_validation_mod<-mclapply(1:length(gam_fus_mod), list_param=list_param_validation, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement
|
|
| 273 |
#gam_fus_validation_mod<-mclapply(1:length(method_mod_obj), calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement
|
|
| 274 |
validation_mod_obj <-mclapply(1:length(method_mod_obj), list_param=list_param_validation, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement
|
|
| 261 | 275 |
|
| 262 | 276 |
#gam_fus_validation_mod<-mclapply(1:1, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement |
| 263 |
save(gam_fus_validation_mod,file= paste("gam_fus_validation_mod_",y_var_name,out_prefix,".RData",sep=""))
|
|
| 277 |
save(validation_mod_obj,file= paste(interpolation_method,"_validation_mod_obj_",y_var_name,out_prefix,".RData",sep=""))
|
|
| 264 | 278 |
t2<-proc.time()-t1 |
| 265 | 279 |
writeLines(as.character(t2),con=log_file,sep="\n") |
| 266 | 280 |
|
| 267 | 281 |
#################### ASSESSMENT OF PREDICTIONS: PLOTS OF ACCURACY METRICS ########### |
| 268 | 282 |
|
| 269 | 283 |
##Create data.frame with valiation metrics for a full year |
| 270 |
tb_diagnostic_v<-extract_from_list_obj(gam_fus_validation_mod,"metrics_v")
|
|
| 284 |
tb_diagnostic_v<-extract_from_list_obj(validation_mod_obj,"metrics_v")
|
|
| 271 | 285 |
rownames(tb_diagnostic_v)<-NULL #remove row names |
| 272 | 286 |
|
| 273 | 287 |
#Call function to create plots of metrics for validation dataset |
| ... | ... | |
| 290 | 304 |
################### PREPARE RETURN OBJECT ############### |
| 291 | 305 |
#Will add more information to be returned |
| 292 | 306 |
|
| 293 |
raster_prediction_obj<-list(gamclim_fus_mod,gam_fus_mod,gam_fus_validation_mod,tb_diagnostic_v,summary_metrics_v)
|
|
| 294 |
names(raster_prediction_obj)<-c("gamclim_fus_mod","gam_fus_mod","gam_fus_validation_mod","tb_diagnostic_v",
|
|
| 307 |
raster_prediction_obj<-list(clim_method_mod_obj,method_mod_obj,validation_mod_obj,tb_diagnostic_v,summary_metrics_v)
|
|
| 308 |
names(raster_prediction_obj)<-c("clim_method_mod_obj","method_mod_obj","validation_mod_obj","tb_diagnostic_v",
|
|
| 295 | 309 |
"summary_metrics_v") |
| 296 |
save(raster_prediction_obj,file= paste("raster_prediction_obj_",y_var_name,out_prefix,".RData",sep=""))
|
|
| 310 |
save(raster_prediction_obj,file= paste("raster_prediction_obj_",interpolation_method,"_", y_var_name,out_prefix,".RData",sep=""))
|
|
| 297 | 311 |
|
| 298 | 312 |
return(raster_prediction_obj) |
| 299 | 313 |
} |
Also available in: Unified diff
raster prediction function, modifications to make code general for any interpolation method