Revision a96491e0
Added by Benoit Parmentier almost 12 years ago
climate/research/oregon/interpolation/GAM_fusion_function_multisampling.R | ||
---|---|---|
12 | 12 |
mod <-in_models[[i]] #accessing GAM model ojbect "j" |
13 | 13 |
raster_name<-out_filename[[i]] |
14 | 14 |
if (inherits(mod,"gam")) { #change to c("gam","autoKrige") |
15 |
raster_pred<- predict(object=s_raster,model=mod,na.rm=FALSE) #Using the coeff to predict new values.
|
|
15 |
raster_pred<- predict(object=r_stack,model=mod,na.rm=FALSE) #Using the coeff to predict new values.
|
|
16 | 16 |
names(raster_pred)<-"y_pred" |
17 | 17 |
writeRaster(raster_pred, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
18 |
print(paste("Interpolation:","mod", j ,sep=" ")) |
|
18 |
#print(paste("Interpolation:","mod", j ,sep=" "))
|
|
19 | 19 |
list_rast_pred[[i]]<-raster_name |
20 | 20 |
} |
21 | 21 |
} |
... | ... | |
132 | 132 |
# |
133 | 133 |
|
134 | 134 |
runClim_KGFusion<-function(j){ |
135 |
|
|
135 | 136 |
#Make this a function with multiple argument that can be used by mcmapply?? |
136 | 137 |
#This creates clim fusion layers... |
138 |
#Parameters: |
|
139 |
|
|
140 |
#1)s_raster: brick of covariates : could pass as argument only specific variables?? |
|
141 |
#2)dst: monthly data (infile_monthly): could pass only the subset from the month?? |
|
142 |
#3)list_models |
|
143 |
#4)brick names covarnames |
|
144 |
#5)out_prefix |
|
145 |
|
|
137 | 146 |
|
147 |
## STEP 1: PARSE PARAMETERS AND ARGUMENTS |
|
138 | 148 |
|
139 | 149 |
#Model and response variable can be changed without affecting the script |
140 | 150 |
prop_month<-0 #proportion retained for validation |
141 | 151 |
run_samp<-1 |
142 | 152 |
|
143 |
list_formulas<-lapply(list_models,as.formula,env=.GlobalEnv) #mulitple arguments passed to lapply!!
|
|
144 |
|
|
153 |
## STEP 2: PREPARE DATA
|
|
154 |
|
|
145 | 155 |
data_month<-dst[dst$month==j,] #Subsetting dataset for the relevant month of the date being processed |
146 | 156 |
LST_name<-lst_avg[j] # name of LST month to be matched |
147 | 157 |
data_month$LST<-data_month[[LST_name]] |
148 | 158 |
|
149 |
#LST bias to model... |
|
150 |
data_month$LSTD_bias<-data_month$LST-data_month$TMax |
|
151 |
data_month$y_var<-data_month$LSTD_bias #Adding bias as the variable modeled |
|
152 |
mod_list<-fit_models(list_formulas,data_month) #only gam at this stage |
|
153 |
cname<-paste("mod",1:length(mod_list),sep="") #change to more meaningful name? |
|
154 |
names(mod_list)<-cname |
|
155 | 159 |
#Adding layer LST to the raster stack |
156 |
pos<-match("elev",names(s_raster)) |
|
157 |
layerNames(s_raster)[pos]<-"elev_1" |
|
158 |
|
|
160 |
covar_rast<-s_raster |
|
159 | 161 |
pos<-match("LST",names(s_raster)) #Find the position of the layer with name "LST", if not present pos=NA |
160 | 162 |
s_raster<-dropLayer(s_raster,pos) # If it exists drop layer |
161 | 163 |
LST<-subset(s_raster,LST_name) |
162 | 164 |
names(LST)<-"LST" |
163 |
#Screen for extreme values": this needs more thought, min and max val vary with regions |
|
164 |
#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) |
|
165 |
#r1[r1 < (min_val)]<-NA |
|
166 | 165 |
s_raster<-addLayer(s_raster,LST) #Adding current month |
167 | 166 |
|
167 |
#LST bias to model... |
|
168 |
#if (var==TMAX): will need to modify to take into account TMAX, TMIN and LST_day,LST_night |
|
169 |
data_month$LSTD_bias<-data_month$LST-data_month$TMax |
|
170 |
data_month$y_var<-data_month$LSTD_bias #Adding bias as the variable modeled |
|
171 |
|
|
172 |
#### STEP3: NOW FIT AND PREDICT MODEL |
|
173 |
|
|
174 |
list_formulas<-lapply(list_models,as.formula,env=.GlobalEnv) #mulitple arguments passed to lapply!! |
|
175 |
|
|
176 |
mod_list<-fit_models(list_formulas,data_month) #only gam at this stage |
|
177 |
cname<-paste("mod",1:length(mod_list),sep="") #change to more meaningful name? |
|
178 |
names(mod_list)<-cname |
|
179 |
|
|
168 | 180 |
#Now generate file names for the predictions... |
169 | 181 |
list_out_filename<-vector("list",length(mod_list)) |
170 | 182 |
names(list_out_filename)<-cname |
... | ... | |
195 | 207 |
writeRaster(clim_fus_rast, filename=raster_name,overwrite=TRUE) #Wri |
196 | 208 |
} |
197 | 209 |
|
198 |
#Adding Kriging for Climatology options |
|
210 |
#### STEP 4:Adding Kriging for Climatology options
|
|
199 | 211 |
|
200 | 212 |
bias_xy<-coordinates(data_month) |
201 | 213 |
fitbias<-Krig(bias_xy,data_month$LSTD_bias,theta=1e5) #use TPS or krige |
202 | 214 |
mod_krtmp1<-fitbias |
203 | 215 |
model_name<-"mod_kr" |
216 |
|
|
204 | 217 |
|
205 | 218 |
bias_rast<-interpolate(LST,fitbias) #interpolation using function from raster package |
206 | 219 |
#Saving kriged surface in raster images |
... | ... | |
221 | 234 |
rast_bias_list[[model_name]]<-raster_name_bias |
222 | 235 |
rast_clim_list[[model_name]]<-raster_name_clim |
223 | 236 |
|
224 |
#Prepare object to return |
|
237 |
#### STEP 5: Prepare object and return |
|
238 |
|
|
225 | 239 |
clim_obj<-list(rast_bias_list,rast_clim_list,data_month,mod_list,list_formulas) |
226 | 240 |
names(clim_obj)<-c("bias","clim","data_month","mod","formulas") |
227 | 241 |
|
... | ... | |
232 | 246 |
|
233 | 247 |
## Run function for kriging...? |
234 | 248 |
|
235 |
runGAMFusion <- function(i) { # loop over dates |
|
236 |
#Change this to allow explicitly arguments... |
|
249 |
#runGAMFusion <- function(i) { # loop over dates |
|
250 |
runGAMFusion <- function(i,list_param) { # loop over dates |
|
251 |
#### Change this to allow explicitly arguments... |
|
237 | 252 |
#Arguments: |
238 | 253 |
#1)list of climatology files for all models...(12*nb of models) |
239 |
#2)data_s:training |
|
240 |
#3)data_v:testing |
|
241 |
#4)list of dates?? |
|
242 |
#5)stack of covariates: not needed at this this stage |
|
254 |
#2)sampling_obj$data_day_gcn: ghcn.subsets (data per date ) |
|
255 |
#4)sampling_dat: list of sampling information for every run (proporation, month,sample) |
|
256 |
#5)ampling_index : list of training |
|
243 | 257 |
#6)dst: data at the monthly time scale |
258 |
#7)var: |
|
259 |
#8)y_var_name: |
|
260 |
#9)out_prefix |
|
261 |
|
|
262 |
### PARSING INPUT ARGUMENTS |
|
263 |
|
|
264 |
#list_param_runGAMFusion<-list(i,clim_yearlist,sampling_obj,var,y_var_name, out_prefix) |
|
265 |
rast_clim_yearlist<-list_param$clim_yearlist |
|
266 |
sampling_obj<-list_param$sampling_obj |
|
267 |
ghcn.subsets<-sampling_obj$ghcn_data_day |
|
268 |
sampling_dat <- sampling_obj$sampling_dat |
|
269 |
sampling <- sampling_obj$sampling_index |
|
270 |
var<-list_param$var |
|
271 |
y_var_name<-list_param$y_var_name |
|
272 |
out_prefix<-list_param$out_prefix |
|
273 |
dst<-list_param$dst #monthly station dataset |
|
244 | 274 |
|
245 |
#Function used in the script |
|
275 |
########## |
|
276 |
# STEP 1 - interpolate delta across space |
|
277 |
############# Read in information and get traiing and testing stations |
|
246 | 278 |
|
247 | 279 |
date<-strptime(sampling_dat$date[i], "%Y%m%d") # interpolation date being processed |
248 | 280 |
month<-strftime(date, "%m") # current month of the date being processed |
249 | 281 |
LST_month<-paste("mm_",month,sep="") # name of LST month to be matched |
250 |
proj_str<-proj4string(dst) |
|
282 |
proj_str<-proj4string(dst) #get the local projection information from monthly data
|
|
251 | 283 |
|
252 | 284 |
###Regression part 1: Creating a validation dataset by creating training and testing datasets |
253 | 285 |
data_day<-ghcn.subsets[[i]] |
... | ... | |
269 | 301 |
day<-as.integer(strftime(date_proc, "%d")) |
270 | 302 |
year<-as.integer(strftime(date_proc, "%Y")) |
271 | 303 |
|
304 |
########## |
|
305 |
# STEP 2 - JOIN DAILY AND MONTHLY STATION INFORMATION |
|
306 |
########## |
|
307 |
|
|
272 | 308 |
modst<-dst[dst$month==mo,] #Subsetting dataset for the relevant month of the date being processed |
273 | 309 |
#Change to y_var...could be TMin |
274 | 310 |
#modst$LSTD_bias <- modst$LST-modst$y_var |
275 | 311 |
modst$LSTD_bias <- modst$LST-modst$TMax; #That is the difference between the monthly LST mean and monthly station mean |
276 |
|
|
312 |
|
|
313 |
#Clearn out this part: make this a function call |
|
277 | 314 |
x<-as.data.frame(data_v) |
278 | 315 |
d<-as.data.frame(data_s) |
279 | 316 |
#x[x$value==-999.9]<-NA |
... | ... | |
314 | 351 |
#xmoday contains the daily tmax values for validation with TMax being the monthly station tmax mean |
315 | 352 |
|
316 | 353 |
########## |
317 |
# STEP 7 - interpolate delta across space
|
|
354 |
# STEP 3 - interpolate daily delta across space
|
|
318 | 355 |
########## |
319 | 356 |
|
357 |
#Change to take into account TMin and TMax |
|
320 | 358 |
daily_delta<-dmoday$dailyTmax-dmoday$TMax |
321 | 359 |
daily_delta_xy<-as.matrix(cbind(dmoday$x,dmoday$y)) |
322 | 360 |
fitdelta<-Krig(daily_delta_xy,daily_delta,theta=1e5) #use TPS or krige |
... | ... | |
326 | 364 |
data_s$daily_delta<-daily_delta |
327 | 365 |
|
328 | 366 |
######### |
329 |
# STEP 8 - assemble final answer - T=LST+Bias(interpolated)+delta(interpolated)
|
|
367 |
# STEP 4 - Calculate daily predictions - T(day) = clim(month) + delta(day)
|
|
330 | 368 |
######### |
331 | 369 |
|
332 | 370 |
rast_clim_list<-rast_clim_yearlist[[mo]] #select relevant month |
... | ... | |
354 | 392 |
temp_list[[j]]<-raster_name |
355 | 393 |
} |
356 | 394 |
|
395 |
########## |
|
396 |
# STEP 5 - Prepare output object to return |
|
397 |
########## |
|
398 |
|
|
357 | 399 |
mod_krtmp2<-fitdelta |
358 | 400 |
model_name<-paste("mod_kr","day",sep="_") |
359 | 401 |
names(temp_list)<-names(rast_clim_list) |
Also available in: Unified diff
GAM fusion, modification of runGAMFusion function call to allow mulitple parameters explicitly