Revision 33095a53
Added by Benoit Parmentier over 11 years ago
climate/research/oregon/interpolation/GAM_fusion_analysis_raster_prediction_multisampling.R | ||
---|---|---|
11 | 11 |
#5)GAM fusion: 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: 02/26/2013
|
|
14 |
#DATE: 03/05/2013
|
|
15 | 15 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#568-- |
16 | 16 |
################################################################################################### |
17 | 17 |
|
... | ... | |
30 | 30 |
library(plotrix) |
31 | 31 |
library(maptools) |
32 | 32 |
|
33 |
### Parameters and argument |
|
34 |
|
|
35 |
infile2<-"list_365_dates_04212012.txt" |
|
33 |
### Parameters and arguments |
|
36 | 34 |
|
37 | 35 |
## output param from previous script: Database_stations_covariates_processing_function |
38 |
infile_monthly<-"monthly_covariates_ghcn_data_TMAXy2010_2010_VE_02082013.shp" |
|
39 |
infile_daily<-"daily_covariates_ghcn_data_TMAXy2010_2010_VE_02082013.shp" |
|
40 |
infile_locs<-"stations_venezuela_region_y2010_2010_VE_02082013.shp" |
|
41 |
infile3<-"covariates__venezuela_region__VE_01292013.tif" #this is an output from covariate script
|
|
36 |
#infile_monthly<-"monthly_covariates_ghcn_data_TMAXy2010_2010_VE_02082013.shp"
|
|
37 |
#infile_daily<-"daily_covariates_ghcn_data_TMAXy2010_2010_VE_02082013.shp"
|
|
38 |
#infile_locs<-"stations_venezuela_region_y2010_2010_VE_02082013.shp"
|
|
39 |
infile_covariates<-"covariates__venezuela_region__VE_01292013.tif" #this is an output from covariate script
|
|
42 | 40 |
var<-"TMAX" |
43 | 41 |
out_prefix<-"_365d_GAM_fus5_all_lstd_02202013" #User defined output prefix |
44 | 42 |
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84: same as earlier |
45 |
#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"; |
|
46 |
#infile_monthly<-"monthly_covariates_ghcn_data_TMAXy2010_2010_VE_02082013.shp" #outile4 from database_covar script |
|
47 |
#infile_daily<-"daily_covariates_ghcn_data_TMAXy2010_2010_VE_02082013.shp" #outfile3 from database_covar script |
|
48 |
#infile_locs<-"stations_venezuela_region_y2010_2010_VE_02082013.shp" #outfile2? from database covar script |
|
49 | 43 |
|
44 |
infile_monthly<-list_outfiles$monthly_covar_ghcn_data #outile4 from database_covar script |
|
45 |
infile_daily<-list_outfiles$daily_covar_ghcn_data #outfile3 from database_covar script |
|
46 |
infile_locs<- list_outfiles$loc_stations_ghcn #outfile2? from database covar script |
|
47 |
rnames <-c("x","y","lon","lat","N","E","N_w","E_w","elev","slope","aspect","CANHEIGHT","DISTOC") |
|
48 |
lc_names<-c("LC1","LC2","LC3","LC4","LC5","LC6","LC7","LC8","LC9","LC10","LC11","LC12") |
|
49 |
lst_names<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12", |
|
50 |
"nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08", |
|
51 |
"nobs_09","nobs_10","nobs_11","nobs_12") |
|
52 |
covar_names<-c(rnames,lc_names,lst_names) |
|
50 | 53 |
### |
51 |
|
|
52 |
if (var=="TMAX"){ |
|
53 |
y_var_name<-"dailyTmax" |
|
54 |
} |
|
55 |
if (var=="TMIN"){ |
|
56 |
y_var_name<-"dailyTmin" |
|
57 |
} |
|
58 |
|
|
59 | 54 |
#Input for sampling function... |
60 | 55 |
seed_number<- 100 #if seed zero then no seed? |
61 | 56 |
nb_sample<-1 #number of time random sampling must be repeated for every hold out proportion |
62 |
prop_min<-0.3 #if prop_min=prop_max and step=0 then predicitons are done for the number of dates... |
|
63 |
prop_max<-0.3 |
|
64 | 57 |
step<-0 |
65 | 58 |
constant<-0 #if value 1 then use the same samples as date one for the all set of dates |
59 |
prop_minmax<-c(0.3,0.3) #if prop_min=prop_max and step=0 then predicitons are done for the number of dates... |
|
60 |
infile_dates<-"list_365_dates_04212012.txt" |
|
66 | 61 |
|
67 | 62 |
#Models to run...this can be change for each run |
68 | 63 |
list_models<-c("y_var ~ s(elev_1)", |
... | ... | |
87 | 82 |
script_path<-"/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/" |
88 | 83 |
setwd(in_path) |
89 | 84 |
|
90 |
source(file.path(script_path,"GAM_fusion_function_multisampling_02262013.R")) |
|
85 |
|
|
86 |
list_param_data_prep<-c(infile_monthly,infile_daily,infile_locs,infile_covariates,var,out_prefix,CRS_locs_WGS84) |
|
87 |
list_param_raster_prediction<-c(list_param_data_prep, |
|
88 |
seed_number,nb_sample,step,constant,prop_minmax,infile_dates, |
|
89 |
list_models,lst_avg,in_path,out_path,script_path, |
|
90 |
interpolation_method) |
|
91 |
|
|
92 |
|
|
93 |
source(file.path(script_path,"sampling_script_functions_03052013.R")) |
|
94 |
source(file.path(script_path,"GAM_fusion_function_multisampling_03052013.R")) |
|
91 | 95 |
source(file.path(script_path,"GAM_fusion_function_multisampling_validation_metrics_02262013.R")) |
92 | 96 |
|
97 |
|
|
93 | 98 |
###################### START OF THE SCRIPT ######################## |
94 | 99 |
|
100 |
|
|
101 |
if (var=="TMAX"){ |
|
102 |
y_var_name<-"dailyTmax" |
|
103 |
} |
|
104 |
if (var=="TMIN"){ |
|
105 |
y_var_name<-"dailyTmin" |
|
106 |
} |
|
107 |
|
|
95 | 108 |
################# CREATE LOG FILE ##################### |
96 | 109 |
|
97 | 110 |
#create log file to keep track of details such as processing times and parameters. |
98 | 111 |
|
99 |
log_fname<-paste("R_log_t",out_prefix, ".log",sep="")
|
|
112 |
log_fname<-paste("R_log_raster_prediction",out_prefix, ".log",sep="")
|
|
100 | 113 |
|
101 | 114 |
if (file.exists(log_fname)){ #Stop the script??? |
102 | 115 |
file.remove(log_fname) |
... | ... | |
112 | 125 |
writeLines("Starting script process time:",con=log_file,sep="\n") |
113 | 126 |
writeLines(as.character(time1),con=log_file,sep="\n") |
114 | 127 |
|
115 |
############### Reading the daily station data and other data #################
|
|
128 |
############### READING INPUTS: DAILY STATION DATA AND OTEHR DATASETS #################
|
|
116 | 129 |
|
117 |
ghcn<-readOGR(dsn=in_path,layer=sub(".shp","",infile_daily))
|
|
130 |
ghcn<-readOGR(dsn=in_path,layer=sub(".shp","",basename(infile_daily)))
|
|
118 | 131 |
CRS_interp<-proj4string(ghcn) #Storing projection information (ellipsoid, datum,etc.) |
119 |
stat_loc<-readOGR(dsn=in_path,layer=sub(".shp","",infile_locs))
|
|
120 |
dates <-readLines(file.path(in_path,infile2)) #dates to be predicted
|
|
132 |
stat_loc<-readOGR(dsn=in_path,layer=sub(".shp","",basename(infile_locs)))
|
|
133 |
dates <-readLines(file.path(in_path,infile_dates)) #dates to be predicted
|
|
121 | 134 |
|
122 | 135 |
#Reading of covariate brick covariates can be changed... |
123 |
rnames <-c("x","y","lon","lat","N","E","N_w","E_w","elev","slope","aspect","CANHEIGHT","DISTOC") |
|
124 |
lc_names<-c("LC1","LC2","LC3","LC4","LC5","LC6","LC7","LC8","LC9","LC10","LC11","LC12") |
|
125 |
lst_names<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12", |
|
126 |
"nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08", |
|
127 |
"nobs_09","nobs_10","nobs_11","nobs_12") |
|
128 |
covar_names<-c(rnames,lc_names,lst_names) |
|
129 |
s_raster<-brick(infile3) #read in the data brck |
|
136 |
|
|
137 |
s_raster<-brick(infile_covariates) #read in the data brck |
|
130 | 138 |
names(s_raster)<-covar_names #Assigning names to the raster layers: making sure it is included in the extraction |
139 |
pos<-match("elev",names(s_raster)) |
|
140 |
names(s_raster)[pos]<-"elev_1" |
|
141 |
|
|
142 |
#Screen for extreme values": this needs more thought, min and max val vary with regions |
|
143 |
#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) |
|
144 |
#r1[r1 < (min_val)]<-NA |
|
131 | 145 |
|
132 | 146 |
#Reading monthly data |
133 |
data3<-readOGR(dsn=in_path,layer=sub(".shp","",infile_monthly))
|
|
147 |
data3<-readOGR(dsn=in_path,layer=sub(".shp","",basename(infile_monthly)))
|
|
134 | 148 |
dst_all<-data3 |
135 | 149 |
dst<-data3 |
136 | 150 |
|
137 | 151 |
### TO DO -important ### |
138 |
#Cleaning/sceerniging functions for daily stations, monthly stations and covariates?? do this at the preparation stage!!! |
|
139 |
## |
|
140 |
|
|
141 |
##### Create sampling: select training and testing sites ### |
|
142 |
#Make this a a function!!!! |
|
143 |
|
|
144 |
if (seed_number>0) { |
|
145 |
set.seed(seed_number) #Using a seed number allow results based on random number to be compared... |
|
146 |
} |
|
147 |
nel<-length(dates) |
|
148 |
dates_list<-vector("list",nel) #list of one row data.frame |
|
149 |
|
|
150 |
prop_range<-(seq(from=prop_min,to=prop_max,by=step))*100 #range of proportion to run |
|
151 |
sn<-length(dates)*nb_sample*length(prop_range) #Number of samples to run |
|
152 |
|
|
153 |
for(i in 1:length(dates)){ |
|
154 |
d_tmp<-rep(dates[i],nb_sample*length(prop_range)) #repeating same date |
|
155 |
s_nb<-rep(1:nb_sample,length(prop_range)) #number of random sample per proportion |
|
156 |
prop_tmp<-sort(rep(prop_range, nb_sample)) |
|
157 |
tab_run_tmp<-cbind(d_tmp,s_nb,prop_tmp) |
|
158 |
dates_list[[i]]<-tab_run_tmp |
|
159 |
} |
|
152 |
#Cleaning/sceerniging functions for daily stations, monthly stations and covariates?? do this during the preparation stage!!!?? |
|
153 |
### |
|
160 | 154 |
|
161 |
sampling_dat<-as.data.frame(do.call(rbind,dates_list)) |
|
162 |
names(sampling_dat)<-c("date","run_samp","prop") |
|
155 |
########### CREATE SAMPLING -TRAINING AND TESTING STATIONS ########### |
|
163 | 156 |
|
164 |
for(i in 2:3){ # start of the for loop #1 |
|
165 |
sampling_dat[,i]<-as.numeric(as.character(sampling_dat[,i])) |
|
166 |
} |
|
157 |
#Input for sampling function... |
|
167 | 158 |
|
168 |
sampling_dat$date<- as.character(sampling_dat[,1]) |
|
169 |
#ghcn.subsets <-lapply(dates, function(d) subset(ghcn, date==d)) #this creates a list of 10 or 365 subsets dataset based on dates |
|
170 |
ghcn.subsets <-lapply(as.character(sampling_dat$date), function(d) subset(ghcn, date==d)) #this creates a list of 10 or 365 subsets dataset based on dates |
|
159 |
#dates #list of dates for prediction |
|
160 |
#ghcn_name<-"ghcn" #infile daily data |
|
171 | 161 |
|
172 |
#Make this a function?? |
|
173 |
## adding choice of constant sample |
|
174 |
if (seed_number>0) { |
|
175 |
set.seed(seed_number) #Using a seed number allow results based on random number to be compared... |
|
176 |
} |
|
162 |
list_param_sampling<-list(seed_number,nb_sample,step,constant,prop_minmax,dates,ghcn) |
|
163 |
#list_param_sampling<-list(seed_number,nb_sample,step,constant,prop_minmax,dates,ghcn_name) |
|
164 |
names(list_param_sampling)<-c("seed_number","nb_sample","step","constant","prop_minmax","dates","ghcn") |
|
165 |
|
|
166 |
#run function |
|
167 |
sampling_obj<-sampling_training_testing(list_param_sampling) |
|
177 | 168 |
|
178 |
sampling<-vector("list",length(ghcn.subsets)) |
|
179 |
sampling_station_id<-vector("list",length(ghcn.subsets)) |
|
180 |
for(i in 1:length(ghcn.subsets)){ |
|
181 |
n<-nrow(ghcn.subsets[[i]]) |
|
182 |
prop<-(sampling_dat$prop[i])/100 |
|
183 |
ns<-n-round(n*prop) #Create a sample from the data frame with 70% of the rows |
|
184 |
nv<-n-ns #create a sample for validation with prop of the rows |
|
185 |
ind.training <- sample(nrow(ghcn.subsets[[i]]), size=ns, replace=FALSE) #This selects the index position for 70% of the rows taken randomly |
|
186 |
ind.testing <- setdiff(1:nrow(ghcn.subsets[[i]]), ind.training) |
|
187 |
#Find the corresponding |
|
188 |
data_sampled<-ghcn.subsets[[i]][ind.training,] #selected the randomly sampled stations |
|
189 |
station_id.training<-data_sampled$station #selected id for the randomly sampled stations (115) |
|
190 |
#Save the information |
|
191 |
sampling[[i]]<-ind.training |
|
192 |
sampling_station_id[[i]]<- station_id.training |
|
193 |
} |
|
194 |
## Use same samples across the year... |
|
195 |
if (constant==1){ |
|
196 |
sampled<-sampling[[1]] |
|
197 |
data_sampled<-ghcn.subsets[[1]][sampled,] #selected the randomly sampled stations |
|
198 |
station_sampled<-data_sampled$station #selected id for the randomly sampled stations (115) |
|
199 |
list_const_sampling<-vector("list",sn) |
|
200 |
list_const_sampling_station_id<-vector("list",sn) |
|
201 |
for(i in 1:sn){ |
|
202 |
station_id.training<-intersect(station_sampled,ghcn.subsets[[i]]$station) |
|
203 |
ind.training<-match(station_id.training,ghcn.subsets[[i]]$station) |
|
204 |
list_const_sampling[[i]]<-ind.training |
|
205 |
list_const_sampling_station_id[[i]]<-station_id.training |
|
206 |
} |
|
207 |
sampling<-list_const_sampling |
|
208 |
sampling_station_id<-list_const_sampling_station_id |
|
209 |
} |
|
210 |
#return() |
|
211 | 169 |
|
212 |
########### PREDICT FOR MONTHLY SCALE ############# |
|
170 |
########### PREDICT FOR MONTHLY SCALE ##################
|
|
213 | 171 |
|
214 | 172 |
#First predict at the monthly time scale: climatology |
215 | 173 |
writeLines("Predictions at monthly scale:",con=log_file,sep="\n") |
216 | 174 |
t1<-proc.time() |
217 | 175 |
gamclim_fus_mod<-mclapply(1:12, runClim_KGFusion,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement |
218 |
#gamclim_fus_mod<-mclapply(1:12, runClim_KGFusion,mc.preschedule=FALSE,mc.cores = 4) #This is the end bracket from mclapply(...) statement
|
|
176 |
#gamclim_fus_mod<-mclapply(1:6, runClim_KGFusion,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement
|
|
219 | 177 |
save(gamclim_fus_mod,file= paste("gamclim_fus_mod",out_prefix,".RData",sep="")) |
220 | 178 |
t2<-proc.time()-t1 |
221 | 179 |
writeLines(as.character(t2),con=log_file,sep="\n") |
... | ... | |
231 | 189 |
################## PREDICT AT DAILY TIME SCALE ################# |
232 | 190 |
|
233 | 191 |
#put together list of clim models per month... |
234 |
rast_clim_yearlist<-list_tmp |
|
192 |
#rast_clim_yearlist<-list_tmp |
|
193 |
clim_yearlist<-list_tmp |
|
235 | 194 |
#Second predict at the daily time scale: delta |
236 | 195 |
|
237 | 196 |
#gam_fus_mod<-mclapply(1:1, runGAMFusion,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement |
238 | 197 |
writeLines("Predictions at the daily scale:",con=log_file,sep="\n") |
239 | 198 |
t1<-proc.time() |
240 |
gam_fus_mod<-mclapply(1:length(ghcn.subsets), runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement |
|
199 |
|
|
200 |
#input a list:note that ghcn.subsets is not sampling_obj$data_day_ghcn |
|
201 |
list_param_runGAMFusion<-list(i,clim_yearlist,sampling_obj,dst,var,y_var_name, out_prefix) |
|
202 |
names(list_param_runGAMFusion)<-c("list_index","clim_yearlist","sampling_obj","dst","var","y_var_name","out_prefix") |
|
203 |
#test<-mclapply(1:18, runGAMFusion,list_param=list_param_runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) |
|
204 |
|
|
205 |
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 |
|
206 |
|
|
207 |
#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 |
|
208 |
#gam_fus_mod<-mclapply(1:length(ghcn.subsets), runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement |
|
241 | 209 |
save(gam_fus_mod,file= paste("gam_fus_mod",out_prefix,".RData",sep="")) |
242 | 210 |
t2<-proc.time()-t1 |
243 | 211 |
writeLines(as.character(t2),con=log_file,sep="\n") |
Also available in: Unified diff
Raster prediction, now calling sampling function and cleaning out input parameters