Project

General

Profile

« Previous | Next » 

Revision fc7c30c1

Added by Benoit Parmentier over 11 years ago

raster prediction function, modifications related to kriging daily method

View differences:

climate/research/oregon/interpolation/GAM_fusion_analysis_raster_prediction_multisampling.R
11 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: 05/06/2013
14 
#DATE: 06/05/2013
15 15 
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#568--
16 16 
#
17 17 
# TO DO:
......
68 68 
  library(plotrix)
69 69 
  library(maptools)
70 70 
  library(gdata) #Nesssary to use cbindX
71
71 
  library(automap)
72 72 
  ### Parameters and arguments
73 73 
  #PARSING INPUTS/ARGUMENTS
74 74 
#
......
126 126 
  #log_fname<-paste("R_log_raster_prediction",out_prefix, ".log",sep="")
127 127 
  log_fname<-paste("R_log_raster_prediction",out_prefix, ".log",sep="")
128 128 
  #sink(log_fname) #create new log file
129 
  file.create(file.path(path,log_fname)) #create new log file
129 
  file.create(file.path(out_path,log_fname)) #create new log file
130 130
131 131 
  time1<-proc.time()    #Start stop watch
132 132
......
148 148 
    dates<-dates_selected #dates to be predicted
149 149 
  }
150 150
151 
  #Reaccccding of covariate brick covariates can be changed...
151 
  #Reading in covariate brickcan be changed...
152 152
153 153 
  s_raster<-brick(infile_covariates)                   #read in the data brck
154 154 
  names(s_raster)<-covar_names               #Assigning names to the raster layers: making sure it is included in the extraction
155 155
156 156 
  #Reading monthly data
157 157 
  dst<-readOGR(dsn=out_path,layer=sub(".shp","",basename(infile_monthly)))
158
159 
  ### TO DO -important ###
160 
  #SCREENING IN COVARIATE SCRIPT AND DATA PREP SCRIPT !!! Only perform predictions here
161 
  #Screen for extreme values": this needs more thought, min and max val vary with regions
162 
  #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)
163 
  #r1[r1 < (min_val)]<-NA
164
158
165 159 
  ########### CREATE SAMPLING -TRAINING AND TESTING STATIONS ###########
166 160
167 161 
  #Input for sampling function...
......
215 209 
    }
216 210 
    clim_yearlist<-list_tmp
217 211 
  }
218
219 
  cat(as.character(t2),file=log_fname,sep="\n", append=TRUE)
220 212 
  t2<-proc.time()-t1
221
213 
  cat(as.character(t2),file=log_fname,sep="\n", append=TRUE)
214 

  		




  222
  215
  
    
  ################## PREDICT AT DAILY TIME SCALE #################


  		




  223
  216
  
    
  #Predict at daily time scale from single time scale or multiple time scale methods: 2 methods availabe now


  		




  224
  217
  
    
  


  		




  ......




  250
  243
  
    
    list_param_run_prediction_gam_daily <-list(i,s_raster,covar_names,lst_avg,list_models,dst,var,y_var_name, sampling_obj,interpolation_method,out_prefix,out_path)


  		




  251
  244
  
    
    names(list_param_run_prediction_gam_daily)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","var","y_var_name","sampling_obj","interpolation_method","out_prefix","out_path")


  		




  252
  245
  
    
    #test <- runGAM_day_fun(1,list_param_run_prediction_gam_daily)


  		




  
  246
  
    
    


  		




  253
  247
  
    
    method_mod_obj<-mclapply(1:length(sampling_obj$ghcn_data_day),list_param=list_param_run_prediction_gam_daily,runGAM_day_fun,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement


  		




  254
  248
  
    
    #method_mod_obj<-mclapply(1:18,list_param=list_param_run_prediction_gam_daily,runGAM_day_fun,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement


  		




  255
  249
  
    
    


  		




  ......




  257
  251
  
    
    


  		




  258
  252
  
    
  }


  		




  259
  253
  
    
  


  		




  
  254
  
    
  if (interpolation_method=="kriging_daily"){


  		




  
  255
  
    
    #input a list:note that ghcn.subsets is not sampling_obj$data_day_ghcn


  		




  
  256
  
    
    i<-1


  		




  
  257
  
    
    list_param_run_prediction_kriging_daily <-list(i,s_raster,covar_names,lst_avg,list_models,dst,var,y_var_name, sampling_obj,interpolation_method,out_prefix,out_path)


  		




  
  258
  
    
    names(list_param_run_prediction_kriging_daily)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","var","y_var_name","sampling_obj","interpolation_method","out_prefix","out_path")


  		




  
  259
  
    
    #test <- runKriging_day_fun(1,list_param_run_prediction_kriging_daily)


  		




  
  260
  
    
    method_mod_obj<-mclapply(1:length(sampling_obj$ghcn_data_day),list_param=list_param_run_prediction_kriging_daily,runKriging_day_fun,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement


  		




  
  261
  
    
    #method_mod_obj<-mclapply(1:18,list_param=list_param_run_prediction_kriging_daily,runKriging_day_fun,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement


  		




  
  262
  
    
    


  		




  
  263
  
    
    save(method_mod_obj,file= file.path(out_path,paste("method_mod_obj_",interpolation_method,"_",y_var_name,out_prefix,".RData",sep="")))


  		




  
  264
  
    
    


  		




  
  265
  
    
  }


  		




  260
  266
  
    
  t2<-proc.time()-t1


  		




  261
  267
  
    
  cat(as.character(t2),file=log_fname,sep="\n", append=TRUE)


  		




  262
  268
  
    
  #browser()


  		




  ......




  321
  327
  
    
    


  		




  322
  328
  
    
  }


  		




  323
  329
  
    
  


  		




  324
  
  
    
  if (interpolation_method=="gam_daily"){


  		




  
  330
  
    
  if (interpolation_method=="gam_daily" | interpolation_method=="kriging_daily" ){


  		




  325
  331
  
    
    raster_prediction_obj<-list(method_mod_obj,validation_mod_obj,tb_diagnostic_v,


  		




  326
  332
  
    
                                summary_metrics_v,summary_month_metrics_v)


  		




  327
  333
  
    
    names(raster_prediction_obj)<-c("method_mod_obj","validation_mod_obj","tb_diagnostic_v",


  		












Also available in:  Unified diff