Revision 5588b17c
Added by Benoit Parmentier over 11 years ago
climate/research/oregon/interpolation/GAM_fusion_function_multisampling.R | ||
---|---|---|
14 | 14 |
#This functions performs predictions on a raster grid given input models. |
15 | 15 |
#Arguments: list of fitted models, raster stack of covariates |
16 | 16 |
#Output: spatial grid data frame of the subset of tiles |
17 |
#s_sgdf<-as(r_stack,"SpatialGridDataFrame") #Conversion to spatial grid data frame |
|
18 | 17 |
list_rast_pred<-vector("list",length(in_models)) |
19 | 18 |
for (i in 1:length(in_models)){ |
20 | 19 |
mod <-in_models[[i]] #accessing GAM model ojbect "j" |
21 | 20 |
raster_name<-out_filename[[i]] |
22 |
if (inherits(mod,"gam")) { |
|
23 |
#rpred<- predict(mod, newdata=s_sgdf, se.fit = TRUE) #Using the coeff to predict new values. |
|
24 |
#rast_pred2<- predict(object=s_raster,model=mod,na.rm=TRUE) #Using the coeff to predict new values. |
|
21 |
if (inherits(mod,"gam")) { #change to c("gam","autoKrige") |
|
25 | 22 |
raster_pred<- predict(object=s_raster,model=mod,na.rm=FALSE) #Using the coeff to predict new values. |
26 | 23 |
names(raster_pred)<-"y_pred" |
27 | 24 |
writeRaster(raster_pred, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI) |
... | ... | |
30 | 27 |
} |
31 | 28 |
} |
32 | 29 |
if (inherits(mod,"try-error")) { |
33 |
print(paste("no gam model fitted:",mod[1],sep=" ")) |
|
30 |
print(paste("no gam model fitted:",mod[1],sep=" ")) #change message for any model type...
|
|
34 | 31 |
} |
35 | 32 |
return(list_rast_pred) |
36 | 33 |
} |
... | ... | |
43 | 40 |
list_fitted_models<-vector("list",length(list_formulas)) |
44 | 41 |
for (k in 1:length(list_formulas)){ |
45 | 42 |
formula<-list_formulas[[k]] |
46 |
mod<- try(gam(formula, data=data_training)) |
|
43 |
mod<- try(gam(formula, data=data_training)) #change to any model!! |
|
44 |
#mod<- try(autoKrige(formula, input_data=data_s,new_data=s_sgdf,data_variogram=data_s)) |
|
47 | 45 |
model_name<-paste("mod",k,sep="") |
48 | 46 |
assign(model_name,mod) |
49 | 47 |
list_fitted_models[[k]]<-mod |
... | ... | |
53 | 51 |
#Model and response variable can be changed without affecting the script |
54 | 52 |
prop_month<-0 #proportion retained for validation |
55 | 53 |
run_samp<-1 |
56 |
list_formulas<-vector("list",nmodels) |
|
57 |
|
|
58 |
list_formulas[[1]] <- as.formula("y_var ~ s(elev_1)", env=.GlobalEnv) |
|
59 |
list_formulas[[2]] <- as.formula("y_var ~ s(LST)", env=.GlobalEnv) |
|
60 |
list_formulas[[3]] <- as.formula("y_var ~ s(elev_1,LST)", env=.GlobalEnv) |
|
61 |
list_formulas[[4]] <- as.formula("y_var ~ s(lat) + s(lon)+ s(elev_1)", env=.GlobalEnv) |
|
62 |
list_formulas[[5]] <- as.formula("y_var ~ s(lat,lon,elev_1)", env=.GlobalEnv) |
|
63 |
list_formulas[[6]] <- as.formula("y_var ~ s(lat,lon) + s(elev_1) + s(N_w,E_w) + s(LST)", env=.GlobalEnv) |
|
64 |
list_formulas[[7]] <- as.formula("y_var ~ s(lat,lon) + s(elev_1) + s(N_w,E_w) + s(LST) + s(LC2)", env=.GlobalEnv) |
|
65 |
list_formulas[[8]] <- as.formula("y_var ~ s(lat,lon) + s(elev_1) + s(N_w,E_w) + s(LST) + s(LC6)", env=.GlobalEnv) |
|
66 |
list_formulas[[9]] <- as.formula("y_var ~ s(lat,lon) + s(elev_1) + s(N_w,E_w) + s(LST) + s(DISTOC)", env=.GlobalEnv) |
|
67 |
lst_avg<-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") |
|
68 | 54 |
|
55 |
list_formulas<-lapply(list_models,as.formula,env=.GlobalEnv) #mulitple arguments passed to lapply!! |
|
56 |
|
|
69 | 57 |
data_month<-dst[dst$month==j,] #Subsetting dataset for the relevant month of the date being processed |
70 | 58 |
LST_name<-lst_avg[j] # name of LST month to be matched |
71 | 59 |
data_month$LST<-data_month[[LST_name]] |
... | ... | |
223 | 211 |
names(x)[pos]<-c("id") |
224 | 212 |
names(modst)[1]<-c("id") #modst contains the average tmax per month for every stations... |
225 | 213 |
|
226 |
dmoday <-merge(modst,d,by="id",suffixes=c("",".y2")) #LOOSING DATA HERE!!! from 113 t0 103
|
|
227 |
xmoday <-merge(modst,x,by="id",suffixes=c("",".y2")) #LOOSING DATA HERE!!! from 48 t0 43
|
|
214 |
dmoday <-merge(modst,d,by="id",suffixes=c("",".y2")) |
|
215 |
xmoday <-merge(modst,x,by="id",suffixes=c("",".y2")) |
|
228 | 216 |
mod_pat<-glob2rx("*.y2") |
229 | 217 |
var_pat<-grep(mod_pat,names(dmoday),value=FALSE) # using grep with "value" extracts the matching names |
230 | 218 |
dmoday<-dmoday[,-var_pat] |
... | ... | |
280 | 268 |
|
281 | 269 |
mod_krtmp2<-fitdelta |
282 | 270 |
model_name<-paste("mod_kr","day",sep="_") |
283 |
names(tmp_list)<-names(rast_clim_list) |
|
271 |
names(temp_list)<-names(rast_clim_list)
|
|
284 | 272 |
coordinates(data_s)<-cbind(data_s$x,data_s$y) |
285 | 273 |
proj4string(data_s)<-proj_str |
286 | 274 |
coordinates(data_v)<-cbind(data_v$x,data_v$y) |
Also available in: Unified diff
GAM fusion, VE slight modifications of functions such as model formulas input