Revision 87a06ee7
Added by Benoit Parmentier over 11 years ago
climate/research/oregon/interpolation/GAM_fusion_function_multisampling.R | ||
---|---|---|
4 | 4 |
# 1)predict_raster_model<-function(in_models,r_stack,out_filename) |
5 | 5 |
# 2)fit_models<-function(list_formulas,data_training) |
6 | 6 |
# 3)runClimCAI<-function(j) : not working yet |
7 |
# 4)runClim_KGFusion<-function(j,list_param) |
|
8 |
# 5)runGAMFusion <- function(i,list_param) |
|
7 |
# 4)runClim_KGFusion<-function(j,list_param) function for monthly step (climatology) in the fusion method
|
|
8 |
# 5)runGAMFusion <- function(i,list_param) : daily step for fusion method, perform daily prediction
|
|
9 | 9 |
# |
10 | 10 |
#AUTHOR: Benoit Parmentier |
11 |
#DATE: 03/12/2013
|
|
11 |
#DATE: 03/29/2013
|
|
12 | 12 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#363-- |
13 | 13 |
|
14 | 14 |
##Comments and TODO: |
15 | 15 |
#This script is meant to be for general processing tile by tile or region by region. |
16 | 16 |
# Note that the functions are called from GAM_fusion_analysis_raster_prediction_mutlisampling.R. |
17 | 17 |
# This will be expanded to other methods. |
18 |
|
|
18 |
# Change name of output tif to include the variable!!! (TMIN or TMAX) |
|
19 | 19 |
################################################################################################## |
20 | 20 |
|
21 | 21 |
|
... | ... | |
111 | 111 |
data_month$y_var<-data_month$TMax #Adding TMax as the variable modeled |
112 | 112 |
} |
113 | 113 |
if (var=="TMIN"){ |
114 |
data_month$y_var<-data_month$TMin #Adding TMi as the variable modeled |
|
114 |
data_month$y_var<-data_month$TMin #Adding TMin as the variable modeled
|
|
115 | 115 |
} |
116 | 116 |
|
117 | 117 |
mod_list<-fit_models(list_formulas,data_month) #only gam at this stage |
... | ... | |
359 | 359 |
dst<-list_param$dst #monthly station dataset |
360 | 360 |
|
361 | 361 |
########## |
362 |
# STEP 1 - interpolate delta across space
|
|
363 |
############# Read in information and get traiing and testing stations
|
|
362 |
# STEP 1 - Read in information and get traing and testing stations
|
|
363 |
############# |
|
364 | 364 |
|
365 | 365 |
date<-strptime(sampling_dat$date[i], "%Y%m%d") # interpolation date being processed |
366 | 366 |
month<-strftime(date, "%m") # current month of the date being processed |
... | ... | |
400 | 400 |
modst$LSTD_bias <- modst$LST-modst$TMax; #That is the difference between the monthly LST mean and monthly station mean |
401 | 401 |
} |
402 | 402 |
#This may be unnecessary since LSTD_bias is already in dst?? check the info |
403 |
#Some loss of observations: LSTD_bias for January has only 56 out of 66 possible TMIN!!! We may need to look into this issue |
|
404 |
#to avoid some losses of station data... |
|
403 | 405 |
|
404 | 406 |
#Clearn out this part: make this a function call |
405 | 407 |
x<-as.data.frame(data_v) |
406 | 408 |
d<-as.data.frame(data_s) |
407 |
#x[x$value==-999.9]<-NA |
|
408 | 409 |
for (j in 1:nrow(x)){ |
409 | 410 |
if (x$value[j]== -999.9){ |
410 | 411 |
x$value[j]<-NA |
... | ... | |
415 | 416 |
d$value[j]<-NA |
416 | 417 |
} |
417 | 418 |
} |
418 |
#x[x$value==-999.9]<-NA |
|
419 |
#d[d$value==-999.9]<-NA |
|
420 | 419 |
pos<-match("value",names(d)) #Find column with name "value" |
421 | 420 |
#names(d)[pos]<-c("dailyTmax") |
422 | 421 |
names(d)[pos]<-y_var_name |
422 |
pos<-match("value",names(x)) #Find column with name "value" |
|
423 | 423 |
names(x)[pos]<-y_var_name |
424 |
#names(x)[pos]<-c("dailyTmax") |
|
425 |
pos<-match("station",names(d)) #Find column with name "value" |
|
424 |
pos<-match("station",names(d)) #Find column with station ID |
|
426 | 425 |
names(d)[pos]<-c("id") |
426 |
pos<-match("station",names(x)) #Find column with name station ID |
|
427 | 427 |
names(x)[pos]<-c("id") |
428 |
names(modst)[1]<-c("id") #modst contains the average tmax per month for every stations... |
|
428 |
pos<-match("station",names(modst)) #Find column with name station ID |
|
429 |
names(modst)[pos]<-c("id") #modst contains the average tmax per month for every stations... |
|
429 | 430 |
|
430 | 431 |
dmoday <-merge(modst,d,by="id",suffixes=c("",".y2")) |
431 | 432 |
xmoday <-merge(modst,x,by="id",suffixes=c("",".y2")) |
432 |
mod_pat<-glob2rx("*.y2") |
|
433 |
mod_pat<-glob2rx("*.y2") #remove duplicate columns that have ".y2" in their names
|
|
433 | 434 |
var_pat<-grep(mod_pat,names(dmoday),value=FALSE) # using grep with "value" extracts the matching names |
434 |
dmoday<-dmoday[,-var_pat] |
|
435 |
dmoday<-dmoday[,-var_pat] #dropping relevant columns
|
|
435 | 436 |
mod_pat<-glob2rx("*.y2") |
436 | 437 |
var_pat<-grep(mod_pat,names(xmoday),value=FALSE) # using grep with "value" extracts the matching names |
437 | 438 |
xmoday<-xmoday[,-var_pat] #Removing duplicate columns |
438 | 439 |
|
439 | 440 |
data_v<-xmoday |
440 | 441 |
|
441 |
#dmoday contains the daily tmax values for training with TMax being the monthly station tmax mean
|
|
442 |
#xmoday contains the daily tmax values for validation with TMax being the monthly station tmax mean
|
|
442 |
#dmoday contains the daily tmax values for training with TMax/TMin being the monthly station tmax/tmin mean
|
|
443 |
#xmoday contains the daily tmax values for validation with TMax/TMin being the monthly station tmax/tmin mean
|
|
443 | 444 |
|
444 | 445 |
########## |
445 | 446 |
# STEP 3 - interpolate daily delta across space |
... | ... | |
447 | 448 |
|
448 | 449 |
#Change to take into account TMin and TMax |
449 | 450 |
if (var=="TMIN"){ |
450 |
daily_delta<-dmoday$dailyTmin-dmoday$TMin |
|
451 |
daily_delta<-dmoday$dailyTmin-dmoday$TMin #daily detl is the difference between monthly and daily temperatures
|
|
451 | 452 |
} |
452 | 453 |
if (var=="TMAX"){ |
453 | 454 |
daily_delta<-dmoday$dailyTmax-dmoday$TMax |
... | ... | |
472 | 473 |
#Saving kriged surface in raster images |
473 | 474 |
data_name<-paste("daily_delta_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
474 | 475 |
"_",sampling_dat$run_samp[i],sep="") |
475 |
raster_name_delta<-paste("fusion_",data_name,out_prefix,".tif", sep="") |
|
476 |
raster_name_delta<-paste("fusion_",var,"_",data_name,out_prefix,".tif", sep="")
|
|
476 | 477 |
writeRaster(daily_delta_rast, filename=raster_name_delta,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
477 | 478 |
|
478 | 479 |
#Now predict daily after having selected the relevant month |
... | ... | |
506 | 507 |
|
507 | 508 |
obj_names<-c(y_var_name,"clim","delta","data_s","data_v", |
508 | 509 |
"sampling_dat",model_name) |
509 |
names(delta_obj)<-obj_names |
|
510 |
save(delta_obj,file= paste("delta_obj_",sampling_dat$date[i],"_",sampling_dat$prop[i], |
|
510 |
names(delta_obj)<-obj_names #add TMIN or TMAX name in saving obj
|
|
511 |
save(delta_obj,file= paste("delta_obj_",var,"_",sampling_dat$date[i],"_",sampling_dat$prop[i],
|
|
511 | 512 |
"_",sampling_dat$run_samp[i],out_prefix,".RData",sep="")) |
512 | 513 |
return(delta_obj) |
513 | 514 |
|
Also available in: Unified diff
raster prediction and gam fusion modifications for TMIN predictions