Revision c28e9f58
Added by Benoit Parmentier about 10 years ago
| climate/research/oregon/interpolation/multi_timescale_paper_interpolation.R | ||
|---|---|---|
| 7 | 7 |
#Analyses, figures, tables and data for the paper are also produced in the script. |
| 8 | 8 |
#AUTHOR: Benoit Parmentier |
| 9 | 9 |
#CREATED ON: 10/31/2013 |
| 10 |
#MODIFIED ON: 08/12/2014
|
|
| 10 |
#MODIFIED ON: 10/06/2014
|
|
| 11 | 11 |
#Version: 6 |
| 12 | 12 |
#PROJECT: Environmental Layers project |
| 13 | 13 |
################################################################################################# |
| ... | ... | |
| 42 | 42 |
#### FUNCTION USED IN SCRIPT |
| 43 | 43 |
|
| 44 | 44 |
function_analyses_paper1 <-"contribution_of_covariates_paper_interpolation_functions_07182014.R" #first interp paper |
| 45 |
function_analyses_paper2 <-"multi_timescales_paper_interpolation_functions_08132014.R"
|
|
| 45 |
function_analyses_paper2 <-"multi_timescales_paper_interpolation_functions_10062014.R"
|
|
| 46 | 46 |
|
| 47 | 47 |
############################## |
| 48 | 48 |
#### Parameters and constants |
| ... | ... | |
| 111 | 111 |
"gam_fss","kriging_fss","gwr_fss") |
| 112 | 112 |
|
| 113 | 113 |
y_var_name <- "dailyTmax" |
| 114 |
out_prefix<-"analyses_08112014"
|
|
| 114 |
out_prefix<-"analyses_10062014"
|
|
| 115 | 115 |
out_dir<-"/home/parmentier/Data/IPLANT_project/paper_multitime_scale__analyses_tables" |
| 116 | 116 |
create_out_dir_param = TRUE |
| 117 | 117 |
|
| ... | ... | |
| 127 | 127 |
infile_reg_outline <- "/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/OR83M_state_outline.shp" #input region outline defined by polygon: Oregon |
| 128 | 128 |
met_stations_outfiles_obj_file<-"/data/project/layers/commons/data_workflow/output_data_365d_gam_fus_lst_test_run_07172013/met_stations_outfiles_obj_gam_fusion__365d_gam_fus_lst_test_run_07172013.RData" |
| 129 | 129 |
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| 130 |
CRS_WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
|
| 131 |
|
|
| 130 | 132 |
ref_rast_name<- "/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/mean_day244_rescaled.rst" #This is the shape file of outline of the study area. #local raster name defining resolution, exent, local projection--. set on the fly?? |
| 131 | 133 |
infile_reg_outline <- "/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/OR83M_state_outline.shp" #input region outline defined by polygon: Oregon |
| 132 | 134 |
ref_rast_name <-"/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/mean_day244_rescaled.rst" #local raster name defining resolution, exent: oregon |
| ... | ... | |
| 370 | 372 |
file_name<-paste("table6_correlation_multi_timescale_paper","_",out_prefix,".txt",sep="")
|
| 371 | 373 |
write.table(table6,file=file_name,sep=",") |
| 372 | 374 |
|
| 373 |
### Additiona information: |
|
| 375 |
### Additional information:
|
|
| 374 | 376 |
|
| 375 | 377 |
#model with LST and elev |
| 376 | 378 |
raster_obj <- load_obj(list_raster_obj_files$gam_CAI) |
| ... | ... | |
| 402 | 404 |
list_tb_mod_month[[i]] <- test |
| 403 | 405 |
|
| 404 | 406 |
} |
| 405 |
|
|
| 407 |
names(list_tb_mod_month) <- interp_method_selected |
|
| 408 |
list_formulas |
|
| 406 | 409 |
list_tb_mod_month[[1]][,c(3,6,7)] |
| 407 | 410 |
list_tb_mod_month[[2]][,c(3,6,7)] |
| 408 | 411 |
list_tb_mod_month[[3]][,c(3,6,7)] |
| 409 | 412 |
#xyplot(test) |
| 410 | 413 |
|
| 411 |
list_tb_mod_month[[1]][,c(6,7)]-list_tb_mod_month[[1]][,c(3)] |
|
| 412 |
list_tb_mod_month[[2]][,c(6,7)]-list_tb_mod_month[[1]][,c(3)] |
|
| 413 |
list_tb_mod_month[[3]][,c(6,7)]-list_tb_mod_month[[1]][,c(3)] |
|
| 414 |
|
|
| 415 |
#model with lev only |
|
| 416 |
tb1_mod1_month <- raster_prediction_obj_1$summary_month_metrics_v[[1]] #note that this is for model1 |
|
| 417 |
|
|
| 418 |
tb1_mod1_month<-raster_prediction_obj_1$summary_month_metrics_v[[1]] #note that this is for model1 |
|
| 419 |
rmse_dif <- (tb1_mod4_month$rmse) - tb1_mod1_month$rmse |
|
| 420 |
plot(rmse_dif) |
|
| 414 |
rmse_dif_gam_CAI <- list_tb_mod_month[[1]][,c(6,7)]-list_tb_mod_month[[1]][,c(3)] |
|
| 415 |
rmse_dif_kriging_CAI <- list_tb_mod_month[[2]][,c(6,7)]-list_tb_mod_month[[2]][,c(3)] |
|
| 416 |
rmse_dif_gwr_CAI <- list_tb_mod_month[[3]][,c(6,7)]-list_tb_mod_month[[3]][,c(3)] |
|
| 417 |
|
|
| 418 |
rmse_dif_gam_CAI$interp_method <- interp_method_selected[1] |
|
| 419 |
rmse_dif_kriging_CAI$interp_method <- interp_method_selected[2] |
|
| 420 |
rmse_dif_gwr_CAI$interp_method <- interp_method_selected[3] |
|
| 421 |
|
|
| 422 |
mydata <- do.call(rbind,list(rmse_dif_gam_CAI,rmse_dif_kriging_CAI,rmse_dif_gwr_CAI)) |
|
| 423 |
dd <- do.call(make.groups, mydata[,-ncol(mydata)]) |
|
| 424 |
dd$interp_method <- mydata$interp_method |
|
| 425 |
dd$month <- 1:12 |
|
| 426 |
|
|
| 427 |
layout_m<-c(1,3) #one row two columns |
|
| 428 |
|
|
| 429 |
png(paste("Figure_11_improvement_dure to LST and models_tmax_prediction_",out_prefix,".png", sep=""),
|
|
| 430 |
height=500*layout_m[1],width=500*layout_m[2]) |
|
| 431 |
|
|
| 432 |
p<-xyplot(data ~ month | interp_method, groups=which,dd,type="b", |
|
| 433 |
par.settings = list(axis.text = list(font = 2, cex = 1.3),layout=layout_m, |
|
| 434 |
par.main.text=list(font=2,cex=2),strip.background=list(col="white")),par.strip.text=list(font=2,cex=1.5), |
|
| 435 |
#strip=strip.custom(factor.levels=names_layers), |
|
| 436 |
xlab=list(label="Month", cex=2,font=2), |
|
| 437 |
ylab=list(label="difference to mod3", cex=2, font=2), |
|
| 438 |
auto.key=list(columns=2,space="top",cex=2.5,font=2), |
|
| 439 |
|
|
| 440 |
#list(label="\u0394RMSE between mod1 and mod4",cex=1.5) |
|
| 441 |
as.table=TRUE) #as.table controls the order of the pannels!! |
|
| 442 |
p <- update(p, panel = function(...) {
|
|
| 443 |
panel.abline(h = 0, lty = 2, col = "gray") |
|
| 444 |
panel.xyplot(...) |
|
| 445 |
}) |
|
| 446 |
print(p) |
|
| 447 |
|
|
| 448 |
#p_dif <- xyplot(rmse_dif ~ 1:12, |
|
| 449 |
# col=c("black"),
|
|
| 450 |
# type="b", |
|
| 451 |
# ylab=list(label="\u0394RMSE between mod1 and mod4",cex=1.5), |
|
| 452 |
# xlab=list(label="Month",cex=1.5), |
|
| 453 |
# main=list(label="(b) Monthly \u0394RMSE betwen mod1 and mod4",cex=1.8), |
|
| 454 |
# par.settings = list(axis.text = list(font = 2, cex = 1.3), |
|
| 455 |
# par.main.text=list(font=2,cex=2),strip.background=list(col="white")),par.strip.text=list(font=2,cex=1.5)) |
|
| 456 |
#p_dif <- update(p_dif, panel = function(...) {
|
|
| 457 |
# panel.abline(h = 0, lty = 2, col = "gray") |
|
| 458 |
# panel.xyplot(...) |
|
| 459 |
# }) |
|
| 421 | 460 |
|
| 461 |
dev.off() |
|
| 422 | 462 |
|
| 423 | 463 |
####################################################### |
| 424 | 464 |
####### PART 2: generate figures for paper ############# |
| ... | ... | |
| 680 | 720 |
list_plots_spt <- vector("list",length=length(lf))
|
| 681 | 721 |
for (i in 1:length(lf)){
|
| 682 | 722 |
pred_temp_s <-stack(lf[[i]]) |
| 683 |
|
|
| 723 |
pred_temp_s <- projectRaster(pred_temp_s,crs=CRS_WGS84) |
|
| 684 | 724 |
#lf2 <- raster_prediction_obj_2$method_mod_obj[[index]][[y_var_name]] |
| 685 | 725 |
#lf2 #contains the models for gam |
| 686 | 726 |
|
| ... | ... | |
| 699 | 739 |
#temp.colors <- colorRampPalette(c('blue', 'lightgoldenrodyellow', 'red'))
|
| 700 | 740 |
#temp.colors <- matlab.like(no_brks) |
| 701 | 741 |
temp.colors <- colorRampPalette(c('blue', 'khaki', 'red'))
|
| 742 |
#temp.colors <- colorRampPalette(c('khaki', 'darkred'))
|
|
| 702 | 743 |
|
| 703 | 744 |
png(paste("Figure_",nb_fig[i],"_spatial_pattern_tmax_prediction_models_gam_levelplot_",date_selected,out_prefix,".png", sep=""),
|
| 704 | 745 |
height=480*layout_m[1],width=480*layout_m[2]) |
| 705 | 746 |
|
| 706 |
p <- levelplot(pred_temp_s,main=methods_names[i],
|
|
| 747 |
p <- levelplot(pred_temp_s,main=paste(methods_names[i],"temperature (°C)",sep=" "),
|
|
| 707 | 748 |
ylab=NULL,xlab=NULL, |
| 708 | 749 |
par.settings = list(axis.text = list(font = 2, cex = 2),layout=layout_m, |
| 709 |
par.main.text=list(font=2,cex=2.5),strip.background=list(col="white")),par.strip.text=list(font=2,cex=2), |
|
| 750 |
par.main.text=list(font=2,cex=2.5),strip.background=list(col="white")), |
|
| 751 |
par.strip.text=list(font=2,cex=2), |
|
| 710 | 752 |
names.attr=names_layers, |
| 711 |
col.regions=temp.colors(no_brks),at=seq(min_val,max_val,by=0.1)) |
|
| 753 |
col.regions=temp.colors(no_brks), |
|
| 754 |
at=seq(min_val,max_val,by=0.1), |
|
| 755 |
colorkey = list(space = "right",legend), |
|
| 756 |
#main = "Temp (C)" |
|
| 757 |
) |
|
| 758 |
|
|
| 712 | 759 |
#col.regions=temp.colors(25)) |
| 713 | 760 |
print(p) |
| 714 | 761 |
dev.off() |
| ... | ... | |
| 913 | 960 |
p2 <- list_plot_std_moranI[[2]][[2]] #std for gam CAI |
| 914 | 961 |
|
| 915 | 962 |
p1$main <- "Moran I in daily predictions averaged by month" |
| 916 |
p2$main <- "Standard devation in daily predictions averaged by month" |
|
| 917 |
p3$main <- "Moran I devation in daily predictions average by month"
|
|
| 918 |
p4$main <- "Standard devation in daily predictions average by month"
|
|
| 963 |
p2$main <- "Standard deviation in daily predictions averaged by month"
|
|
| 964 |
p3$main <- "Moran I in daily predictions averaged by month"
|
|
| 965 |
p4$main <- "Standard deviation in daily predictions averaged by month"
|
|
| 919 | 966 |
|
| 920 | 967 |
p1$legend$right$args$title <- "gam_CAI" |
| 921 | 968 |
p2$legend$right$args$title <- "gam_CAI" |
| ... | ... | |
| 941 | 988 |
index <- 289 |
| 942 | 989 |
temp_day_s <- stack(lf_clim_gam_CAI[[index]]$mod7,lf_delta_gam_CAI[[index]]$mod7,lf_daily_gam_CAI[[index]]$mod7) |
| 943 | 990 |
names(temp_day_s) <- c("Monhtly Climatology","Daily Devation","Daily Prediction")
|
| 944 |
|
|
| 991 |
temp_day_s <- projectRaster(temp_day_s,crs=CRS_WGS84) |
|
| 945 | 992 |
|
| 946 | 993 |
no_brks=25 |
| 947 | 994 |
temp.colors <- colorRampPalette(c('blue', 'khaki', 'red'))
|
| ... | ... | |
| 955 | 1002 |
ylab=NULL,xlab=NULL, |
| 956 | 1003 |
par.settings = list(axis.text = list(font = 2, cex = 1.5), |
| 957 | 1004 |
par.main.text=list(font=2,cex=2.5),strip.background=list(col="white")),par.strip.text=list(font=2,cex=2), |
| 958 |
main=names_layers[i],col.regions=temp.colors(no_brks)) |
|
| 1005 |
main=paste(names_layers[i],"temp. (°C)",sep=" "), |
|
| 1006 |
col.regions=temp.colors(no_brks)) |
|
| 959 | 1007 |
|
| 960 | 1008 |
png(paste("Figure",fig_nb[i],"_paper_climatology_daily_deviation_daily_prediction_model7_levelplot_",out_prefix,".png", sep=""),
|
| 961 | 1009 |
height=480*1,width=480*1) |
| ... | ... | |
| 1019 | 1067 |
data_month_all <-convert_spdf_to_df_from_list(list_data_month_s) #long rownames |
| 1020 | 1068 |
#LSTD_bias_avgm<-aggregate(LSTD_bias~month,data=data_month_all,mean) |
| 1021 | 1069 |
#LSTD_bias_sdm<-aggregate(LSTD_bias~month,data=data_month_all,sd) |
| 1022 |
LST_avgm <- aggregate(LST~month,data=data_month_all,mean) |
|
| 1023 | 1070 |
LST_avgm_min <- aggregate(LST~month,data=data_month_all,min) |
| 1024 | 1071 |
LST_avgm_max <- aggregate(LST~month,data=data_month_all,max) |
| 1072 |
TMax_avgm_min <- aggregate(TMax~month,data=data_month_all,min) |
|
| 1073 |
TMax_avgm_max <- aggregate(TMax~month,data=data_month_all,max) |
|
| 1074 |
|
|
| 1075 |
LST_avgm <- aggregate(LST~month,data=data_month_all,mean) |
|
| 1025 | 1076 |
TMax_avgm <- aggregate(TMax~month,data=data_month_all,mean) |
| 1026 | 1077 |
|
| 1027 | 1078 |
#plot(LST_avgm,type="b") |
| ... | ... | |
| 1031 | 1082 |
statistics_LST_s<- LST_stat$avg #extracted earlier!!! |
| 1032 | 1083 |
|
| 1033 | 1084 |
##Now plot Figure 8a |
| 1034 |
plot(TMax_avgm,type="b",ylim=c(-10,50),col="red",xlab="month", |
|
| 1035 |
ylab="tmax (degree C)",cex=2,cex.lab=1.5) #,cex.axis=1.8) |
|
| 1036 |
lines(1:12,LST_avgm$LST,type="b",col="green",lty="dashed",pch=2,cex=2) |
|
| 1085 |
plot(TMax_avgm,type="b",ylim=c(-10,50),col="red",lty="dotted",xlab="month", |
|
| 1086 |
ylab="tmax (degree C)",pch=2,cex=2,cex.lab=1.5) #,cex.axis=1.8) |
|
| 1087 |
|
|
| 1088 |
lines(1:12,TMax_avgm_min$TMax,type="b",col="red",lty="dotted",pch=3,cex=2) |
|
| 1089 |
lines(1:12,TMax_avgm_max$TMax,type="b",col="red",lty="dotted",pch=5,cex=2) |
|
| 1090 |
|
|
| 1091 |
lines(1:12,LST_avgm$LST,type="b",col="blue",lty="dashed",pch=2,cex=2) |
|
| 1037 | 1092 |
#lines(1:12,statistics_LST_s$mean,type="b",col="darkgreen") |
| 1038 |
lines(1:12,LST_avgm_min$LST,type="b",col="black",lty="dashed",pch=3,cex=2) |
|
| 1039 |
lines(1:12,LST_avgm_max$LST,type="b",col="black",lty="dashed",pch=4,cex=2) |
|
| 1093 |
lines(1:12,LST_avgm_min$LST,type="b",col="blue",lty="dashed",pch=3,cex=2) |
|
| 1094 |
lines(1:12,LST_avgm_max$LST,type="b",col="blue",lty="dashed",pch=5,cex=2) |
|
| 1095 |
|
|
| 1040 | 1096 |
#text(1:12,LST_avgm$LST,month.abb,cex=0.8,pos=2) |
| 1041 |
legend("topleft",legend=c("TMax","LST","LST min","LST max"), cex=1.5,
|
|
| 1042 |
col=c("red","green","black","black"),
|
|
| 1043 |
lty=c("solid","dashed","dashed","dashed"),pch=1:4,bty="n")
|
|
| 1044 |
#legend("topright",legend=c("TMax","LST","LST min","LST max"), cex=1.5,
|
|
| 1097 |
#legend("topleft",legend=c("TMax_avg","LST_avg","LST_min","LST_max"), cex=1.5,
|
|
| 1045 | 1098 |
# col=c("red","green","black","black"),
|
| 1046 | 1099 |
# lty=c("solid","dashed","dashed","dashed"),pch=1:4,bty="n")
|
| 1100 |
legend("topright",legend=c("LST_mean","LST_min","LST_max"), cex=1.5,
|
|
| 1101 |
col=c("blue","blue","blue"),
|
|
| 1102 |
lty=c("dashed","dashed","dashed"),
|
|
| 1103 |
pch=c(2,3,5),bty="n") |
|
| 1104 |
legend("topleft",legend=c("TMax_mean","TMax_min","TMax_max"), cex=1.5,
|
|
| 1105 |
col=c("red","red","red"),
|
|
| 1106 |
lty=c("dotted","dotted","dotted"),
|
|
| 1107 |
pch=c(2,3,5),bty="n") |
|
| 1108 |
|
|
| 1047 | 1109 |
title("Monthly average LST and tmax at stations in Oregon",cex=2.4)
|
| 1048 | 1110 |
|
| 1049 | 1111 |
############################## |
| ... | ... | |
| 1060 | 1122 |
LST_avgm_grass <-aggregate(LST~month,data=data_month_all_grass,mean) |
| 1061 | 1123 |
LST_avgm_urban <-aggregate(LST~month,data=data_month_all_urban,mean) |
| 1062 | 1124 |
|
| 1125 |
TMax_avgm_forest <-aggregate(TMax~month,data=data_month_all_forest,mean) |
|
| 1126 |
TMax_avgm_grass <-aggregate(TMax~month,data=data_month_all_grass,mean) |
|
| 1127 |
TMAx_avgm_urban <-aggregate(TMax~month,data=data_month_all_urban,mean) |
|
| 1128 |
|
|
| 1063 | 1129 |
##Now plot Figure 8b |
| 1064 | 1130 |
plot(TMax_avgm,type="b",ylim=c(-7,50),col="red",xlab="month", |
| 1065 |
ylab="tmax (degree C)",cex=2,cex.lab=1.5)#,cex.axis=1.8) |
|
| 1066 |
lines(1:12,LST_avgm$LST,type="b",col="green",lty="dashed",pch=2,cex=2) |
|
| 1067 |
lines(LST_avgm_forest,type="b",col="black",lty="dotted",pch=3,cex=2) |
|
| 1068 |
lines(LST_avgm_grass,type="b",col="black",lty="dotdash",pch=4,cex=2) |
|
| 1131 |
ylab="tmax (degree C)",pch=2,cex=2,cex.lab=1.5,lty="dashed")#,cex.axis=1.8) |
|
| 1132 |
lines(TMax_avgm_forest,type="b",col="red",lty="dashed",pch=3,cex=2) |
|
| 1133 |
lines(TMax_avgm_grass,type="b",col="red",lty="dashed",pch=5,cex=2) |
|
| 1134 |
|
|
| 1135 |
lines(1:12,LST_avgm$LST,type="b",col="blue",lty="dotted",pch=2,cex=2) |
|
| 1136 |
lines(LST_avgm_forest,type="b",col="blue",lty="dotted",pch=3,cex=2) |
|
| 1137 |
lines(LST_avgm_grass,type="b",col="blue",lty="dotted",pch=5,cex=2) |
|
| 1069 | 1138 |
#lines(LST_avgm_urban,type="b",col="black",lty="longdash",pch=4) |
| 1070 |
legend("topleft",legend=c("TMax","LST","LST_forest","LST_grass"), cex=1.5,
|
|
| 1071 |
col=c("red","green","black","black","black"),
|
|
| 1072 |
lty=1:4,pch=1:4,bty="n") |
|
| 1139 |
#legend("topleft",legend=c("TMax","LST","LST_forest","LST_grass"), cex=1.5,
|
|
| 1140 |
# col=c("red","green","black","black","black"),
|
|
| 1141 |
# lty=1:4,pch=1:4,bty="n") |
|
| 1142 |
|
|
| 1143 |
legend("topleft",legend=c("TMax_mean","TMax_forest","TMax_grass"), cex=1.5,
|
|
| 1144 |
col=c("red","red","red"),
|
|
| 1145 |
lty=c("dotted","dotted","dotted"),
|
|
| 1146 |
pch=c(2,3,5),bty="n") |
|
| 1147 |
|
|
| 1148 |
legend("topright",legend=c("LST_mean","LST_forest","LST_grass"), cex=1.5,
|
|
| 1149 |
col=c("blue","blue","blue"),
|
|
| 1150 |
lty=c("dashed","dashed","dashed"),
|
|
| 1151 |
pch=c(2,3,5),bty="n") |
|
| 1152 |
|
|
| 1073 | 1153 |
title("Monthly average LST and land cover types for stations ",cex=2.4)
|
| 1074 | 1154 |
|
| 1075 | 1155 |
######################################## |
| ... | ... | |
| 1089 | 1169 |
names(r_tmp2) <- 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","LC1")
|
| 1090 | 1170 |
x_cor_tmp2 <-layerStats(r_tmp2,stat="pearson",na.rm=T) |
| 1091 | 1171 |
|
| 1092 |
##Now plot Figure 8d
|
|
| 1172 |
##Now plot Figure 8c
|
|
| 1093 | 1173 |
plot(1:12,x_cor_tmp1[[1]][1:12,13],ylim=c(-0.9,0.9), |
| 1094 | 1174 |
type="b",lty="solid",pch=1, |
| 1095 | 1175 |
ylab="Pearson corelation",xlab="month", |
| ... | ... | |
| 1133 | 1213 |
|
| 1134 | 1214 |
lst_md <- raster(ref_rast_name) |
| 1135 | 1215 |
projection(lst_md)<-projection(s_raster) |
| 1136 |
lst_mm_09<-subset(s_raster,"mm_09") |
|
| 1216 |
lst_mm_09 <-subset(s_raster,"mm_09") |
|
| 1217 |
lst_mm_09 <- projectRaster(lst_mm_09,crs=CRS_WGS84) |
|
| 1137 | 1218 |
|
| 1138 | 1219 |
lst_md<-raster(ref_rast_d001) |
| 1220 |
projection(lst_md) <- proj4string(elev) |
|
| 1221 |
lst_md <- projectRaster(lst_md,crs=CRS_WGS84) |
|
| 1222 |
|
|
| 1139 | 1223 |
lst_md<- lst_md - 273.16 |
| 1140 | 1224 |
lst_mm_01<-subset(s_raster,"mm_01") |
| 1225 |
lst_mm_01 <- projectRaster(lst_mm_01,crs=CRS_WGS84) |
|
| 1141 | 1226 |
|
| 1142 | 1227 |
png(filename=paste("Figure_annex1_paper_Comparison_daily_monthly_mean_lst",out_prefix,".png",sep=""),width=960,height=480)
|
| 1143 | 1228 |
par(mfrow=c(1,2)) |
| 1144 | 1229 |
plot(lst_md) |
| 1145 |
plot(interp_area,add=TRUE) |
|
| 1146 |
title("Mean for January 1")
|
|
| 1230 |
plot(interp_area_WGS84,add=TRUE)
|
|
| 1231 |
title("Mean temperature (°C) for January 1")
|
|
| 1147 | 1232 |
plot(lst_mm_01) |
| 1148 |
plot(interp_area,add=TRUE) |
|
| 1149 |
title("Mean for month of January")
|
|
| 1233 |
plot(interp_area_WGS84,add=TRUE)
|
|
| 1234 |
title("Mean temperature (°C) for month of January")
|
|
| 1150 | 1235 |
dev.off() |
| 1151 | 1236 |
|
| 1152 | 1237 |
############################################ |
| ... | ... | |
| 1156 | 1241 |
layout_m<-c(1,1) |
| 1157 | 1242 |
png(paste("Figure_annex2_paper_elevation_transect_paths_",out_prefix,".png", sep=""),
|
| 1158 | 1243 |
height=480*layout_m[1],width=480*layout_m[2]) |
| 1159 |
vx_text <- c(395770.1,395770.1,395770.1) #location of labels |
|
| 1160 |
vy_text <-c(473000,297045.9,112611.5) |
|
| 1161 |
|
|
| 1162 |
plot(elev) |
|
| 1244 |
#vx_text <- c(395770.1,395770.1,395770.1) #location of labels |
|
| 1245 |
#vy_text <-c(473000,297045.9,112611.5) |
|
| 1246 |
vx_text <- c(-122,-122,-122) |
|
| 1247 |
vy_text <- c(45.9,44.4,42.7) |
|
| 1248 |
|
|
| 1249 |
#plot(elev) |
|
| 1250 |
#plot(elev_WGS84,cex.axis=2.1,cex.z=2.1) |
|
| 1251 |
plot(elev_WGS84,cex.axis=1.4,ylab="latitude (degree)",xlab="longitude (degree)", |
|
| 1252 |
cex.lab=1.4, |
|
| 1253 |
axis.args=list(at=seq(0,3500,500), |
|
| 1254 |
labels=seq(0, 3500, 500), |
|
| 1255 |
cex.axis=1.4)) |
|
| 1256 |
title_text <-c("Elevation (m) and spatial transects"," in Oregon")
|
|
| 1257 |
title(title_text) |
|
| 1163 | 1258 |
for(i in 1:length(transect_list)){
|
| 1164 | 1259 |
filename<-sub(".shp","",transect_list[i]) #Removing the extension from file.
|
| 1165 | 1260 |
transect<-readOGR(dirname(filename), basename(filename)) #reading shapefile |
| 1261 |
transect<- spTransform(transect,CRS=CRS_WGS84) |
|
| 1166 | 1262 |
#transect2 <- elide(transect, shift=c(0, -24000)) |
| 1167 | 1263 |
#plot(transect2,add=TRUE) |
| 1168 | 1264 |
#writeOGR(transect2,dsn=".",layer="transect_OR_1_12282013",driver="ESRI Shapefile") |
| 1169 | 1265 |
plot(transect,add=TRUE) |
| 1170 | 1266 |
} |
| 1171 | 1267 |
text(x=vx_text,y=vy_text,labels=c("Transect 1","Transect 2","Transect 3"))
|
| 1172 |
title("Transect Oregon")
|
|
| 1268 |
#title("Transects in Oregon with elevation (m) in the background")
|
|
| 1173 | 1269 |
dev.off() |
| 1174 | 1270 |
|
| 1175 | 1271 |
############################## |
| ... | ... | |
| 1248 | 1344 |
trans_data2 <- plot_transect_m2(list_transect2,rast_pred2,title_plot2,disp=FALSE,m_layers_sc) |
| 1249 | 1345 |
trans_data3 <- plot_transect_m2(list_transect3,rast_pred3,title_plot2,disp=FALSE,m_layers_sc) |
| 1250 | 1346 |
|
| 1251 |
|
|
| 1252 |
### REVISIONS 1 PAPER ##### |
|
| 1347 |
############################## |
|
| 1348 |
####### REVISIONS 1 PAPER #####
|
|
| 1253 | 1349 |
|
| 1254 | 1350 |
## DATA INFORMATION ON STATIONS |
| 1255 | 1351 |
|
| ... | ... | |
| 1280 | 1376 |
|
| 1281 | 1377 |
p <- histogram(dsts_nb1) |
| 1282 | 1378 |
|
| 1283 |
hist(ghcnd_day_dat$elev_s) |
|
| 1379 |
hist(ghcn_day_dat$elev_s) |
|
| 1380 |
histogram(ghcn_day_dat$elev_s) |
|
| 1284 | 1381 |
|
| 1285 | 1382 |
mean(ghcn_day_dat) |
| 1286 | 1383 |
unique(ghcn_day_dat$date) # |
| ... | ... | |
| 1288 | 1385 |
hist(table(ghcn_day_dat$date)) |
| 1289 | 1386 |
range(tb_n_day$Freq) #range 134 to 159 |
| 1290 | 1387 |
median(tb_n_day$Freq) #range 134 to 159 |
| 1388 |
mean(tb_n_day$Freq) #range 134 to 159 |
|
| 1291 | 1389 |
|
| 1292 | 1390 |
|
| 1293 | 1391 |
## EXAMINE MONTHLY DAta |
| ... | ... | |
| 1342 | 1440 |
# "gam_CAI","kriging_CAI","gwr_CAI", |
| 1343 | 1441 |
# "gam_fss","kriging_fss","gwr_fss") |
| 1344 | 1442 |
|
| 1345 |
list_data_v <-lapply(list_raster_obj_files[5:7],FUN=function(x){x<-load_obj(x);extract_list_from_list_obj(x$validation_mod_obj,"data_v")})
|
|
| 1443 |
#get data for gam_CAI,kriging_CAI, gwr_CAI |
|
| 1444 |
list_data_v_all <-lapply(list_raster_obj_files[5:7],FUN=function(x){x<-load_obj(x);extract_list_from_list_obj(x$validation_mod_obj,"data_v")})
|
|
| 1346 | 1445 |
#list_data_s <-lapply(list_raster_obj_files[5:7],FUN=function(x){x<-load_obj(x);extract_list_from_list_obj(x$validation_mod_obj,"data_s")})
|
| 1347 | 1446 |
|
| 1348 | 1447 |
met_obj <-load_obj(file.path(in_dir1,met_obj_file_1)) |
| ... | ... | |
| 1357 | 1456 |
interp_method <- "gam_CAI" |
| 1358 | 1457 |
list_data_v <- list_data_v_all$gam_CAI |
| 1359 | 1458 |
data_mae_gam_CAI_obj <- plot_MAE_per_station_fun(list_data_v,names_var,interp_method,var_background,stat_loc,out_suffix) |
| 1459 |
#debug(plot_MAE_per_station_fun) |
|
| 1460 |
#data_mae_gam_CAI_obj <- plot_MAE_per_station_fun(list_data_v,names_var,interp_method,var_background,stat_loc,out_suffix) |
|
| 1360 | 1461 |
|
| 1361 | 1462 |
interp_method <- "kriging_CAI" |
| 1362 | 1463 |
list_data_v <- list_data_v_all$kriging_CAI |
| ... | ... | |
| 1371 | 1472 |
|
| 1372 | 1473 |
list_p_mae <- data_mae_gam_CAI_obj$list_p_mae |
| 1373 | 1474 |
interp_method <- "gam_CAI" |
| 1374 |
layout_m<-c(1,4) # works if set to this?? ok set the resolution... |
|
| 1475 |
#layout_m<-c(1,4) # works if set to this?? ok set the resolution... |
|
| 1476 |
#layout_m<-c(4,1) # works if set to this?? ok set the resolution... |
|
| 1477 |
|
|
| 1478 |
#This is for model 1, model2, mod3 and mod |
|
| 1479 |
#png(paste("Figure10_paper_","average_MAE_",interp_method,out_prefix,".png", sep=""),
|
|
| 1480 |
# height=480*layout_m[1],width=480*layout_m[2]) |
|
| 1481 |
#grid.arrange(list_p_mae[[1]],list_p_mae[[2]],list_p_mae[[3]],list_p_mae[[7]], ncol=1) |
|
| 1482 |
#dev.off() |
|
| 1483 |
|
|
| 1484 |
#This is for model 1, model2, mod3 and mod7 |
|
| 1485 |
|
|
| 1486 |
layout_m<-c(1,3) # works if set to this?? ok set the resolution... |
|
| 1375 | 1487 |
|
| 1376 |
png(paste("Figure10_paper_","average_MAE_",interp_method,out_prefix,".png", sep=""),
|
|
| 1488 |
png(paste("Figure10_paper_","average_MAE_","mod1_",out_prefix,".png", sep=""),
|
|
| 1377 | 1489 |
height=480*layout_m[1],width=480*layout_m[2]) |
| 1378 |
grid.arrange(list_p_mae[[1]],list_p_mae[[2]],list_p_mae[[3]],list_p_mae[[7]], ncol=4) |
|
| 1490 |
grid.arrange(data_mae_gam_CAI_obj$list_p_mae[[1]],data_mae_kriging_CAI_obj$list_p_mae[[1]],data_mae_gwr_CAI_obj$list_p_mae[[1]] |
|
| 1491 |
, ncol=3) |
|
| 1379 | 1492 |
dev.off() |
| 1380 | 1493 |
|
| 1381 |
list_p_mae <- data_mae_kriging_CAI_obj$list_p_mae |
|
| 1382 |
layout_m<-c(1,4) # works if set to this?? ok set the resolution... |
|
| 1383 |
interp_method <- "kriging_CAI" |
|
| 1384 |
|
|
| 1385 |
png(paste("Figure10_paper_","average_MAE_",interp_method,out_prefix,".png", sep=""),
|
|
| 1494 |
#This is for model 1, model2, mod3 and mod |
|
| 1495 |
png(paste("Figure10_paper_","average_MAE_","mod2_",out_prefix,".png", sep=""),
|
|
| 1386 | 1496 |
height=480*layout_m[1],width=480*layout_m[2]) |
| 1387 |
grid.arrange(list_p_mae[[1]],list_p_mae[[2]],list_p_mae[[3]],list_p_mae[[7]], ncol=4) |
|
| 1497 |
grid.arrange(data_mae_gam_CAI_obj$list_p_mae[[2]],data_mae_kriging_CAI_obj$list_p_mae[[2]],data_mae_gwr_CAI_obj$list_p_mae[[2]] |
|
| 1498 |
, ncol=3) |
|
| 1388 | 1499 |
dev.off() |
| 1389 | 1500 |
|
| 1390 |
list_p_mae <- data_mae_gwr_CAI_obj$list_p_mae |
|
| 1391 |
layout_m<-c(1,4) # works if set to this?? ok set the resolution... |
|
| 1392 |
interp_method <- "gwr_CAI" |
|
| 1501 |
#This is for model 1, model2, mod3 and mod |
|
| 1502 |
png(paste("Figure10_paper_","average_MAE_","mod3_",out_prefix,".png", sep=""),
|
|
| 1503 |
height=480*layout_m[1],width=480*layout_m[2]) |
|
| 1504 |
grid.arrange(data_mae_gam_CAI_obj$list_p_mae[[3]],data_mae_kriging_CAI_obj$list_p_mae[[3]],data_mae_gwr_CAI_obj$list_p_mae[[3]] |
|
| 1505 |
, ncol=3) |
|
| 1506 |
dev.off() |
|
| 1393 | 1507 |
|
| 1394 |
png(paste("Figure10_paper_","average_MAE_",interp_method,out_prefix,".png", sep=""),
|
|
| 1508 |
#This is for model 1, model2, mod3 and mod |
|
| 1509 |
png(paste("Figure10_paper_","average_MAE_","mod7_",out_prefix,".png", sep=""),
|
|
| 1395 | 1510 |
height=480*layout_m[1],width=480*layout_m[2]) |
| 1396 |
grid.arrange(list_p_mae[[1]],list_p_mae[[2]],list_p_mae[[3]],list_p_mae[[7]], ncol=4) |
|
| 1511 |
grid.arrange(data_mae_gam_CAI_obj$list_p_mae[[7]],data_mae_kriging_CAI_obj$list_p_mae[[7]],data_mae_gwr_CAI_obj$list_p_mae[[7]] |
|
| 1512 |
, ncol=3) |
|
| 1397 | 1513 |
dev.off() |
| 1398 | 1514 |
|
| 1515 |
#list_p_mae <- data_mae_kriging_CAI_obj$list_p_mae |
|
| 1516 |
#layout_m<-c(1,4) # works if set to this?? ok set the resolution... |
|
| 1517 |
#interp_method <- "kriging_CAI" |
|
| 1518 |
|
|
| 1519 |
#png(paste("Figure10_paper_","average_MAE_",interp_method,out_prefix,".png", sep=""),
|
|
| 1520 |
# height=480*layout_m[1],width=480*layout_m[2]) |
|
| 1521 |
#grid.arrange(list_p_mae[[1]],list_p_mae[[2]],list_p_mae[[3]],list_p_mae[[7]], ncol=4) |
|
| 1522 |
#dev.off() |
|
| 1523 |
|
|
| 1524 |
#list_p_mae <- data_mae_gwr_CAI_obj$list_p_mae |
|
| 1525 |
#layout_m<-c(1,4) # works if set to this?? ok set the resolution... |
|
| 1526 |
#interp_method <- "gwr_CAI" |
|
| 1527 |
|
|
| 1528 |
#png(paste("Figure10_paper_","average_MAE_",interp_method,out_prefix,".png", sep=""),
|
|
| 1529 |
# height=480*layout_m[1],width=480*layout_m[2]) |
|
| 1530 |
#grid.arrange(list_p_mae[[1]],list_p_mae[[2]],list_p_mae[[3]],list_p_mae[[7]], ncol=4) |
|
| 1531 |
#dev.off() |
|
| 1532 |
|
|
| 1399 | 1533 |
##### |
| 1400 | 1534 |
## MAE Histogram plots: Show histo for mod1,mod2,mod3 and mod7 for gam, kriging and gwr |
| 1401 | 1535 |
|
| ... | ... | |
| 1502 | 1636 |
### |
| 1503 | 1637 |
|
| 1504 | 1638 |
|
| 1639 |
#http://www.carlislerainey.com/2012/12/28/how-to-add-an-extra-vertical-axis-to-r-plots/ |
|
| 1640 |
# I make the following changes. |
|
| 1641 |
# |
|
| 1642 |
# 1) Adjust the mar option in the graphical parameters to give me a little more room for constructing the vertical axis on the right side. |
|
| 1643 |
# 2) Add the xlab="x" argument to the first plot() function call. This improves the label of the vertical axis on the left. |
|
| 1644 |
# 3) Add the axes = F, xlab = NA, and ylab = NA arguments to the second plot() function call. This removes all the axes and axis notation from this call. |
|
| 1645 |
# 4) Add the vertical axis on the right with a simple call to the axis() function and setting side = 4. |
|
| 1646 |
# 5) Add a label to the vertical axis on the right with the mtext() command. Use the options side = 4 and line = 3 to position the text "y" correctly. |
|
| 1647 |
|
|
| 1648 |
#Adding right side label |
|
| 1649 |
#par(mar = c(5,5,2,5)) |
|
| 1650 |
#plot(ts(x), ylab = "x") |
|
| 1651 |
#par(new = T) |
|
| 1652 |
#plot(ts(y), col = "red", axes = F, xlab = NA, ylab = NA) |
|
| 1653 |
#axis(side = 4) |
|
| 1654 |
#mtext(side = 4, line = 3, "y") |
|
Also available in: Unified diff
revisions2 multitime scale paper, figures and difference analysis