Revision 75151566
Added by Benoit Parmentier about 12 years ago
| climate/research/oregon/interpolation/GAM_LST_reg.R | ||
|---|---|---|
| 1 | 1 |
####################Interpolation of Tmax for 10 dates.##################### |
| 2 |
#This script interpolates station values for the Oregon case study. This program loads the station data from a csv file
|
|
| 3 |
#and perform two types of regression: multiple linear model and general additive model (GAM). Note that this program:
|
|
| 4 |
#1)assumes that the csv file is in the current working
|
|
| 2 |
#This script interpolates station values for the Oregon case study. This program loads the station data from a shapefile
|
|
| 3 |
#and perform 8 regressions using the general additive model (GAM). Note that this program:
|
|
| 4 |
#1)assumes that the shapefile in the current working
|
|
| 5 | 5 |
#2)extract relevant variables from raster images before performing the regressions. |
| 6 | 6 |
#This scripts predicts tmas xsing GAM and LST derived from MOD11A1. |
| 7 | 7 |
#Interactions terms are also included and assessed using the RMSE from validation dataset. |
| ... | ... | |
| 202 | 202 |
|
| 203 | 203 |
results_table_RMSE<-as.data.frame(results_RMSEnum) |
| 204 | 204 |
results_table_AIC<-as.data.frame(results_AICnum) |
| 205 |
colnames(results_table_RMSE)<-c("dates","ns","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7")
|
|
| 206 |
colnames(results_table_AIC)<-c("dates","ns","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7")
|
|
| 205 |
colnames(results_table_RMSE)<-c("dates","ns","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7", "mod8")
|
|
| 206 |
colnames(results_table_AIC)<-c("dates","ns","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7", "mod8")
|
|
| 207 | 207 |
|
| 208 | 208 |
#results_table_RMSE |
| 209 |
#write.csv(results_table_RMSE, file= paste(path,"/","results_GAM_Assessment",out_prefix,".txt",sep="")) |
|
| 210 |
#write.csv(results_table_AIC, file= paste(path,"/","results_GAM_Assessment",out_prefix,".txt",sep=""),append=TRUE) |
|
| 211 |
|
|
| 212 | 209 |
write.table(results_table_RMSE, file= paste(path,"/","results_GAM_Assessment",out_prefix,".txt",sep=""), sep=",") |
| 213 | 210 |
write.table(results_table_AIC, file= paste(path,"/","results_GAM_Assessment",out_prefix,".txt",sep=""),sep=",", append=TRUE) |
| 214 | 211 |
|
| 215 |
#Can also use file connection |
|
| 216 |
f<-file(paste(path,"/","results_GAM_Assessment",out_prefix,"2.txt",sep=""),"w") |
|
| 217 |
write(results_table_RMSE, file=f) |
|
| 218 |
close(f) |
|
| 219 |
# End of script########## |
|
| 212 |
###Analysing the results from the 365 days run: Summarize by month |
|
| 213 |
|
|
| 214 |
for(i in 1:nrow(results_table_RMSE)){
|
|
| 215 |
date<-results_table_RMSE$dates[i] |
|
| 216 |
date<-strptime(date, "%Y%m%d") |
|
| 217 |
results_table_RMSE$month[i]<-as.integer(strftime(date, "%m")) |
|
| 218 |
} |
|
| 219 |
|
|
| 220 |
average<-aggregate(cbind(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8)~month,data=results_table_RMSE,mean, na.rm=TRUE) |
|
| 221 |
average<-aggregate(cbind(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8)~month,data=results_table_RMSE, FUN=mean) |
|
| 222 |
#average on all the data.frame |
|
| 223 |
averaget<-aggregate(results_table_RMSE, by=list(results_table_RMSE$month),FUN=mean, na.rm=TRUE) |
|
| 224 |
#mediant<-aggregate(results_table_RMSE, by=list(results_table_RMSE$month),FUN=median, na.rm=TRUE) |
|
| 225 |
#average_lowt<-aggregate(results_table_RMSE, by=list(results_table_RMSE$month), FUN=function(v) t.test(v)$conf.int[1]) |
|
| 226 |
#average_up<-aggregate(cbind(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8)~month,data=results_table_RMSE, function(v) t.test(v)$conf.int[2]) |
|
| 227 |
|
|
| 228 |
median<-aggregate(cbind(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8)~month,data=results_table_RMSE, median, na.rm=TRUE) |
|
| 229 |
average_low<-aggregate(cbind(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8)~month,data=results_table_RMSE, function(v) t.test(v)$conf.int[1]) |
|
| 230 |
average_up<-aggregate(cbind(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8)~month,data=results_table_RMSE, function(v) t.test(v)$conf.int[2]) |
|
| 220 | 231 |
|
| 221 |
#Selecting dates and files based on names |
|
| 222 |
#cor_LST_LC<-matrix(1,10,1) |
|
| 223 |
# for(i in 1:length(dates)){
|
|
| 224 |
# cor_LST_LC1[i]<-cor(ghcn.subsets[[i]]$LST,ghcn.subsets[[i]]$LC1) |
|
| 225 |
# } |
|
| 226 |
# for(i in 1:length(dates)){
|
|
| 227 |
# cor_LST_LC3[i]<-cor(ghcn.subsets[[i]]$LST,ghcn.subsets[[i]]$LC3) |
|
| 228 |
# } |
|
| 232 |
mod<-names(averaget) |
|
| 233 |
mod<-mod[4:11] |
|
| 234 |
#Saving graphic plots |
|
| 235 |
for(i in 1:length(mod)){
|
|
| 236 |
X11(width=14,height=10) |
|
| 237 |
name<-mod[i] |
|
| 238 |
barplot2(average[[name]],plot.ci=TRUE, ci.l=average_low[[name]], ci.u=average_up[[name]],main="Mean RMSE per month", names.arg=c("Jan", "Feb", "Mar", "Apr", "May", "Jun","Jul", "Aug", "Sep","Oct", "Nov", "Dec"),ylim=c(20,30),ylab="RMSE in tenth deg C",xlab=name)
|
|
| 239 |
#title(paste("Sampling RMSE for mod",i,sep=""))
|
|
| 240 |
savePlot(paste("barplot_results_RMSE_month_",name,out_prefix,".png", sep=""), type="png")
|
|
| 241 |
dev.off() |
|
| 242 |
} |
|
| 243 |
|
|
| 244 |
|
|
| 245 |
for(i in 1:length(mod)){
|
|
| 246 |
X11(width=14,height=10) |
|
| 247 |
name<-mod[i] |
|
| 248 |
barplot2(average[[name]],plot.ci=TRUE, ci.l=average_low[[name]], ci.u=average_up[[name]],main=paste(" Mean RMSE per month ",name, sep=""), names.arg=c("Jan", "Feb", "Mar", "Apr", "May", "Jun","Jul", "Aug", "Sep","Oct", "Nov", "Dec"),ylim=c(20,30),ylab="RMSE in tenth deg C",xlab=name)
|
|
| 249 |
#title(paste("Sampling RMSE for mod",i,sep=""))
|
|
| 250 |
savePlot(paste("barplot_results_RMSE_month_",name,out_prefix,".png", sep=""), type="png")
|
|
| 251 |
dev.off() |
|
| 252 |
|
|
| 253 |
X11(width=14,height=10) |
|
| 254 |
name<-mod[i] |
|
| 255 |
hist(results_table_RMSE[[name]],breaks=15, main=paste(" Histogram RMSE_",name, sep=""),xlab=paste("RMSE ",name, sep=""))
|
|
| 256 |
savePlot(paste("Hist_results_RMSE_365_",name,out_prefix,".png", sep=""), type="png")
|
|
| 257 |
dev.off() |
|
| 258 |
|
|
| 259 |
} |
|
| 260 |
|
|
| 261 |
for(i in 1:length(mod)){
|
|
| 262 |
X11(width=14,height=10) |
|
| 263 |
name<-mod[i] |
|
| 264 |
hist(results_table_RMSE[[name]],breaks=15, main=paste(" Histogram RMSE_",name, sep=""),xlab=paste("RMSE ",name, sep=""))
|
|
| 265 |
savePlot(paste("Hist_results_RMSE_365_",name,out_prefix,".png", sep=""), type="png")
|
|
| 266 |
dev.off() |
|
| 267 |
} |
|
| 268 |
|
|
| 269 |
r<-(results_table_RMSE[,3:10]) #selecting only the columns related to models... |
|
| 270 |
|
|
| 271 |
mean_r<-mean(r) |
|
| 272 |
median_r<-sapply(r, median) |
|
| 273 |
sd_r<-sapply(r, sd) |
|
| 274 |
|
|
| 275 |
barplot(mean_r,ylim=c(23,26),ylab="RMSE in tenth deg C") |
|
| 276 |
barplot(median_r,ylim=c(23,26),ylab="RMSE in tenth deg C",add=TRUE,inside=FALSE,beside=TRUE) # put both on the same plot |
|
| 277 |
barplot(sd_r,ylim=c(6,8),ylab="RMSE in tenth deg C") # put both on the same plot |
|
| 278 |
|
|
| 279 |
height<-rbind(mean_r,median_r) |
|
| 280 |
barplot(height,ylim=c(23,26),ylab="RMSE in tenth deg C",beside=TRUE,legend=rownames(height)) |
|
| 281 |
barplot(height,ylim=c(23,26),ylab="RMSE in tenth deg C",beside=TRUE, col=c("darkblue","red"),legend=rownames(height)) # put both on the same plot
|
|
| 282 |
|
|
| 283 |
barplot2(mean_r,median_r,ylim=c(23,26),ylab="RMSE in tenth deg C") # put both on the same plot |
|
| 284 |
#Collect var explained and p values for each var... |
|
| 285 |
|
|
| 286 |
# End of script########## |
|
| 229 | 287 |
|
| 230 |
#mod9<- gam(tmax~ te(lat,lon,ELEV_SRTM), data=ghcn_s_20101031) |
|
| 288 |
names(results_table_RMSE) |
|
Also available in: Unified diff
GAM LST, adding lines to assess RMSE per month over 365 dates, task #406