Revision b535f939
Added by Benoit Parmentier over 12 years ago
| climate/research/oregon/interpolation/fusion_reg.R | ||
|---|---|---|
| 4 | 4 |
#interpolation area. It requires the text file of stations and a shape file of the study area. # |
| 5 | 5 |
#Note that the projection for both GHCND and study area is lonlat WGS84. # |
| 6 | 6 |
#AUTHOR: Brian McGill # |
| 7 |
#DATE: 06/09/212 #
|
|
| 7 |
#DATE: 06/19/212 #
|
|
| 8 | 8 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#363-- # |
| 9 | 9 |
################################################################################################### |
| 10 | 10 |
|
| ... | ... | |
| 12 | 12 |
library(gtools) # loading some useful tools |
| 13 | 13 |
library(mgcv) # GAM package by Simon Wood |
| 14 | 14 |
library(sp) # Spatial pacakge with class definition by Bivand et al. |
| 15 |
library(spdep) # Spatial pacakge with methods and spatial stat. by Bivand et al. |
|
| 15 |
library(spdep) # Spatial pacakge with methods and spatial stat. by Bivand et al.
|
|
| 16 | 16 |
library(rgdal) # GDAL wrapper for R, spatial utilities |
| 17 | 17 |
library(gstat) # Kriging and co-kriging by Pebesma et al. |
| 18 |
library(fields) # Spatial Interpolation methods such as kriging, splines |
|
| 19 |
library(raster) |
|
| 18 |
library(fields) # NCAR Spatial Interpolation methods such as kriging, splines
|
|
| 19 |
library(raster) # Hijmans et al. package for raster processing
|
|
| 20 | 20 |
### Parameters and argument |
| 21 | 21 |
|
| 22 | 22 |
infile1<- "ghcn_or_tmax_b_04142012_OR83M.shp" #GHCN shapefile containing variables for modeling 2010 |
| ... | ... | |
| 36 | 36 |
|
| 37 | 37 |
prop<-0.3 #Proportion of testing retained for validation |
| 38 | 38 |
seed_number<- 100 #Seed number for random sampling |
| 39 |
out_prefix<-"_06142012_10d_fusion2" #User defined output prefix
|
|
| 39 |
out_prefix<-"_06192012_10d_fusion5" #User defined output prefix
|
|
| 40 | 40 |
setwd(path) |
| 41 | 41 |
############ START OF THE SCRIPT ################## |
| 42 | 42 |
|
| ... | ... | |
| 66 | 66 |
models <-readLines(paste(path,"/",infile4, sep="")) |
| 67 | 67 |
|
| 68 | 68 |
#Model assessment: specific diagnostic/metrics for GAM |
| 69 |
# results_AIC<- matrix(1,length(dates),length(models)+3)
|
|
| 70 |
# results_GCV<- matrix(1,length(dates),length(models)+3)
|
|
| 71 |
# results_DEV<- matrix(1,length(dates),length(models)+3)
|
|
| 72 |
# results_RMSE_f<- matrix(1,length(dates),length(models)+3)
|
|
| 69 |
results_AIC<- matrix(1,length(dates),length(models)+3) |
|
| 70 |
results_GCV<- matrix(1,length(dates),length(models)+3) |
|
| 71 |
results_DEV<- matrix(1,length(dates),length(models)+3) |
|
| 72 |
results_RMSE_f<- matrix(1,length(dates),length(models)+3) |
|
| 73 | 73 |
|
| 74 | 74 |
#Model assessment: general diagnostic/metrics |
| 75 |
results_RMSE <- matrix(1,length(dates),length(models)+3)
|
|
| 76 |
results_MAE <- matrix(1,length(dates),length(models)+3)
|
|
| 77 |
results_ME <- matrix(1,length(dates),length(models)+3)
|
|
| 78 |
results_R2 <- matrix(1,length(dates),length(models)+3) #Coef. of determination for the validation dataset
|
|
| 79 |
results_RMSE_f<- matrix(1,length(dates),length(models)+3) #RMSE fit, RMSE for the training dataset
|
|
| 80 |
results_RMSE_f_kr<- matrix(1,length(dates),length(models)+3)
|
|
| 75 |
results_RMSE <- matrix(1,length(dates),length(models)+4)
|
|
| 76 |
results_MAE <- matrix(1,length(dates),length(models)+4)
|
|
| 77 |
results_ME <- matrix(1,length(dates),length(models)+4) #There are 8+1 models
|
|
| 78 |
results_R2 <- matrix(1,length(dates),length(models)+4) #Coef. of determination for the validation dataset
|
|
| 79 |
results_RMSE_f<- matrix(1,length(dates),length(models)+4) #RMSE fit, RMSE for the training dataset
|
|
| 80 |
results_RMSE_f_kr<- matrix(1,length(dates),length(models)+4)
|
|
| 81 | 81 |
|
| 82 | 82 |
# #Tracking relationship between LST AND LC |
| 83 | 83 |
# cor_LST_LC1<-matrix(1,10,1) #correlation LST-LC1 |
| ... | ... | |
| 132 | 132 |
#year=2010 |
| 133 | 133 |
datelabel=format(ISOdate(year,mo,day),"%b %d, %Y") |
| 134 | 134 |
|
| 135 |
# |
|
| 136 |
# Step 1 - 10 year monthly averages
|
|
| 137 |
# |
|
| 135 |
###########
|
|
| 136 |
# STEP 1 - 10 year monthly averages
|
|
| 137 |
###########
|
|
| 138 | 138 |
|
| 139 |
library(raster) |
|
| 139 |
#library(raster)
|
|
| 140 | 140 |
#old<-getwd() |
| 141 | 141 |
#setwd("c:/data/benoit/data_Oregon_stations_Brian_04242012")
|
| 142 | 142 |
#l=list.files(pattern="mean_month.*rescaled.tif") |
| 143 | 143 |
l=list.files(pattern="mean_month.*rescaled.rst") |
| 144 |
molst<-stack(l) |
|
| 144 |
molst<-stack(l) #Creating a raster stack...
|
|
| 145 | 145 |
#setwd(old) |
| 146 | 146 |
molst=molst-273.16 #K->C |
| 147 | 147 |
idx <- seq(as.Date('2010-01-15'), as.Date('2010-12-15'), 'month')
|
| 148 | 148 |
molst <- setZ(molst, idx) |
| 149 | 149 |
layerNames(molst) <- month.abb |
| 150 |
themolst<-raster(molst,mo) |
|
| 150 |
themolst<-raster(molst,mo) #current month being processed saved in a raster image
|
|
| 151 | 151 |
plot(themolst) |
| 152 |
# |
|
| 153 |
# Step 2 - Weather station means across same days |
|
| 154 |
# |
|
| 152 |
|
|
| 153 |
########### |
|
| 154 |
# STEP 2 - Weather station means across same days: Monthly mean calculation |
|
| 155 |
########### |
|
| 156 |
|
|
| 155 | 157 |
# ??? which years & what quality flags??? |
| 156 | 158 |
#select ghcn.id, lat,lon, elev, month, avg(value/10) as "TMax", count(*) as "NumDays" from ghcn, stations where ghcn.id in (select id from stations where state=='OR') and ghcn.id==stations.id and value<>-9999 and year>=2000 and element=='TMAX' group by stations.id, month;select ghcn.id, lat,lon, elev, month, avg(value/10) as "TMax", count(*) as "NumDays" from ghcn, stations where ghcn.id in (select id from stations where state=='OR') and ghcn.id==stations.id and value<>-9999 and year>=2000 and element=='TMAX' group by stations.id, month; |
| 157 | 159 |
#below table from above SQL query |
| 158 | 160 |
#dst<-read.csv('/data/benoit/data_oregon_stations_brian_04242012/station_means.csv',h=T)
|
| 159 | 161 |
|
| 160 |
|
|
| 161 | 162 |
##Added by Benoit ###### |
| 162 | 163 |
date1<-ISOdate(data3$year,data3$month,data3$day) #Creating a date object from 3 separate column |
| 163 | 164 |
date2<-as.POSIXlt(as.Date(date1)) |
| 164 | 165 |
data3$date<-date2 |
| 165 |
d<-subset(data3,year>=2000 & mflag=="0" ) |
|
| 166 |
d<-subset(data3,year>=2000 & mflag=="0" ) #Selecting dataset 2000-2010 with good quality
|
|
| 166 | 167 |
d1<-aggregate(value~station+month, data=d, mean) #Calculate monthly mean for every station in OR |
| 167 |
id<-as.data.frame(unique(d1$station)) |
|
| 168 |
id<-as.data.frame(unique(d1$station)) #Unique station in OR for year 2000-2010
|
|
| 168 | 169 |
|
| 169 | 170 |
dst<-merge(d1, stat_loc, by.x="station", by.y="STAT_ID") |
| 170 | 171 |
#This allows to change only one name of the data.frame |
| 171 | 172 |
names(dst)[3]<-c("TMax")
|
| 172 |
dst$TMax<-dst$TMax/10 |
|
| 173 |
dst$TMax<-dst$TMax/10 #TMax is the average max temp for months
|
|
| 173 | 174 |
#dstjan=dst[dst$month==9,] #dst contains the monthly averages for tmax for every station over 2000-2010 |
| 174 | 175 |
############## |
| 175 | 176 |
|
| 176 |
modst=dst[dst$month==mo,] |
|
| 177 |
# |
|
| 178 |
# Step 3 - get LST at stations |
|
| 179 |
# |
|
| 180 |
sta_lola=modst[,c("lon","lat")]
|
|
| 177 |
modst=dst[dst$month==mo,] #Subsetting dataset for the relevnat month of the date being processed |
|
| 178 |
|
|
| 179 |
########## |
|
| 180 |
# STEP 3 - get LST at stations |
|
| 181 |
########## |
|
| 182 |
|
|
| 183 |
sta_lola=modst[,c("lon","lat")] #Extracting locations of stations for the current month...
|
|
| 181 | 184 |
library(rgdal) |
| 182 | 185 |
proj_str="+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"; |
| 183 | 186 |
lookup<-function(r,lat,lon) {
|
| ... | ... | |
| 185 | 188 |
cidx<-cellFromXY(r,xy); |
| 186 | 189 |
return(r[cidx]) |
| 187 | 190 |
} |
| 188 |
sta_tmax_from_lst=lookup(themolst,sta_lola$lat,sta_lola$lon) |
|
| 189 |
# |
|
| 190 |
# Step 4 - bias at stations |
|
| 191 |
# |
|
| 192 |
sta_bias=sta_tmax_from_lst-modst$TMax; |
|
| 191 |
sta_tmax_from_lst=lookup(themolst,sta_lola$lat,sta_lola$lon) #Extracted values of LST for the stations |
|
| 192 |
|
|
| 193 |
######### |
|
| 194 |
# STEP 4 - bias at stations |
|
| 195 |
######### |
|
| 196 |
|
|
| 197 |
sta_bias=sta_tmax_from_lst-modst$TMax; #That is the difference between the monthly LST mean and monthly station mean |
|
| 198 |
|
|
| 193 | 199 |
bias_xy=project(as.matrix(sta_lola),proj_str) |
| 194 | 200 |
# windows() |
| 195 |
plot(modst$TMax,sta_tmax_from_lst,xlab="Station mo Tmax",ylab="LST mo Tmax") |
|
| 201 |
X11() |
|
| 202 |
plot(modst$TMax,sta_tmax_from_lst,xlab="Station mo Tmax",ylab="LST mo Tmax",main=paste("LST vs TMax for",datelabel,sep=" "))
|
|
| 196 | 203 |
abline(0,1) |
| 197 |
# |
|
| 198 |
# Step 5 - interpolate bias |
|
| 199 |
# |
|
| 204 |
savePlot(paste("LST_TMax_scatterplot_",dates[i],out_prefix,".png", sep=""), type="png")
|
|
| 205 |
dev.off() |
|
| 206 |
|
|
| 207 |
######## |
|
| 208 |
# STEP 5 - interpolate bias |
|
| 209 |
######## |
|
| 210 |
|
|
| 200 | 211 |
# ?? include covariates like elev, distance to coast, cloud frequency, tree height |
| 201 |
library(fields) |
|
| 212 |
#library(fields)
|
|
| 202 | 213 |
#windows() |
| 203 | 214 |
quilt.plot(sta_lola,sta_bias,main="Bias at stations",asp=1) |
| 204 | 215 |
US(add=T,col="magenta",lwd=2) |
| 205 | 216 |
#fitbias<-Tps(bias_xy,sta_bias) #use TPS or krige |
| 206 |
fitbias<-Krig(bias_xy,sta_bias,theta=1e5) #use TPS or krige |
|
| 207 |
# windows() |
|
| 208 |
surface(fitbias,col=rev(terrain.colors(100)),asp=1,main="Interpolated bias") |
|
| 217 |
fitbias<-Krig(bias_xy,sta_bias,theta=1e5) #use TPS or krige |
|
| 218 |
#The output is a krig object using fields |
|
| 219 |
# Creating plot of bias surface and saving it |
|
| 220 |
X11() |
|
| 221 |
datelabel2=format(ISOdate(year,mo,day),"%B ") #added by Benoit, label |
|
| 222 |
surface(fitbias,col=rev(terrain.colors(100)),asp=1,main=paste("Interpolated bias for",datelabel2,sep=" "))
|
|
| 223 |
savePlot(paste("Bias_surface_LST_TMax_",dates[i],out_prefix,".png", sep=""), type="png")
|
|
| 224 |
dev.off() |
|
| 225 |
|
|
| 209 | 226 |
#US(add=T,col="magenta",lwd=2) |
| 210 |
# |
|
| 211 |
# Step 6 - return to daily station data & calcualate delta=daily T-monthly T from stations |
|
| 227 |
|
|
| 228 |
########## |
|
| 229 |
# STEP 6 - return to daily station data & calcualate delta=daily T-monthly T from stations |
|
| 230 |
########## |
|
| 231 |
|
|
| 212 | 232 |
#Commmented out by Benoit 06/14 |
| 213 | 233 |
# library(RSQLite) |
| 214 | 234 |
# m<-dbDriver("SQLite")
|
| ... | ... | |
| 236 | 256 |
names(dmoday)[4]<-c("lat")
|
| 237 | 257 |
names(dmoday)[5]<-c("lon")
|
| 238 | 258 |
### |
| 259 |
#dmoday contains the daily tmax values with TMax being the monthly station tmax mean |
|
| 239 | 260 |
|
| 240 | 261 |
# windows() |
| 262 |
X11() |
|
| 241 | 263 |
plot(dailyTmax~TMax,data=dmoday,xlab="Mo Tmax",ylab=paste("Daily for",datelabel),main="across stations in OR")
|
| 242 |
# |
|
| 243 |
# Step 7 - interpolate delta across space |
|
| 244 |
# |
|
| 264 |
savePlot(paste("Daily_tmax_monthly_TMax_scatterplot_",dates[i],out_prefix,".png", sep=""), type="png")
|
|
| 265 |
dev.off() |
|
| 266 |
########## |
|
| 267 |
# STEP 7 - interpolate delta across space |
|
| 268 |
########## |
|
| 269 |
|
|
| 245 | 270 |
daily_sta_lola=dmoday[,c("lon","lat")] #could be same as before but why assume merge does this - assume not
|
| 246 | 271 |
daily_sta_xy=project(as.matrix(daily_sta_lola),proj_str) |
| 247 | 272 |
daily_delta=dmoday$dailyTmax-dmoday$TMax |
| ... | ... | |
| 250 | 275 |
US(add=T,col="magenta",lwd=2) |
| 251 | 276 |
#fitdelta<-Tps(daily_sta_xy,daily_delta) #use TPS or krige |
| 252 | 277 |
fitdelta<-Krig(daily_sta_xy,daily_delta,theta=1e5) #use TPS or krige |
| 253 |
# windows() |
|
| 254 |
surface(fitdelta,col=rev(terrain.colors(100)),asp=1,main="Interpolated delta") |
|
| 278 |
#Kriging using fields package |
|
| 279 |
|
|
| 280 |
# Creating plot of bias surface and saving it |
|
| 281 |
X11() |
|
| 282 |
surface(fitdelta,col=rev(terrain.colors(100)),asp=1,main=paste("Interpolated delta for",datelabel,sep=" "))
|
|
| 283 |
savePlot(paste("Delta_surface_LST_TMax_",dates[i],out_prefix,".png", sep=""), type="png")
|
|
| 284 |
dev.off() |
|
| 255 | 285 |
#US(add=T,col="magenta",lwd=2) |
| 256 | 286 |
# |
| 257 |
# Step 8 - assemble final answer - T=LST+Bias(interpolated)+delta(interpolated) |
|
| 258 |
# |
|
| 259 |
bias_rast=interpolate(themolst,fitbias) |
|
| 260 |
plot(bias_rast,main="Raster bias") |
|
| 261 |
daily_delta_rast=interpolate(themolst,fitdelta) |
|
| 287 |
|
|
| 288 |
#### Added by Benoit on 06/19 |
|
| 289 |
data_s<-dmoday #put the |
|
| 290 |
data_s$daily_delta<-daily_delta |
|
| 291 |
data_s$y_var<-daily_delta #y_var is the variable currently being modeled, may be better with BIAS!! |
|
| 292 |
|
|
| 293 |
#Model and response variable can be changed without affecting the script |
|
| 294 |
|
|
| 295 |
mod1<- gam(y_var~ s(lat) + s (lon) + s (ELEV_SRTM), data=data_s) |
|
| 296 |
mod2<- gam(y_var~ s(lat,lon)+ s(ELEV_SRTM), data=data_s) #modified nesting....from 3 to 2 |
|
| 297 |
mod3<- gam(y_var~ s(lat) + s (lon) + s (ELEV_SRTM) + s (Northness)+ s (Eastness) + s(DISTOC), data=data_s) |
|
| 298 |
mod4<- gam(y_var~ s(lat) + s (lon) + s(ELEV_SRTM) + s(Northness) + s (Eastness) + s(DISTOC) + s(LST), data=data_s) |
|
| 299 |
mod5<- gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST), data=data_s) |
|
| 300 |
mod6<- gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC1), data=data_s) |
|
| 301 |
mod7<- gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC3), data=data_s) |
|
| 302 |
mod8<- gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST) + s(LC1), data=data_s) |
|
| 303 |
#### Added by Benoit ends |
|
| 304 |
|
|
| 305 |
######### |
|
| 306 |
# STEP 8 - assemble final answer - T=LST+Bias(interpolated)+delta(interpolated) |
|
| 307 |
######### |
|
| 308 |
|
|
| 309 |
bias_rast=interpolate(themolst,fitbias) #interpolation using function from raster package |
|
| 310 |
#themolst is raster layer, fitbias is "Krig" object from bias surface |
|
| 311 |
plot(bias_rast,main="Raster bias") #This not displaying... |
|
| 312 |
|
|
| 313 |
daily_delta_rast=interpolate(themolst,fitdelta) #Interpolation of the bias surface... |
|
| 314 |
|
|
| 262 | 315 |
plot(daily_delta_rast,main="Raster Daily Delta") |
| 263 |
tmax_predicted=themolst+daily_delta_rast-bias_rast #Final surface?? |
|
| 316 |
|
|
| 317 |
tmax_predicted=themolst+daily_delta_rast-bias_rast #Final surface?? but daily_rst |
|
| 318 |
#tmax_predicted=themolst+daily_delta_rast+bias_rast #Added by Benoit, why is it -bias_rast |
|
| 264 | 319 |
plot(tmax_predicted,main="Predicted daily") |
| 265 |
# |
|
| 266 |
# check |
|
| 267 |
# |
|
| 320 |
|
|
| 321 |
######## |
|
| 322 |
# check: assessment of results: validation |
|
| 323 |
######## |
|
| 324 |
RMSE<-function(x,y) {return(mean((x-y)^2)^0.5)}
|
|
| 325 |
|
|
| 326 |
#FIT ASSESSMENT |
|
| 327 |
sta_pred_data_s=lookup(tmax_predicted,data_s$lat,data_s$lon) |
|
| 328 |
rmse_fit=RMSE(sta_pred_data_s,data_s$dailyTmax) |
|
| 329 |
|
|
| 268 | 330 |
sta_pred=lookup(tmax_predicted,data_v$lat,data_v$lon) |
| 269 | 331 |
#sta_pred=lookup(tmax_predicted,daily_sta_lola$lat,daily_sta_lola$lon) |
| 270 |
RMSE<-function(x,y) {return(mean((x-y)^2)^0.5)}
|
|
| 271 |
rmse=RMSE(sta_pred,dmoday$dailyTmax) |
|
| 272 |
#plot(sta_pred~dmoday$dailyTmax,xlab=paste("Actual daily for",datelabel),ylab="Pred daily",main=paste("RMSE=",rmse))
|
|
| 332 |
#rmse=RMSE(sta_pred,dmoday$dailyTmax) |
|
| 273 | 333 |
tmax<-data_v$tmax/10 |
| 334 |
rmse=RMSE(sta_pred,tmax) |
|
| 335 |
#plot(sta_pred~dmoday$dailyTmax,xlab=paste("Actual daily for",datelabel),ylab="Pred daily",main=paste("RMSE=",rmse))
|
|
| 336 |
X11() |
|
| 274 | 337 |
plot(sta_pred~tmax,xlab=paste("Actual daily for",datelabel),ylab="Pred daily",main=paste("RMSE=",rmse))
|
| 275 | 338 |
abline(0,1) |
| 276 |
resid=sta_pred-dmoday$dailyTmax |
|
| 339 |
savePlot(paste("Predicted_tmax_versus_observed_scatterplot_",dates[i],out_prefix,".png", sep=""), type="png")
|
|
| 340 |
dev.off() |
|
| 341 |
#resid=sta_pred-dmoday$dailyTmax |
|
| 342 |
resid=sta_pred-tmax |
|
| 277 | 343 |
quilt.plot(daily_sta_lola,resid) |
| 278 | 344 |
|
| 279 | 345 |
### END OF BRIAN's code |
| 280 | 346 |
|
| 281 | 347 |
### Added by benoit |
| 282 |
j=1 |
|
| 348 |
#Store results using TPS |
|
| 349 |
j=9 |
|
| 283 | 350 |
results_RMSE[i,1]<- dates[i] #storing the interpolation dates in the first column |
| 284 | 351 |
results_RMSE[i,2]<- ns #number of stations used in the training stage |
| 285 | 352 |
results_RMSE[i,3]<- "RMSE" |
| 286 | 353 |
results_RMSE[i,j+3]<- rmse #Storing RMSE for the model j |
| 287 | 354 |
|
| 355 |
results_RMSE_f[i,1]<- dates[i] #storing the interpolation dates in the first column |
|
| 356 |
results_RMSE_f[i,2]<- ns #number of stations used in the training stage |
|
| 357 |
results_RMSE_f[i,3]<- "RMSE" |
|
| 358 |
results_RMSE_f[i,j+3]<- rmse_fit #Storing RMSE for the model j |
|
| 359 |
|
|
| 360 |
ns<-nrow(data_s) |
|
| 361 |
|
|
| 362 |
|
|
| 363 |
for (j in 1:length(models)){
|
|
| 364 |
|
|
| 365 |
##Model assessment: specific diagnostic/metrics for GAM |
|
| 366 |
|
|
| 367 |
name<-paste("mod",j,sep="") #modj is the name of The "j" model (mod1 if j=1)
|
|
| 368 |
mod<-get(name) #accessing GAM model ojbect "j" |
|
| 369 |
results_AIC[i,1]<- dates[i] #storing the interpolation dates in the first column |
|
| 370 |
results_AIC[i,2]<- ns #number of stations used in the training stage |
|
| 371 |
results_AIC[i,3]<- "AIC" |
|
| 372 |
results_AIC[i,j+3]<- AIC (mod) |
|
| 373 |
|
|
| 374 |
results_GCV[i,1]<- dates[i] #storing the interpolation dates in the first column |
|
| 375 |
results_GCV[i,2]<- ns #number of stations used in the training |
|
| 376 |
results_GCV[i,3]<- "GCV" |
|
| 377 |
results_GCV[i,j+3]<- mod$gcv.ubre |
|
| 378 |
|
|
| 379 |
results_DEV[i,1]<- dates[i] #storing the interpolation dates in the first column |
|
| 380 |
results_DEV[i,2]<- ns #number of stations used in the training stage |
|
| 381 |
results_DEV[i,3]<- "DEV" |
|
| 382 |
results_DEV[i,j+3]<- mod$deviance |
|
| 383 |
|
|
| 384 |
results_RMSE_f[i,1]<- dates[i] #storing the interpolation dates in the first column |
|
| 385 |
results_RMSE_f[i,2]<- ns #number of stations used in the training stage |
|
| 386 |
results_RMSE_f[i,3]<- "RSME" |
|
| 387 |
results_RMSE_f[i,j+3]<- sqrt(sum((mod$residuals)^2)/nv) |
|
| 388 |
|
|
| 389 |
##Model assessment: general diagnostic/metrics |
|
| 390 |
##validation: using the testing data |
|
| 391 |
|
|
| 392 |
#This was modified on 06192012 |
|
| 393 |
y_mod<- predict(mod, newdata=data_v, se.fit = TRUE) #Using the coeff to predict new values. |
|
| 394 |
sta_LST=lookup(themolst,data_v$lat,data_v$lon) |
|
| 395 |
sta_bias=lookup(bias_rast,data_v$lat,data_v$lon) |
|
| 396 |
tmax_predicted=sta_LST+sta_bias-y_mod$fit |
|
| 397 |
|
|
| 398 |
data_v$tmax<-(data_v$tmax)/10 |
|
| 399 |
res_mod<- data_v$tmax - tmax_predicted #Residuals for the model |
|
| 400 |
#res_mod<- data_v$tmax - y_mod$fit #Residuals for the model |
|
| 401 |
|
|
| 402 |
RMSE_mod <- sqrt(sum(res_mod^2)/nv) #RMSE FOR REGRESSION STEP 1: GAM |
|
| 403 |
MAE_mod<- sum(abs(res_mod))/nv #MAE, Mean abs. Error FOR REGRESSION STEP 1: GAM |
|
| 404 |
ME_mod<- sum(res_mod)/nv #ME, Mean Error or bias FOR REGRESSION STEP 1: GAM |
|
| 405 |
R2_mod<- cor(data_v$tmax,y_mod$fit)^2 #R2, coef. of var FOR REGRESSION STEP 1: GAM |
|
| 406 |
|
|
| 407 |
results_RMSE[i,1]<- dates[i] #storing the interpolation dates in the first column |
|
| 408 |
results_RMSE[i,2]<- ns #number of stations used in the training stage |
|
| 409 |
results_RMSE[i,3]<- "RMSE" |
|
| 410 |
results_RMSE[i,j+3]<- RMSE_mod #Storing RMSE for the model j |
|
| 411 |
results_MAE[i,1]<- dates[i] #storing the interpolation dates in the first column |
|
| 412 |
results_MAE[i,2]<- ns #number of stations used in the training stage |
|
| 413 |
results_MAE[i,3]<- "MAE" |
|
| 414 |
results_MAE[i,j+3]<- MAE_mod #Storing MAE for the model j |
|
| 415 |
results_ME[i,1]<- dates[i] #storing the interpolation dates in the first column |
|
| 416 |
results_ME[i,2]<- ns #number of stations used in the training stage |
|
| 417 |
results_ME[i,3]<- "ME" |
|
| 418 |
results_ME[i,j+3]<- ME_mod #Storing ME for the model j |
|
| 419 |
results_R2[i,1]<- dates[i] #storing the interpolation dates in the first column |
|
| 420 |
results_R2[i,2]<- ns #number of stations used in the training stage |
|
| 421 |
results_R2[i,3]<- "R2" |
|
| 422 |
results_R2[i,j+3]<- R2_mod #Storing R2 for the model j |
|
| 423 |
|
|
| 424 |
#Saving residuals and prediction in the dataframes: tmax predicted from GAM |
|
| 425 |
pred<-paste("pred_mod",j,sep="")
|
|
| 426 |
data_v[[pred]]<-as.numeric(y_mod$fit) |
|
| 427 |
data_s[[pred]]<-as.numeric(mod$fit) #Storing model fit values (predicted on training sample) |
|
| 428 |
|
|
| 429 |
name2<-paste("res_mod",j,sep="")
|
|
| 430 |
data_v[[name2]]<-as.numeric(res_mod) |
|
| 431 |
data_s[[name2]]<-as.numeric(mod$residuals) |
|
| 432 |
#end of loop calculating RMSE |
|
| 433 |
|
|
| 434 |
} |
|
| 435 |
|
|
| 288 | 436 |
# end of the for loop1 |
| 289 | 437 |
|
| 290 | 438 |
} |
| ... | ... | |
| 299 | 447 |
results_table_ME<-as.data.frame(results_ME) |
| 300 | 448 |
results_table_R2<-as.data.frame(results_R2) |
| 301 | 449 |
|
| 302 |
cname<-c("dates","ns","metric","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7","mod8")
|
|
| 450 |
results_table_RMSE_f<-as.data.frame(results_RMSE_f) |
|
| 451 |
|
|
| 452 |
cname<-c("dates","ns","metric","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7","mod8","mod9")
|
|
| 303 | 453 |
colnames(results_table_RMSE)<-cname |
| 454 |
colnames(results_table_RMSE_f)<-cname |
|
| 304 | 455 |
|
| 305 | 456 |
#tb_diagnostic1<-rbind(results_table_RMSE,results_table_MAE, results_table_ME, results_table_R2) # |
| 306 |
tb_diagnostic1<-results_table_RMSE |
|
| 457 |
tb_diagnostic1<-results_table_RMSE #measures of validation |
|
| 458 |
tb_diagnostic2<-results_table_RMSE_f #measures of fit |
|
| 307 | 459 |
|
| 308 | 460 |
write.table(tb_diagnostic1, file= paste(path,"/","results_fusion_Assessment_measure1",out_prefix,".txt",sep=""), sep=",") |
| 461 |
write.table(tb_diagnostic2, file= paste(path,"/","results_fusion_Assessment_measure2",out_prefix,".txt",sep=""), sep=",") |
|
| 309 | 462 |
|
| 310 | 463 |
#### END OF SCRIPT |
Also available in: Unified diff
FUSION, modifications saving plots and adding GAM models for delta surface