Revision 4306add2
Added by Benoit Parmentier about 12 years ago
| climate/research/oregon/interpolation/methods_comparison_assessment_part7.R | ||
|---|---|---|
| 1 | 1 |
##################################### METHODS COMPARISON part 7 ########################################## |
| 2 | 2 |
#################################### Spatial Analysis: validation CAI-fusion ############################################ |
| 3 | 3 |
#This script utilizes the R ojbects created during the interpolation phase. # |
| 4 |
#We use the SNOTEL dataset and the |
|
| 5 |
#This scripts focuses on a detailed studay of differences in the predictions of CAI_kr and FUsion_Kr #
|
|
| 4 |
#We use the SNOTEL dataset and the GHCN network to assess the prediction accuracy.
|
|
| 5 |
#This scripts focuses on a detailed study of differences in the predictions of CAI_kr and FUsion_Kr # |
|
| 6 | 6 |
#AUTHOR: Benoit Parmentier # |
| 7 | 7 |
#DATE: 12/03/2012 # |
| 8 | 8 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#491 -- # |
| ... | ... | |
| 27 | 27 |
library(reshape) |
| 28 | 28 |
library(RCurl) |
| 29 | 29 |
######### Functions used in the script |
| 30 |
|
|
| 31 |
load_obj <- function(f) |
|
| 32 |
{
|
|
| 33 |
env <- new.env() |
|
| 34 |
nm <- load(f, env)[1] |
|
| 35 |
env[[nm]] |
|
| 36 |
} |
|
| 37 |
|
|
| 38 |
format_padding_month<-function(date_str){
|
|
| 39 |
date_trans<-character(length=length(date_str)) |
|
| 40 |
for (i in 1:length(date_str)){
|
|
| 41 |
tmp_date<-date_str[i] |
|
| 42 |
nc<-nchar(tmp_date) |
|
| 43 |
nstart<-nc-1 |
|
| 44 |
year<-substr(tmp_date,start=nstart,stop=nc) |
|
| 45 |
md<-substr(tmp_date,start=1,stop=(nstart-1)) |
|
| 46 |
if (nchar(md)==3){
|
|
| 47 |
md<-paste("0",md,sep="")
|
|
| 48 |
} |
|
| 49 |
date_trans[i]<-paste(md,year,sep="") |
|
| 50 |
} |
|
| 51 |
return(date_trans) |
|
| 52 |
} |
|
| 53 |
|
|
| 54 |
merge_multiple_df<-function(df_list,by_name){
|
|
| 55 |
for (i in 1:(length(df_list)-1)){
|
|
| 56 |
if (i==1){
|
|
| 57 |
df1=df_list[[i]] |
|
| 58 |
} |
|
| 59 |
if (i!=1){
|
|
| 60 |
df1=df_m |
|
| 61 |
} |
|
| 62 |
df2<-df_list[[i+1]] |
|
| 63 |
df_m<-merge(df1,df2,by=by_name,all=T) |
|
| 64 |
} |
|
| 65 |
return(df_m) |
|
| 66 |
} |
|
| 67 |
|
|
| 68 |
reclassify_df<-function(df,var_name,brks,lab_brks,suffix,summary_var){
|
|
| 69 |
var_tab<-vector("list",length(var_name))
|
|
| 70 |
for (i in 1:length(var_name)){
|
|
| 71 |
var_rec_name<-paste(var_name[i],suffix,sep="_") |
|
| 72 |
var_rcstat<-cut(df[[var_name[i]]],breaks=brks,labels=lab_brks,right=T) |
|
| 73 |
df[[var_rec_name]]<-var_rcstat |
|
| 74 |
tmp<-aggregate(df[[summary_var]]~df[[var_rec_name]],data=df,FUN=mean) |
|
| 75 |
names(tmp)<-c(suffix,var_rec_name) |
|
| 76 |
var_tab[[i]]<-tmp |
|
| 77 |
} |
|
| 78 |
obj<-list(var_tab,df) |
|
| 79 |
names(obj)<-c("agg_df","df")
|
|
| 80 |
return(list(var_tab,df)) |
|
| 81 |
} |
|
| 82 |
|
|
| 83 |
station_data_interp<-function(date_str,obj_mod_interp_str,training=TRUE,testing=TRUE){
|
|
| 84 |
date_selected<-date_str |
|
| 85 |
#load interpolation object |
|
| 86 |
obj_mod_interp<-load_obj(obj_mod_interp_str) |
|
| 87 |
sampling_date_list<-obj_mod_interp$sampling_obj$sampling_dat$date |
|
| 88 |
k<-match(date_selected,sampling_date_list) |
|
| 89 |
names(obj_mod_interp[[1]][[k]]) #Show the name structure of the object/list |
|
| 90 |
|
|
| 91 |
#Extract the training and testing information for the given date... |
|
| 92 |
data_s<-obj_mod_interp[[1]][[k]]$data_s #object for the first date...20100103 |
|
| 93 |
data_v<-obj_mod_interp[[1]][[k]]$data_v #object for the first date...20100103 |
|
| 94 |
if (testing==TRUE & training==FALSE){
|
|
| 95 |
return(data_v) |
|
| 96 |
} |
|
| 97 |
if (training==TRUE & testing==FALSE){
|
|
| 98 |
return(data_s) |
|
| 99 |
} |
|
| 100 |
if (training==TRUE & testing==TRUE ){
|
|
| 101 |
dataset_stat<-list(data_v,data_s) |
|
| 102 |
names(dataset_stat)<-c("testing","training")
|
|
| 103 |
return(dataset_stat) |
|
| 104 |
} |
|
| 105 |
} |
|
| 106 |
|
|
| 107 |
### Caculate accuracy metrics |
|
| 108 |
calc_accuracy_metrics<-function(x,y){
|
|
| 109 |
residuals<-x-y |
|
| 110 |
mae<-mean(abs(residuals),na.rm=T) |
|
| 111 |
rmse<-sqrt(mean((residuals)^2,na.rm=T)) |
|
| 112 |
me<-mean(residuals,na.rm=T) |
|
| 113 |
r<-cor(x,y,use="complete") |
|
| 114 |
avg<-mean(residuals,na.rm=T) |
|
| 115 |
m50<-median(residuals,na.rm=T) |
|
| 116 |
metrics_dat<-as.data.frame(cbind(mae,rmse,me,r,avg,m50)) |
|
| 117 |
names(metrics_dat)<-c("mae","rmse","me","r","avg","m50")
|
|
| 118 |
return(metrics_dat) |
|
| 119 |
} |
|
| 120 |
|
|
| 121 |
#MODIFY LATER |
|
| 122 |
# raster_pred_interp<-function(date_str,rast_file_name_list,path_data,data_sp){
|
|
| 123 |
# date_selected<-date_str |
|
| 124 |
# #load interpolation object |
|
| 125 |
# setwd(path_data) |
|
| 126 |
# file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_CAI2_const_all_10312012.rst",sep="")) #Search for files in relation to fusion
|
|
| 127 |
# lf_pred<-list.files(pattern=file_pat) #Search for files in relation to fusion |
|
| 128 |
# |
|
| 129 |
# rast_cai2c<-stack(lf_cai2c) #lf_cai2c CAI results with constant sampling over 365 dates |
|
| 130 |
# rast_cai2c<-mask(rast_cai2c,mask_ELEV_SRTM) |
|
| 131 |
# |
|
| 132 |
# obj_mod_interp<-load_obj(obj_mod_interp_str) |
|
| 133 |
# sampling_date_list<-obj_mod_interp$sampling_obj$sampling_dat$date |
|
| 134 |
# k<-match(date_selected,sampling_date_list) |
|
| 135 |
# names(obj_mod_interp[[1]][[k]]) #Show the name structure of the object/list |
|
| 136 |
# |
|
| 137 |
# #Extract the training and testing information for the given date... |
|
| 138 |
# data_s<-obj_mod_interp[[1]][[k]]$data_s #object for the first date...20100103 |
|
| 139 |
# data_v<-obj_mod_interp[[1]][[k]]$data_v #object for the first date...20100103 |
|
| 140 |
# if (testing==TRUE & training==FALSE){
|
|
| 141 |
# return(data_v) |
|
| 142 |
# } |
|
| 143 |
# if (training==TRUE & testing==FALSE){
|
|
| 144 |
# return(data_s) |
|
| 145 |
# } |
|
| 146 |
# if (training==TRUE & testing==TRUE ){
|
|
| 147 |
# dataset_stat<-list(data_v,data_s) |
|
| 148 |
# names(dataset_stat)<-c("testing","training")
|
|
| 149 |
# return(dataset_stat) |
|
| 150 |
# } |
|
| 151 |
# } |
|
| 152 |
|
|
| 153 |
######### |
|
| 30 | 154 |
#loading R objects that might have similar names |
| 31 | 155 |
|
| 32 | 156 |
out_prefix<-"_method_comp7_12042012_" |
| 33 |
dates_input<-c("20100103","20100901")
|
|
| 157 |
infile2<-"list_365_dates_04212012.txt" |
|
| 158 |
|
|
| 34 | 159 |
i=2 |
| 35 | 160 |
##### LOAD USEFUL DATA |
| 36 | 161 |
|
| ... | ... | |
| 44 | 169 |
obj_mod_fus_name<-"results_mod_obj__365d_GAM_fusion_const_all_lstd_11022012.RData" |
| 45 | 170 |
obj_mod_cai_name<-"results_mod_obj__365d_GAM_CAI2_const_all_10312012.RData" |
| 46 | 171 |
|
| 47 |
gam_fus<-load_obj(file.path(path_data_fus,obj_mod_fus_name)) |
|
| 48 |
gam_cai<-load_obj(file.path(path_data_cai,obj_mod_cai_name)) #This contains all the info |
|
| 49 |
sampling_date_list<-gam_fus$sampling_obj$sampling_dat$date |
|
| 50 | 172 |
|
| 51 | 173 |
### Projection for the current region |
| 52 | 174 |
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"; |
| ... | ... | |
| 62 | 184 |
layerNames(s_raster)<-covar_names #Assigning names to the raster layers |
| 63 | 185 |
projection(s_raster)<-proj_str |
| 64 | 186 |
|
| 65 |
|
|
| 66 | 187 |
########## Load Snotel data |
| 67 | 188 |
infile_snotname<-"snot_OR_2010_sp2_methods_11012012_.shp" #load Snotel data |
| 68 | 189 |
snot_OR_2010_sp<-readOGR(".",sub(".shp","",infile_snotname))
|
| 69 | 190 |
snot_OR_2010_sp$date<-as.character(snot_OR_2010_sp$date) |
| 70 |
tmp_date<-snot_OR_2010_sp$date[1] |
|
| 71 |
#date<-strptime(dates[i], "%Y%m%d") # interpolation date being processed |
|
| 72 | 191 |
|
| 73 |
#change format of dates... |
|
| 74 |
#date_test<-strptime(tmp_date, "%e%m%y") # interpolation date being processed |
|
| 75 |
#date_test<-strptime(tmp_date, "%D") # interpolation date being processed |
|
| 76 |
#month<-strftime(date, "%m") # current month of the date being processed |
|
| 77 |
#LST_month<-paste("mm_",month,sep="") # name of LST month to be matched
|
|
| 192 |
#dates<-c("20100103","20100901")
|
|
| 193 |
#dates_snot<-c("10310","90110")
|
|
| 194 |
#dates<-c("20100101","20100103","20100301","20100302","20100501","20100502","20100801","20100802","20100901","20100902")
|
|
| 195 |
#dates_snot<-c("10110","10310","30110","30210","50110","50210","80110","80210","90110","90210")
|
|
| 78 | 196 |
|
| 197 |
#Use file with date |
|
| 198 |
dates<-readLines(file.path(path,infile2)) |
|
| 199 |
#Or use list of date in string |
|
| 200 |
#dates<-c("20100103","20100901")
|
|
| 79 | 201 |
|
| 202 |
dates_snot_tmp<-snot_OR_2010_sp$date |
|
| 203 |
dates_snot_formatted<-format_padding_month(dates_snot_tmp) |
|
| 204 |
date_test<-strptime(dates_snot_formatted, "%m%d%y") # interpolation date being processed |
|
| 205 |
snot_OR_2010_sp$date_formatted<-date_test |
|
| 80 | 206 |
#Load GHCN data used in modeling: training and validation site |
| 81 | 207 |
|
| 82 | 208 |
### load specific date...and plot: make a function to extract the diff and prediction... |
| 83 | 209 |
rast_diff_fc<-rast_fus_pred-rast_cai_pred |
| 84 | 210 |
layerNames(rast_diff)<-paste("diff",date_selected,sep="_")
|
| 85 | 211 |
|
| 86 |
|
|
| 87 | 212 |
####COMPARE WITH LOCATION OF GHCN and SNOTEL NETWORK |
| 88 | 213 |
|
| 89 |
X11(width=12,height=12) |
|
| 214 |
|
|
| 215 |
i=1 |
|
| 216 |
date_selected<-dates[i] |
|
| 217 |
|
|
| 218 |
X11(width=16,height=9) |
|
| 219 |
par(mfrow=c(1,2)) |
|
| 220 |
|
|
| 90 | 221 |
plot(rast_diff_fc) |
| 91 | 222 |
plot(snot_OR_2010_sp,pch=2,col="red",add=T) |
| 92 | 223 |
plot(data_stat,add=T) #This is the GHCN network |
| ... | ... | |
| 94 | 225 |
cex=0.8, col=c("red","black"),
|
| 95 | 226 |
pch=c(2,1)) |
| 96 | 227 |
title(paste("SNOTEL and GHCN networks on ", date_selected, sep=""))
|
| 97 |
savePlot(paste("fig1a_map_SNOT_GHCN_network_diff_bckgd",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 98 | 228 |
|
| 99 | 229 |
plot(ELEV_SRTM) |
| 100 | 230 |
plot(snot_OR_2010_sp,pch=2,col="red",add=T) |
| 101 | 231 |
plot(data_stat,add=T) |
| 102 | 232 |
legend("bottom",legend=c("SNOTEL", "GHCN"),
|
| 103 |
cex=0.8, col=c("red","black"),
|
|
| 104 |
pch=c(2,1))
|
|
| 105 |
title(paste("SNOTEL and GHCN networks on ", date_selected, sep=""))
|
|
| 106 |
savePlot(paste("fig1b_map_SNOT_GHCN_network_elev_bckgd",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 233 |
cex=0.8, col=c("red","black"),
|
|
| 234 |
pch=c(2,1)) |
|
| 235 |
title(paste("SNOTEL and GHCN networks", sep=""))
|
|
| 236 |
savePlot(paste("fig1_map_SNOT_GHCN_network_diff_elev_bckgd",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 107 | 237 |
|
| 238 |
#add histogram of elev for SNOT and GHCN |
|
| 239 |
hist(snot_data_selected$ELEV_SRTM,main="") |
|
| 240 |
title(paste("SNOT stations and Elevation",date_selected,sep=" "))
|
|
| 241 |
hist(data_vc$ELEV_SRTM,main="") |
|
| 242 |
title(paste("GHCN stations and Elevation",date_selected,sep=" "))
|
|
| 243 |
savePlot(paste("fig2_hist_elev_SNOT_GHCN_",out_prefix,".png", sep=""), type="png")
|
|
| 244 |
dev.off() |
|
| 108 | 245 |
## Select date from SNOT |
| 109 | 246 |
#not_selected<-subset(snot_OR_2010_sp, date=="90110" ) |
| 110 |
|
|
| 111 |
dates<-c("20100103","20100901")
|
|
| 112 |
dates_snot<-c("10310","90110")
|
|
| 113 |
i=2 |
|
| 114 |
|
|
| 115 |
for(i in 1:length(dates)){
|
|
| 116 |
|
|
| 247 |
list_ac_tab <-vector("list", length(dates)) #storing the accuracy metric data.frame in a list...
|
|
| 248 |
names(list_ac_tab)<-paste("date",1:length(dates),sep="")
|
|
| 249 |
X11(width=16,height=9) |
|
| 250 |
par(mfrow=c(1,2)) |
|
| 251 |
#for(i in 1:length(dates)){
|
|
| 252 |
for(i in 163:length(dates)){
|
|
| 117 | 253 |
date_selected<-dates[i] |
| 118 |
date_selected_snot<-as.integer(dates_snot[i]) #Change format of date at a later stage... |
|
| 119 |
snot_selected <-snot_OR_2010_sp[snot_OR_2010_sp$date==date_selected_snot,] |
|
| 120 |
#snot_selected<-na.omit(as.data.frame(snot_OR_2010_sp[snot_OR_2010_sp$date==90110,])) |
|
| 121 | 254 |
|
| 255 |
## Get the relevant raster layers with prediction for fusion and CAI |
|
| 256 |
oldpath<-getwd() |
|
| 257 |
setwd(path_data_cai) |
|
| 258 |
file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_CAI2_const_all_10312012.rst",sep="")) #Search for files in relation to fusion
|
|
| 259 |
lf_cai2c<-list.files(pattern=file_pat) #Search for files in relation to fusion |
|
| 260 |
rast_cai2c<-stack(lf_cai2c) #lf_cai2c CAI results with constant sampling over 365 dates |
|
| 261 |
rast_cai2c<-mask(rast_cai2c,mask_ELEV_SRTM) |
|
| 262 |
|
|
| 263 |
oldpath<-getwd() |
|
| 264 |
setwd(path_data_fus) |
|
| 265 |
file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_fusion_const_all_lstd_11022012.rst",sep="")) #Search for files in relation to fusion
|
|
| 266 |
lf_fus1c<-list.files(pattern=file_pat) #Search for files in relation to fusion |
|
| 267 |
rast_fus1c<-stack(lf_fus1c) |
|
| 268 |
rast_fus1c<-mask(rast_fus1c,mask_ELEV_SRTM) |
|
| 269 |
|
|
| 270 |
#PLOT ALL MODELS |
|
| 271 |
#Prepare for plotting |
|
| 272 |
|
|
| 273 |
setwd(path) #set path to the output path |
|
| 274 |
|
|
| 275 |
rast_fus_pred<-raster(rast_fus1c,1) # Select the first model from the stack i.e fusion with kriging for both steps |
|
| 276 |
rast_cai_pred<-raster(rast_cai2c,1) |
|
| 277 |
layerNames(rast_cai_pred)<-paste("cai",date_selected,sep="_")
|
|
| 278 |
layerNames(rast_fus_pred)<-paste("fus",date_selected,sep="_")
|
|
| 279 |
rast_pred2<-stack(rast_fus_pred,rast_cai_pred) |
|
| 280 |
#function to extract training and test from object from object models created earlier during interpolation... |
|
| 281 |
|
|
| 282 |
#load training and testing date for the specified date for fusion and CAI |
|
| 283 |
data_vf<-station_data_interp(date_selected,file.path(path_data_fus,obj_mod_fus_name),training=FALSE,testing=TRUE) |
|
| 284 |
#data_sf<-station_data_interp(date_selected,file.path(path_data_fus,obj_mod_fus_name),training=TRUE,testing=FALSE) |
|
| 285 |
data_vc<-station_data_interp(date_selected,file.path(path_data_cai,obj_mod_cai_name),training=FALSE,testing=TRUE) |
|
| 286 |
#data_sc<-station_data_interp(date_selected,file.path(path_data_cai,obj_mod_cai_name),training=TRUE,testing=FALSE) |
|
| 287 |
|
|
| 288 |
date_selected_snot<-strptime(date_selected,"%Y%m%d") |
|
| 289 |
snot_selected <-snot_OR_2010_sp[snot_OR_2010_sp$date_formatted==date_selected_snot,] |
|
| 290 |
#snot_selected<-na.omit(as.data.frame(snot_OR_2010_sp[snot_OR_2010_sp$date==90110,])) |
|
| 291 |
rast_diff_fc<-rast_fus_pred-rast_cai_pred |
|
| 122 | 292 |
LC_stack<-stack(LC1,LC2,LC3,LC4,LC6,LC7) |
| 123 |
rast_pred3<-stack(rast_diff_fc,rast_pred2,LC_stack) |
|
| 293 |
rast_pred3<-stack(rast_diff_fc,rast_pred2,ELEV_SRTM,LC_stack)
|
|
| 124 | 294 |
layerNames(rast_pred3)<-c("diff_fc","fus","CAI","ELEV_SRTM","LC1","LC2","LC3","LC4","LC6","LC7") #extract amount of veg...
|
| 125 | 295 |
|
| 126 |
#extract info |
|
| 127 |
extract_snot<-extract(rast_pred3,snot_selected) # |
|
| 128 |
snot_data_selected<-cbind(as.data.frame(snot_selected),extract_snot) |
|
| 129 |
snot_data_selected$res_f<-snot_data_selected$fus-snot_data_selected$tmax |
|
| 130 |
snot_data_selected$res_c<-snot_data_selected$CAI-snot_data_selected$tmax |
|
| 131 |
snot_data_selected<-(na.omit(as.data.frame(snot_data_selected))) |
|
| 132 |
|
|
| 133 |
#Plot predicted vs observed |
|
| 134 |
y_range<-range(c(data_vf$pred_mod7,snot_data_selected$fus),na.rm=T) |
|
| 135 |
x_range<-range(c(data_vf$dailyTmax,snot_data_selected$tmax),na.rm=T) |
|
| 136 |
plot(data_vf$pred_mod7,data_vf$dailyTmax, ylab="Observed daily tmax (C)", xlab="Predicted daily tmax (C)", |
|
| 296 |
#extract predicted tmax corresponding to |
|
| 297 |
extract_snot<-extract(rast_pred3,snot_selected) #return value from extract is a matrix (with input SPDF) |
|
| 298 |
snot_data_selected<-cbind(as.data.frame(snot_selected),extract_snot) #bind data together |
|
| 299 |
snot_data_selected$res_f<-snot_data_selected$fus-snot_data_selected$tmax #calculate the residuals for Fusion |
|
| 300 |
snot_data_selected$res_c<-snot_data_selected$CAI-snot_data_selected$tmax #calculate the residuals for CAI |
|
| 301 |
#snot_data_selected<-(na.omit(as.data.frame(snot_data_selected))) #remove rows containing NA, this may need to be modified later. |
|
| 302 |
|
|
| 303 |
###fig3: Plot predicted vs observed tmax |
|
| 304 |
#fig3a: FUS |
|
| 305 |
x_range<-range(c(data_vf$pred_mod7,snot_data_selected$fus,data_vc$pred_mod9,snot_data_selected$CAI),na.rm=T) |
|
| 306 |
y_range<-range(c(data_vf$dailyTmax,snot_data_selected$tmax,data_vc$dailyTmax,snot_data_selected$tmax),na.rm=T) |
|
| 307 |
plot(data_vf$pred_mod7,data_vf$dailyTmax, ylab="Observed daily tmax (C)", xlab="Fusion predicted daily tmax (C)", |
|
| 137 | 308 |
ylim=y_range,xlim=x_range) |
| 138 |
text(data_vf$pred_mod7,data_vf$dailyTmax,labels=data_vf$idx,pos=3) |
|
| 309 |
#text(data_vf$pred_mod7,data_vf$dailyTmax,labels=data_vf$idx,pos=3)
|
|
| 139 | 310 |
abline(0,1) #takes intercept at 0 and slope as 1 so display 1:1 ine |
| 140 | 311 |
grid(lwd=0.5,col="black") |
| 141 | 312 |
points(snot_data_selected$fus,snot_data_selected$tmax,pch=2,co="red") |
| 142 | 313 |
title(paste("Testing stations tmax fusion vs daily tmax",date_selected,sep=" "))
|
| 143 |
savePlot(paste("fig2_testing_scatterplot_pred_fus_observed_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 144 |
|
|
| 145 |
###### CAI |
|
| 146 |
y_range<-range(c(data_vc$pred_mod9,snot_data_selected$CAI),na.rm=T) |
|
| 147 |
x_range<-range(c(data_vc$dailyTmax,snot_data_selected$tmax),na.rm=T) |
|
| 148 |
plot(data_vc$pred_mod9,data_vc$dailyTmax, ylab="Observed daily tmax (C)", xlab="Predicted daily tmax (C)", |
|
| 314 |
legend("topleft",legend=c("GHCN", "SNOT"),
|
|
| 315 |
cex=1.2, col=c("black","red"),
|
|
| 316 |
pch=c(1,2)) |
|
| 317 |
#fig 3b: CAI |
|
| 318 |
#x_range<-range(c(data_vc$pred_mod9,snot_data_selected$CAI)) |
|
| 319 |
#y_range<-range(c(data_vc$dailyTmax,snot_data_selected$tmax)) |
|
| 320 |
plot(data_vc$pred_mod9,data_vc$dailyTmax, ylab="Observed daily tmax (C)", xlab="CAI predicted daily tmax (C)", |
|
| 149 | 321 |
ylim=y_range,xlim=x_range) |
| 150 |
text(data_vc$pred_mod9,data_vc$dailyTmax,labels=data_vf$idx,pos=3) |
|
| 322 |
#text(data_vc$pred_mod9,data_vc$dailyTmax,labels=data_vf$idx,pos=3)
|
|
| 151 | 323 |
abline(0,1) #takes intercept at 0 and slope as 1 so display 1:1 ine |
| 152 | 324 |
grid(lwd=0.5,col="black") |
| 153 |
points(snot_data_selected$CAI,snot_data_selected$tmax,pch=2,co="red") |
|
| 154 |
title(paste("Testing stations tmax CAI vs daily tmax",date_selected,sep=" "))
|
|
| 155 |
savePlot(paste("fig3_testing_scatterplot_pred_CAI_observed_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 156 |
|
|
| 157 |
##### ELEV CAI |
|
| 158 |
y_range<-range(c(data_vc$dailyTmax,snot_data_selected$tmax),na.rm=T) |
|
| 325 |
points(snot_data_selected$CAI,snot_data_selected$tmax,pch=2,co="red") |
|
| 326 |
#text(snot_data_selected$CAI,snot_data_selected$tmax,labels=1:nrow(snot_data_selected),pos=3) |
|
| 327 |
#title(paste("Testing stations tmax CAI vs daily tmax",date_selected,sep=" "))
|
|
| 328 |
legend("topleft",legend=c("GHCN", "SNOT"),
|
|
| 329 |
cex=1.2, col=c("black","red"),
|
|
| 330 |
pch=c(1,2)) |
|
| 331 |
savePlot(paste("fig3_testing_scatterplot_pred_fus_CAI_observed_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 332 |
|
|
| 333 |
##### Fig4a: ELEV-CAI |
|
| 159 | 334 |
y_range<-range(c(data_vc$pred_mod9,snot_data_selected$CAI),na.rm=T) |
| 335 |
#y_range<-range(c(data_vc$pred_mod9,snot_data_selected$CAI),na.rm=T) |
|
| 160 | 336 |
x_range<-range(c(data_vc$ELEV_SRTM,snot_data_selected$ELEV_SRTM),na.rm=T) |
| 161 |
|
|
| 162 |
plot(data_vc$ELEV_SRTM,data_vc$pred_mod9,ylab="Observed daily tmax (C)", xlab="Predicted daily tmax (C)", |
|
| 337 |
lm_mod1<-lm(data_vc$pred_mod9~data_vc$ELEV_SRTM) |
|
| 338 |
lm_mod2<-lm(snot_data_selected$CAI~snot_data_selected$ELEV_SRTM) |
|
| 339 |
plot(data_vc$ELEV_SRTM,data_vc$pred_mod9,ylab="Observed daily tmax (C)", xlab="Elevation (m)", |
|
| 163 | 340 |
ylim=y_range,xlim=x_range) |
| 164 |
text(data_vc$ELEV_SRTM,data_vc$pred_mod9,labels=data_vc$idx,pos=3) |
|
| 165 |
abline(0,1) #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
| 341 |
#text(data_vc$ELEV_SRTM,data_vc$pred_mod9,labels=data_vc$idx,pos=3) |
|
| 342 |
abline(lm_mod1) #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
| 343 |
abline(lm_mod2,col="red") #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
| 166 | 344 |
grid(lwd=0.5, col="black") |
| 167 | 345 |
points(snot_data_selected$ELEV_SRTM,snot_data_selected$CAI,pch=2,co="red") |
| 168 |
title(paste("Testing stations tmax CAI vs daily tmax",date_selected,sep=" "))
|
|
| 169 |
savePlot(paste("fig4_testing_scatterplot_pred_CAI_observed_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 346 |
title(paste("Testing stations tmax CAI vs elevation",date_selected,sep=" "))
|
|
| 347 |
legend("topleft",legend=c("GHCN", "SNOT"),
|
|
| 348 |
cex=1.2, col=c("black","red"),
|
|
| 349 |
pch=c(1,2)) |
|
| 170 | 350 |
|
| 171 |
##### ELEV CAI |
|
| 172 |
y_range<-range(c(data_vc$dailyTmax,snot_data_selected$tmax),na.rm=T) |
|
| 173 |
#y_range<-range(c(data_vc$pred_mod9,snot_data_selected$CAI),na.rm=T) |
|
| 174 |
x_range<-range(c(data_vc$ELEV_SRTM,snot_data_selected$ELEV_SRTM),na.rm=T) |
|
| 351 |
#Fig4bELEV-FUS |
|
| 352 |
y_range<-range(c(data_vf$pred_mod7,snot_data_selected$fus),na.rm=T) |
|
| 353 |
x_range<-range(c(data_vf$ELEV_SRTM,snot_data_selected$ELEV_SRTM),na.rm=T) |
|
| 354 |
lm_mod1<-lm(data_vf$pred_mod7~data_vf$ELEV_SRTM) |
|
| 355 |
lm_mod2<-lm(snot_data_selected$fus~snot_data_selected$ELEV_SRTM) |
|
| 356 |
plot(data_vf$ELEV_SRTM,data_vf$pred_mod7,ylab="Observed daily tmax (C)", xlab="Elevation (m)", |
|
| 357 |
ylim=y_range,xlim=x_range) |
|
| 358 |
#text(data_vc$ELEV_SRTM,data_vc$pred_mod9,labels=data_vc$idx,pos=3) |
|
| 359 |
abline(lm_mod1) #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
| 360 |
abline(lm_mod2,col="red") #takes intercept at 0 and slope as 1 so display 1:1 ine |
|
| 361 |
grid(lwd=0.5, col="black") |
|
| 362 |
points(snot_data_selected$ELEV_SRTM,snot_data_selected$fus,pch=2,co="red") |
|
| 363 |
title(paste("Testing stations tmax vs elevation",date_selected,sep=" "))
|
|
| 364 |
legend("topleft",legend=c("GHCN", "SNOT"),
|
|
| 365 |
cex=1.2, col=c("black","red"),
|
|
| 366 |
pch=c(1,2)) |
|
| 367 |
savePlot(paste("fig4_testing_scatterplot_pred_fus_CIA_elev_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 368 |
|
|
| 369 |
############ ACCURACY METRICS AND RESIDUALS ############# |
|
| 370 |
|
|
| 371 |
#START FIG 5 |
|
| 372 |
#####Fig5a: CAI vs FUSION: difference by plotting on in terms of the other |
|
| 373 |
lm_mod<-lm(snot_data_selected$CAI~snot_data_selected$fus) |
|
| 374 |
y_range<-range(c(data_vc$pred_mod9,snot_data_selected$CAI),na.rm=T) |
|
| 375 |
x_range<-range(c(data_vf$pred_mod7,snot_data_selected$fus),na.rm=T) |
|
| 175 | 376 |
|
| 176 |
plot(data_vc$ELEV_SRTM,data_vc$dailyTmax,ylab="Observed daily tmax (C)", xlab="Elevation (m)",
|
|
| 377 |
plot(data_vf$pred_mod7,data_vc$pred_mod9,ylab="Predicted CAI daily tmax (C)", xlab="Predicted fusion daily tmax (C)",
|
|
| 177 | 378 |
ylim=y_range,xlim=x_range) |
| 178 |
text(data_vc$ELEV_SRTM,data_vc$dailyTmax,labels=data_vc$idx,pos=3) |
|
| 379 |
#text(data_vc$ELEV_SRTM,data_vc$dailyTmax,labels=data_vc$idx,pos=3)
|
|
| 179 | 380 |
abline(0,1) #takes intercept at 0 and slope as 1 so display 1:1 ine |
| 381 |
abline(lm_mod,col="red") |
|
| 180 | 382 |
grid(lwd=0.5, col="black") |
| 181 |
points(snot_data_selected$ELEV_SRTM,snot_data_selected$tmax,pch=2,co="red") |
|
| 182 |
title(paste("Testing stations tmax CAI vs daily tmax",date_selected,sep=" "))
|
|
| 183 |
savePlot(paste("fig4_testing_scatterplot_pred_CAI_observed_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 184 |
|
|
| 185 |
brks<-seq(0,100,10) |
|
| 186 |
lab_brks<-seq(1,10,1) |
|
| 187 |
lc1_rcstat<-cut(snot_data_selected$LC1,breaks=brks,labels=lab_brks,right=F) |
|
| 188 |
snot_data_selected$LC1_rec<-lc1_rcstat |
|
| 189 |
lc_rcstat<-cut(snot_data_selected$LC3,breaks=brks,labels=lab_brks,right=F) |
|
| 190 |
snot_data_selected$LC3_rec<-lc_rcstat |
|
| 191 |
lc_rcstat<-cut(snot_data_selected$LC2,breaks=brks,labels=lab_brks,right=F) |
|
| 192 |
snot_data_selected$LC2_rec<-lc_rcstat |
|
| 193 |
lc_rcstat<-cut(snot_data_selected$LC4,breaks=brks,labels=lab_brks,right=F) |
|
| 194 |
snot_data_selected$LC4_rec<-lc_rcstat |
|
| 195 |
lc_rcstat<-cut(snot_data_selected$LC6,breaks=brks,labels=lab_brks,right=F) |
|
| 196 |
snot_data_selected$LC6_rec<-lc_rcstat |
|
| 197 |
lc_rcstat<-cut(snot_data_selected$LC7,breaks=brks,labels=lab_brks,right=F) |
|
| 198 |
snot_data_selected$LC7_rec<-lc_rcstat |
|
| 199 |
|
|
| 200 |
tmp<-aggregate(diff_fc~LC1_rec,data=snot_data_selected,FUN=mean) |
|
| 201 |
plot(tmp, type="l") |
|
| 202 |
tmp<-aggregate(diff_fc~LC2_rec,data=snot_data_selected,FUN=mean) |
|
| 203 |
plot(tmp) |
|
| 204 |
tmp<-aggregate(diff_fc~LC3_rec,data=snot_data_selected,FUN=mean) |
|
| 205 |
plot(tmp) |
|
| 206 |
tmp<-aggregate(diff_fc~LC4_rec,data=snot_data_selected,FUN=mean) |
|
| 207 |
plot(tmp) |
|
| 208 |
tmp<-aggregate(diff_fc~LC6_rec,data=snot_data_selected,FUN=mean) |
|
| 209 |
plot(tmp) |
|
| 210 |
plot(snot_data_selected$LC2_rec,snot_data_selected$diff_fc) |
|
| 211 |
table(snot_data_selected$LC7_rec) |
|
| 212 |
scatterplot(diff_fc~LC1|LC1_rec) |
|
| 213 |
mod_diff_fc_LC2<-lm(diff_fc~LC_rec,data=snot_data_selected) |
|
| 214 |
lc_melt<-melt(snot_data_selected[c("diff_fc","LC1_rec","LC2_rec","LC3_rec","LC4_rec","LC6_rec")],
|
|
| 215 |
measure=c("diff_fc"),
|
|
| 216 |
id=c("LC1_rec","LC2_rec","LC3_rec","LC4_rec","LC6_rec"),
|
|
| 217 |
na.rm=F) |
|
| 218 |
lc_cast<-cast(lc_melt,value~variable,mean) |
|
| 219 |
tabf_resmod9_locs<-as.data.frame(tabf_cast[,,1]) |
|
| 220 |
|
|
| 221 |
### PLOT LAND COVER |
|
| 222 |
X11() |
|
| 223 |
plot(zones_stat$zones,zones_stat$LC1_forest,type="b",ylim=c(-4.5,4.5), |
|
| 224 |
ylab="difference between CAI and fusion",xlab="land cover percent class/10") |
|
| 225 |
lines(lab_brks,snot_data_selected,col="red",type="b") |
|
| 226 |
lines(zones_stat$zones,zones_stat[,4],col="blue",type="b") |
|
| 227 |
lines(zones_stat$zones,zones_stat[,5],col="darkgreen",type="b") |
|
| 228 |
lines(zones_stat$zones,zones_stat[,6],col="purple",type="b") |
|
| 229 |
legend("topleft",legend=c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban"),
|
|
| 230 |
cex=1.2, col=c("black","red","blue","darkgreen","purple"),
|
|
| 231 |
lty=1) |
|
| 232 |
title(paste("Prediction tmax difference and land cover ",sep=""))
|
|
| 233 |
|
|
| 234 |
savePlot(paste("fig6_diff_prediction_tmax_difference_land cover",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 235 |
dev.off() |
|
| 236 |
snot_data_selected$LC2_rec<-lc2_rcstat |
|
| 237 |
|
|
| 238 |
############ ACCURACY METRICS AND RESIDUALS #### |
|
| 239 |
#snot_data_selected<- |
|
| 240 |
calc_accuracy_metrics<-function(x,y){
|
|
| 241 |
residuals<-x-y |
|
| 242 |
mae<-mean(abs(residuals),na.rm=T) |
|
| 243 |
rmse<-sqrt(mean((residuals)^2,na.rm=T)) |
|
| 244 |
me<-mean(residuals,na.rm=T) |
|
| 245 |
#r2<-cor(x,y,na.rm=T) |
|
| 246 |
metrics_dat<-list(mae,rmse,me) |
|
| 247 |
names(metrics_dat)<-c("mae","rmse","me")
|
|
| 248 |
return(metrics_dat) |
|
| 249 |
} |
|
| 250 |
|
|
| 251 |
#change to tmax when fixed problem... |
|
| 252 |
ac_metrics_fus<-calc_accuracy_metrics(snot_data_selected$tmax,snot_data_selected$fus) |
|
| 253 |
ac_metrics_cai<-calc_accuracy_metrics(snot_data_selected$tmax,snot_data_selected$CAI) |
|
| 254 |
|
|
| 255 |
#print(ac_metrics_fus,ac_metrics_cai) |
|
| 256 |
ac_metrics_fus |
|
| 257 |
ac_metrics_cai |
|
| 258 |
|
|
| 259 |
plot(snot_data_selected$E_SRTM,snot_data_selected$diff_fc) |
|
| 260 |
|
|
| 261 |
#DO MAE IN TERM OF ELEVATION CLASSES and LAND COVER CLASSES as well as diff... |
|
| 383 |
points(snot_data_selected$fus,snot_data_selected$CAI,pch=2,co="red") |
|
| 384 |
title(paste("Testing stations predicted tmax fusion vs CAI tmax",date_selected,sep=" "))
|
|
| 385 |
legend("topleft",legend=c("GHCN", "SNOT"),
|
|
| 386 |
cex=1.2, col=c("black","red"),
|
|
| 387 |
pch=c(1,2)) |
|
| 388 |
####Fig5b: diff vs elev: difference by plotting on in terms of elev |
|
| 389 |
diff_fc<-data_vf$pred_mod7-data_vc$pred_mod9 |
|
| 390 |
plot(snot_data_selected$ELEV_SRTM,snot_data_selected$diff_fc,pch=2,col="red") |
|
| 391 |
lm_mod<-lm(snot_data_selected$diff_fc~snot_data_selected$ELEV_SRTM) |
|
| 392 |
abline(lm_mod,col="red") |
|
| 393 |
points(data_vf$ELEV_SRTM,diff_fc) |
|
| 394 |
lm_mod<-lm(diff_fc~data_vf$ELEV_SRTM) |
|
| 395 |
abline(lm_mod) |
|
| 396 |
legend("topleft",legend=c("GHCN", "SNOT"),
|
|
| 397 |
cex=1.2, col=c("black","red"),
|
|
| 398 |
pch=c(1,2)) |
|
| 399 |
title(paste("Prediction tmax difference and elevation ",sep=""))
|
|
| 400 |
savePlot(paste("fig5_testing_scatterplot_pred_fus_CAI_observed_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 401 |
|
|
| 262 | 402 |
#DO diff IN TERM OF ELEVATION CLASSES as well as diff.. |
| 263 |
|
|
| 403 |
|
|
| 404 |
#### START FIG 6: difference fc vs elev |
|
| 405 |
#fig6a |
|
| 264 | 406 |
brks<-c(0,500,1000,1500,2000,2500,4000) |
| 265 | 407 |
lab_brks<-1:6 |
| 266 | 408 |
elev_rcstat<-cut(snot_data_selected$ELEV_SRTM,breaks=brks,labels=lab_brks,right=F) |
| 267 | 409 |
snot_data_selected$elev_rec<-elev_rcstat |
| 268 |
y_range<-range(c(snot_data_selected$diff_fc)) |
|
| 269 |
x_range<-range(c(elev_rcstat)) |
|
| 410 |
y_range<-range(c(snot_data_selected$diff_fc),na.rm=T)
|
|
| 411 |
x_range<-range(c(elev_rcstat),na.rm=T)
|
|
| 270 | 412 |
plot(elev_rcstat,snot_data_selected$diff_fc, ylab="diff_fc", xlab="ELEV_SRTM (m) ", |
| 271 | 413 |
ylim=y_range, xlim=x_range) |
| 272 | 414 |
#text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3) |
| 273 | 415 |
grid(lwd=0.5,col="black") |
| 274 | 416 |
title(paste("SNOT stations diff f vs Elevation",date_selected,sep=" "))
|
| 275 | 417 |
|
| 276 |
brks<-c(0,500,1000,1500,2000,2500,4000) |
|
| 277 |
lab_brks<-1:6 |
|
| 278 |
elev_rcstat<-cut(snot_data_selected$E_SRTM,breaks=brks,labels=lab_brks,right=F) |
|
| 418 |
###With fewer classes...fig6b |
|
| 419 |
brks<-c(0,1000,2000,3000,4000) |
|
| 420 |
lab_brks<-1:4 |
|
| 421 |
elev_rcstat<-cut(snot_data_selected$ELEV_SRTM,breaks=brks,labels=lab_brks,right=F) |
|
| 422 |
snot_data_selected$elev_rec<-elev_rcstat |
|
| 423 |
y_range<-range(c(snot_data_selected$diff_fc),na.rm=T) |
|
| 424 |
x_range<-range(c(elev_rcstat),na.rm=T) |
|
| 425 |
plot(elev_rcstat,snot_data_selected$diff_fc, ylab="diff_fc", xlab="ELEV_SRTM (m) ", |
|
| 426 |
ylim=y_range, xlim=x_range) |
|
| 427 |
#text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3) |
|
| 428 |
grid(lwd=0.5,col="black") |
|
| 429 |
title(paste("SNOT stations diff f vs Elevation",date_selected,sep=" "))
|
|
| 430 |
savePlot(paste("fig6_elevation_classes_diff_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 279 | 431 |
|
| 280 |
y_range<-range(c(snot_data_selected$res_f),na.rm=T) |
|
| 281 |
x_range<-range(c(elev_rcstat)) |
|
| 432 |
#START FIG 7 with residuals |
|
| 433 |
#fig 7a |
|
| 434 |
brks<-c(0,1000,2000,3000,4000) |
|
| 435 |
lab_brks<-1:4 |
|
| 436 |
elev_rcstat<-cut(snot_data_selected$ELEV_SRTM,breaks=brks,labels=lab_brks,right=F) |
|
| 437 |
snot_data_selected$elev_rec<-elev_rcstat |
|
| 438 |
y_range<-range(c(snot_data_selected$res_f,snot_data_selected$res_c),na.rm=T) |
|
| 439 |
x_range<-range(c(elev_rcstat),na.rm=T) |
|
| 282 | 440 |
plot(elev_rcstat,snot_data_selected$res_f, ylab="res_f", xlab="ELEV_SRTM (m) ", |
| 283 | 441 |
ylim=y_range, xlim=x_range) |
| 284 | 442 |
#text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3) |
| 285 | 443 |
grid(lwd=0.5,col="black") |
| 286 | 444 |
title(paste("SNOT stations residuals fusion vs Elevation",date_selected,sep=" "))
|
| 287 |
|
|
| 288 |
brks<-c(0,500,1000,1500,2000,2500,4000) |
|
| 289 |
lab_brks<-1:6 |
|
| 445 |
#fig 7b |
|
| 290 | 446 |
elev_rcstat<-cut(snot_data_selected$ELEV_SRTM,breaks=brks,labels=lab_brks,right=F) |
| 291 |
snot_data_selected$elev_rec<-elev_rcstat |
|
| 292 |
y_range<-range(c(snot_data_selected$res_c),na.rm=T) |
|
| 447 |
y_range<-range(c(snot_data_selected$res_c,snot_data_selected$res_f),na.rm=T) |
|
| 293 | 448 |
x_range<-range(c(elev_rcstat)) |
| 294 | 449 |
plot(elev_rcstat,snot_data_selected$res_c, ylab="res_c", xlab="ELEV_SRTM (m) ", |
| 295 | 450 |
ylim=y_range, xlim=x_range) |
| 296 | 451 |
#text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3) |
| 297 | 452 |
grid(lwd=0.5,col="black") |
| 298 | 453 |
title(paste("SNOT stations residuals CAI vs Elevation",date_selected,sep=" "))
|
| 299 |
|
|
| 300 |
#ADD BOTH |
|
| 301 |
plot(elev_rcstat,snot_data_selected$res_c, ylab="res_c", xlab="ELEV_SRTM (m) ", |
|
| 302 |
ylim=y_range, xlim=x_range) |
|
| 303 |
|
|
| 304 |
#CORRELATION BETWEEN RESIDUALS FUS and CAI |
|
| 305 |
|
|
| 306 |
y_range<-range(c(snot_data_selected$res_f,snot_data_selected$res_c),na.rm=T) |
|
| 307 |
x_range<-range(c(snot_data_selected$res_c,snot_data_selected$res_f),na.rm=T) |
|
| 308 |
plot(snot_data_selected$res_f,snot_data_selected$res_c, ylab="res_c", xlab="res_f ", |
|
| 309 |
ylim=y_range, xlim=x_range) |
|
| 310 |
abline(0,1) |
|
| 311 |
# |
|
| 312 |
|
|
| 313 |
#CORRELATION BETWEEN PREDICTION FUS and CAI |
|
| 314 |
|
|
| 315 |
y_range<-range(c(snot_data_selected$fus),na.rm=T) |
|
| 316 |
x_range<-range(c(snot_data_selected$CAI),na.rm=T) |
|
| 317 |
plot(snot_data_selected$fus,snot_data_selected$CAI, ylab="CAI", xlab="FUS", |
|
| 318 |
ylim=y_range, xlim=x_range) |
|
| 319 |
#### |
|
| 320 |
mae<-function(residuals){
|
|
| 454 |
savePlot(paste("fig7_elevation_classes_residuals_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 455 |
|
|
| 456 |
####### COMPARE CAI FUSION USING SNOTEL DATA WITH ACCURACY METRICS############### |
|
| 457 |
################ RESIDUALS and MAE etc. ##################### |
|
| 458 |
|
|
| 459 |
### Run for full list of date? --365 |
|
| 460 |
ac_tab_snot_fus<-calc_accuracy_metrics(snot_data_selected$tmax,snot_data_selected$fus) |
|
| 461 |
ac_tab_snot_cai<-calc_accuracy_metrics(snot_data_selected$tmax,snot_data_selected$CAI) |
|
| 462 |
ac_tab_ghcn_fus<-calc_accuracy_metrics(data_vf$dailyTmax,data_vf$pred_mod7) |
|
| 463 |
ac_tab_ghcn_cai<-calc_accuracy_metrics(data_vc$dailyTmax,data_vc$pred_mod9) |
|
| 464 |
|
|
| 465 |
ac_tab<-do.call(rbind,list(ac_tab_snot_fus,ac_tab_snot_cai,ac_tab_ghcn_fus,ac_tab_ghcn_cai)) |
|
| 466 |
rownames(ac_tab)<-c("snot_fus","snot_cai","ghcn_fus","ghcn_cai")
|
|
| 467 |
ac_tab$date<-date_selected |
|
| 468 |
list_ac_tab[[i]]<-ac_tab #storing the accuracy metric data.frame in a list... |
|
| 469 |
#save(list_ac_tab,) |
|
| 470 |
save(list_ac_tab,file= paste("list_ac_tab_", date_selected,out_prefix,".RData",sep=""))
|
|
| 471 |
|
|
| 472 |
#FIG8: boxplot of residuals for methods (fus, cai) using SNOT and GHCN data |
|
| 473 |
#fig8a |
|
| 474 |
y_range<-range(c(snot_data_selected$res_f,snot_data_selected$res_c,data_vf$res_mod7,data_vc$res_mod9),na.rm=T) |
|
| 475 |
boxplot(snot_data_selected$res_f,snot_data_selected$res_c,names=c("FUS","CAI"),ylim=y_range,ylab="Residuals tmax degree C")
|
|
| 476 |
title(paste("Residuals for fusion and CAI methods for SNOT data ",date_selected,sep=" "))
|
|
| 477 |
#fig8b |
|
| 478 |
boxplot(data_vf$res_mod7,data_vc$res_mod9,names=c("FUS","CAI"),ylim=y_range,ylab="Residuals tmax degree C")
|
|
| 479 |
title(paste("Residuals for fusion and CAI methods for GHCN data ",date_selected,sep=" "))
|
|
| 480 |
savePlot(paste("fig8_residuals_boxplot_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 481 |
|
|
| 482 |
mae_fun<-function(residuals){
|
|
| 321 | 483 |
mean(abs(residuals),na.rm=T) |
| 322 | 484 |
} |
| 323 | 485 |
|
| 324 | 486 |
mean_diff_fc<-aggregate(diff_fc~elev_rec,data=snot_data_selected,mean) |
| 325 |
mean_mae_c<-aggregate(res_c~elev_rec,data=snot_data_selected,mae) |
|
| 326 |
mean_mae_f<-aggregate(res_c~elev_rec,data=snot_data_selected,mae) |
|
| 327 |
plot(mean_fus,type="o") |
|
| 328 |
plot(as.integer(as.character(mean_fus$elev_rec)),mean_fus$diff_fc,type="b") |
|
| 329 |
plot(as.integer(as.character(mean_fus$elev_rec)),mean_fus$diff_fc,type="b") |
|
| 330 |
hist(snot_data_selected$E_SRTM) |
|
| 331 |
hist(data_vf$ELEV_SRTM) |
|
| 332 |
|
|
| 487 |
mean_mae_c<-aggregate(res_c~elev_rec,data=snot_data_selected,mae_fun) |
|
| 488 |
mean_mae_f<-aggregate(res_f~elev_rec,data=snot_data_selected,mae_fun) |
|
| 489 |
|
|
| 490 |
####FIG 9: plot MAE for fusion and CAI as well as boxplots of both thechnique |
|
| 491 |
#fig 9a: boxplot of residuals for MAE and CAI |
|
| 492 |
height<-cbind(snot_data_selected$res_f,snot_data_selected$res_c) |
|
| 493 |
boxplot(height,names=c("FUS","CAI"),ylab="Residuals tmax degree C")
|
|
| 494 |
title(paste("Residuals for fusion and CAI methods for SNOT data ",date_selected,sep=" "))
|
|
| 495 |
#par(new=TRUE) |
|
| 496 |
#abline(h=ac_tab[1,1],col="red") |
|
| 497 |
points(1,ac_tab[1,1],pch=5,col="red") |
|
| 498 |
points(2,ac_tab[2,1],pch=5,col="black") |
|
| 499 |
legend("bottom",legend=c("FUS_MAE", "CAI_MAE"),
|
|
| 500 |
cex=0.8, col=c("red","black"),
|
|
| 501 |
pch=c(2,1)) |
|
| 502 |
#fig 9b: MAE per 3 elevation classes:0-1000,1000-2000,2000-3000,3000-4000 |
|
| 503 |
y_range<-c(0,max(c(mean_mae_c[,2],mean_mae_f[,2]),na.rm=T)) |
|
| 504 |
plot(1:3,mean_mae_c[,2],ylim=y_range,type="n",ylab="MAE in degree C",xlab="elevation classes") |
|
| 505 |
points(mean_mae_c,ylim=y_range) |
|
| 506 |
lines(1:3,mean_mae_c[,2],col="black") |
|
| 507 |
par(new=TRUE) # key: ask for new plot without erasing old |
|
| 508 |
points(mean_mae_f,ylim=y_range) |
|
| 509 |
lines(1:3,mean_mae_f[,2],col="red") |
|
| 510 |
legend("bottom",legend=c("FUS_MAE", "CAI_MAE"),
|
|
| 511 |
cex=0.8, col=c("red","black"),
|
|
| 512 |
pch=c(2,1)) |
|
| 513 |
title(paste("MAE per elevation classes for SNOT data ",date_selected,sep=" "))
|
|
| 514 |
savePlot(paste("fig9_residuals_boxplot_MAE_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 515 |
|
|
| 516 |
### LM MODELS for difference and elevation categories |
|
| 517 |
## Are the differences plotted on fig 9 significant?? |
|
| 333 | 518 |
diffelev_mod<-lm(diff_fc~elev_rec,data=snot_data_selected) |
| 334 | 519 |
summary(diffelev_mod) |
| 520 |
##LM MODEL MAE PER ELEVATION CLASS: residuals for CAI |
|
| 335 | 521 |
diffelev_mod<-lm(res_c~elev_rec,data=snot_data_selected) |
| 336 | 522 |
summary(diffelev_mod) |
| 337 |
diffelev_mod$fit |
|
| 338 |
table(snot_data_selected$elev_rec) #Number of observation per class |
|
| 339 |
max(snot_data_selected$E_STRM) |
|
| 340 |
|
|
| 341 |
avl<-c(0,10,1,10,20,2,20,30,3,30,40,4,40,50,5,50,60,6,60,70,7,70,80,8,80,90,9,90,100,10)#Note that category 1 does not include 0!! |
|
| 342 |
rclmat<-matrix(avl,ncol=3,byrow=TRUE) |
|
| 343 |
|
|
| 344 |
|
|
| 345 |
#### DO THIS FOR IMAGE STACK...DIFF and LAND COVER... |
|
| 346 |
dat_stack<-stack(rast_diff,rast_fus_pred,rast_cai_pred, ELEV_SRTM) |
|
| 347 |
dat_analysis<-as(dat_stack,"SpatialGridDataFrame") |
|
| 348 |
names(dat_analysis)<-c("diff_fc","pred_fus","pred_cai","E_SRTM")
|
|
| 349 |
brks<-c(0,500,1000,1500,2000,2500,4000) |
|
| 350 |
lab_brks<-1:6 |
|
| 351 |
elev_rcstat<-cut(dat_analysis$E_SRTM,breaks=brks,labels=lab_brks,right=F) |
|
| 352 |
dat_analysis$elev_rec<-elev_rcstat |
|
| 523 |
##LM MODEL MAE PER ELEVATION CLASS: residuals for Fusions |
|
| 524 |
diffelev_mod<-lm(res_f~elev_rec,data=snot_data_selected) |
|
| 525 |
summary(diffelev_mod) |
|
| 353 | 526 |
|
| 354 |
spplot(dat_analysis,"elev_rec") |
|
| 355 |
spplot(dat_analysis,"diff_fc") |
|
| 356 |
mean_diff_fc<-aggregate(diff_fc~elev_rec,data=dat_analysis,mean) |
|
| 357 |
table(dat_analysis$elev_rec) #Number of observation per class |
|
| 527 |
### LM MODELS for RESIDUALS BETWEEN CAI AND FUSION |
|
| 528 |
## Are the differences plotted on fig 9 significant?? |
|
| 529 |
## STORE THE p values...?? overall and per cat? |
|
| 358 | 530 |
|
| 359 |
diffelev_mod<-lm(diff_fc~elev_rec,data=dat_analysis) |
|
| 360 |
summary(diffelev_mod) |
|
| 361 |
mean_rec6_val<-0.65993+(-8.56327) |
|
| 362 |
mean_diff_fc |
|
| 531 |
#diffelev_mod<-lm(res_f~elev_rec,data=snot_data_selected) |
|
| 532 |
#table(snot_data_selected$elev_rec) #Number of observation per class |
|
| 533 |
#max(snot_data_selected$E_STRM) |
|
| 363 | 534 |
|
| 364 |
brks<-c(0,500,1000,1500,2000,2500,4000) |
|
| 365 |
lab_brks<-1:6 |
|
| 366 |
elev_rcstat<-cut(data_vf$ELEV_SRTM,breaks=brks,labels=lab_brks,right=F) |
|
| 367 |
y_range<-range(c(diff_fc)) |
|
| 368 |
x_range<-range(c(elev_rcstat)) |
|
| 369 |
plot(elev_rcstat,diff_fc, ylab="diff_cf", xlab="ELEV_SRTM (m) ", |
|
| 370 |
ylim=y_range, xlim=x_range) |
|
| 371 |
text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3) |
|
| 372 |
grid(lwd=0.5,col="black") |
|
| 373 |
title(paste("Testing stations residuals fusion vs Elevation",date_selected,sep=" "))
|
|
| 374 |
|
|
| 375 |
# Combine both training and testing |
|
| 376 |
pred_fus<-c(data_vf$pred_mod7,data_sf$pred_mod7) |
|
| 377 |
pred_cai<-c(data_vc$pred_mod9,data_sc$pred_mod9) |
|
| 378 |
elev_station<-c(data_vf$ELEV_SRTM,data_sf$ELEV_SRTM) |
|
| 379 |
diff_fc<-pred_fus-pred_cai |
|
| 380 |
|
|
| 381 |
elev_rcstat<-cut(elev_station,breaks=brks,labels=lab_brks,right=F) |
|
| 382 |
y_range<-range(diff_fc) |
|
| 383 |
x_range<-range(elev_station) |
|
| 384 |
plot(elev_station,diff_fc, ylab="diff_fc", xlab="ELEV_SRTM (m) ", |
|
| 385 |
ylim=y_range, xlim=x_range) |
|
| 386 |
text(elev_rcstat,diff_fc,labels=data_vf$idx,pos=3) |
|
| 387 |
grid(lwd=0.5,col="black") |
|
| 388 |
title(paste("Testing stations residuals fusion vs Elevation",date_selected,sep=" "))
|
|
| 535 |
#res |
|
| 389 | 536 |
|
| 537 |
############################################# |
|
| 390 | 538 |
#USING BOTH validation and training |
| 391 |
} |
|
| 539 |
#This part is exploratory.... |
|
| 540 |
################## EXAMINING RESIDUALS AND DIFFERENCES IN LAND COVER......############ |
|
| 541 |
###### |
|
| 542 |
|
|
| 543 |
#LC_names<-c("LC1_rec","LC2_rec","LC3_rec","LC4_rec","LC6_rec")
|
|
| 544 |
suf_name<-c("rec1")
|
|
| 545 |
sum_var<-c("diff_fc")
|
|
| 546 |
LC_names<-c("LC1","LC2","LC3","LC4","LC6")
|
|
| 547 |
brks<-c(-1,20,40,60,80,101) |
|
| 548 |
lab_brks<-seq(1,5,1) |
|
| 549 |
#var_name<-LC_names; suffix<-"rec1"; s_function<-"mean";df<-snot_data_selected;summary_var<-"diff_fc" |
|
| 550 |
#reclassify_df(snot_data_selected,LC_names,var_name,brks,lab_brks,suffix,summary_var) |
|
| 551 |
|
|
| 552 |
#Calculate mean per land cover percentage |
|
| 553 |
data_agg<-reclassify_df(snot_data_selected,LC_names,brks,lab_brks,suf_name,sum_var) |
|
| 554 |
data_lc<-data_agg[[1]] |
|
| 555 |
snot_data_selected<-data_agg[[2]] |
|
| 556 |
|
|
| 557 |
by_name<-"rec1" |
|
| 558 |
df_lc_diff_fc<-merge_multiple_df(data_lc,by_name) |
|
| 559 |
|
|
| 560 |
###### FIG10: PLOT LAND COVER |
|
| 561 |
zones_stat<-df_lc_diff_fc #first land cover |
|
| 562 |
#names(zones_stat)<-c("lab_brks","LC")
|
|
| 563 |
y_range<-range(as.vector(t(zones_stat[,-1])),na.rm=T) |
|
| 564 |
lab_brks_mid<-c(10,30,50,70,90) |
|
| 565 |
plot(lab_brks_mid,zones_stat[,2],type="b",ylim=y_range,col="black", lwd=2, |
|
| 566 |
ylab="difference between fusion and CAI",xlab="land cover percent classes") |
|
| 567 |
lines(lab_brks_mid,zones_stat[,3],col="red",type="b") |
|
| 568 |
lines(lab_brks_mid,zones_stat[,4],col="blue",type="b") |
|
| 569 |
lines(lab_brks_mid,zones_stat[,5],col="darkgreen",type="b") |
|
| 570 |
lines(lab_brks_mid,zones_stat[,6],col="purple",type="b") |
|
| 571 |
legend("topleft",legend=c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban"),
|
|
| 572 |
cex=1.2, col=c("black","red","blue","darkgreen","purple"),
|
|
| 573 |
lty=1,lwd=1.8) |
|
| 574 |
title(paste("Prediction tmax difference and land cover ",date_selected,sep=""))
|
|
| 575 |
|
|
| 576 |
###NOW USE RESIDUALS FOR FUSION |
|
| 577 |
sum_var<-"res_f" |
|
| 578 |
suf_name<-"rec2" |
|
| 579 |
data_agg2<-reclassify_df(snot_data_selected,LC_names,brks,lab_brks,suf_name,sum_var) |
|
| 580 |
data_resf_lc<-data_agg2[[1]] |
|
| 581 |
#snot_data_selected<-data_agg[[2]] |
|
| 582 |
|
|
| 583 |
by_name<-"rec2" |
|
| 584 |
df_lc_resf<-merge_multiple_df(data_resf_lc,by_name) |
|
| 585 |
|
|
| 586 |
zones_stat<-df_lc_resf #first land cover |
|
| 587 |
#names(zones_stat)<-c("lab_brks","LC")
|
|
| 588 |
lab_brks_mid<-c(10,30,50,70,90) |
|
| 589 |
plot(lab_brks_mid,zones_stat[,2],type="b",ylim=y_range,col="black",lwd=2, |
|
| 590 |
ylab="tmax residuals fusion ",xlab="land cover percent classes") |
|
| 591 |
lines(lab_brks_mid,zones_stat[,3],col="red",type="b") |
|
| 592 |
lines(lab_brks_mid,zones_stat[,4],col="blue",type="b") |
|
| 593 |
lines(lab_brks_mid,zones_stat[,5],col="darkgreen",type="b") |
|
| 594 |
lines(lab_brks_mid,zones_stat[,6],col="purple",type="b") |
|
| 595 |
legend("topleft",legend=c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban"),
|
|
| 596 |
cex=1.2, col=c("black","red","blue","darkgreen","purple"),
|
|
| 597 |
lty=1,lwd=1.2) |
|
| 598 |
title(paste("Prediction tmax residuals and land cover ",date_selected,sep=""))
|
|
| 599 |
savePlot(paste("fig10_diff_prediction_tmax_diff_res_f_land cover",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 600 |
|
|
| 601 |
#### FIGURE11: res_f and res_c per land cover |
|
| 602 |
|
|
| 603 |
sum_var<-"res_c" |
|
| 604 |
suf_name<-"rec3" |
|
| 605 |
data_agg3<-reclassify_df(snot_data_selected,LC_names,brks,lab_brks,suf_name,sum_var) |
|
| 606 |
data_resc_lc<-data_agg3[[1]] |
|
| 607 |
snot_data_selected<-data_agg3[[2]] |
|
| 608 |
|
|
| 609 |
by_name<-"rec3" |
|
| 610 |
df_lc_resc<-merge_multiple_df(data_resc_lc,by_name) |
|
| 611 |
|
|
| 612 |
zones_stat<-df_lc_resc #first land cover |
|
| 613 |
#names(zones_stat)<-c("lab_brks","LC")
|
|
| 614 |
y_range<-range(as.vector(t(zones_stat[,-1])),na.rm=T) |
|
| 615 |
lab_brks_mid<-c(10,30,50,70,90) |
|
| 616 |
plot(lab_brks_mid,zones_stat[,2],type="b",ylim=y_range,col="black",lwd=2, |
|
| 617 |
ylab="tmax residuals CAI",xlab="land cover percent classes") |
|
| 618 |
lines(lab_brks_mid,zones_stat[,3],col="red",type="b") |
|
| 619 |
lines(lab_brks_mid,zones_stat[,4],col="blue",type="b") |
|
| 620 |
lines(lab_brks_mid,zones_stat[,5],col="darkgreen",type="b") |
|
| 621 |
lines(lab_brks_mid,zones_stat[,6],col="purple",type="b") |
|
| 622 |
legend("topleft",legend=c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban"),
|
|
| 623 |
cex=1.2, col=c("black","red","blue","darkgreen","purple"),
|
|
| 624 |
lty=1,lwd=1.2) |
|
| 625 |
title(paste("Prediction tmax residuals CAI and land cover ",date_selected,sep=""))
|
|
| 626 |
|
|
| 627 |
#fig11b |
|
| 628 |
zones_stat<-df_lc_resf #first land cover |
|
| 629 |
#names(zones_stat)<-c("lab_brks","LC")
|
|
| 630 |
y_range<-range(as.vector(t(zones_stat[,-1])),na.rm=T) |
|
| 631 |
lab_brks_mid<-c(10,30,50,70,90) |
|
| 632 |
plot(lab_brks_mid,zones_stat[,2],type="b",ylim=y_range,col="black",lwd=2, |
|
| 633 |
ylab="tmax residuals fusion ",xlab="land cover percent classes") |
|
| 634 |
lines(lab_brks_mid,zones_stat[,3],col="red",type="b") |
|
| 635 |
lines(lab_brks_mid,zones_stat[,4],col="blue",type="b") |
|
| 636 |
lines(lab_brks_mid,zones_stat[,5],col="darkgreen",type="b") |
|
| 637 |
lines(lab_brks_mid,zones_stat[,6],col="purple",type="b") |
|
| 638 |
legend("topleft",legend=c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban"),
|
|
| 639 |
cex=1.2, col=c("black","red","blue","darkgreen","purple"),
|
|
| 640 |
lty=1,lwd=1.2) |
|
| 641 |
title(paste("Prediction tmax residuals and land cover ",date_selected,sep=""))
|
|
| 642 |
#savePlot(paste("fig10_diff_prediction_tmax_diff_res_f_land cover",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 643 |
savePlot(paste("fig11_prediction_tmax_res_f_res_c_land cover",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 644 |
|
|
| 645 |
} |
|
| 646 |
|
|
| 647 |
#Collect accuracy information for different dates |
|
| 648 |
ac_data_xdates<-do.call(rbind,list_ac_tab) |
|
| 649 |
|
|
| 650 |
ac_data_xdates$mod_id<-rownames(ac_data_xdates) |
|
| 651 |
|
|
| 652 |
tmp_rownames<-rownames(ac_data_xdates) |
|
| 653 |
rowstr<-strsplit(tmp_rownames,"\\.") |
|
| 654 |
for (i in 1:length(rowstr)){
|
|
| 655 |
ac_data_xdates$mod_id[i]<-rowstr[[i]][[2]] |
|
| 656 |
} |
|
| 657 |
##Now subset for each model... |
|
| 658 |
|
|
| 659 |
mod_names<-unique(ac_data_xdates$mod_id) |
|
| 660 |
for (i in 1:length(rowstr)){
|
|
| 661 |
data_ac<-subset(ac_data_xdates,mod_id==mod_names[i]) |
|
| 662 |
data_name<-paste("data_ac_",mod_names[i],sep="")
|
|
| 663 |
assign(data_name,data_ac) |
|
| 664 |
} |
|
| 665 |
|
|
| 666 |
X11(12,12) |
|
| 667 |
boxplot(data_ac_ghcn_fus$mae) |
|
| 668 |
boxplot(data_ac_snot_fus$mae) |
|
| 669 |
boxplot(data_ac_ghcn_cai$mae) |
|
| 670 |
boxplot(data_ac_snot_cai$mae) |
|
| 671 |
boxplot(data_ac_snot_fus$mae,data_ac_snot_cai$mae,names=c("fus","CAI"))
|
|
| 672 |
boxplot(data_ac_ghcn_fus$mae,data_ac_ghcn_cai$mae,names=c("fus","CAI"))
|
|
| 673 |
boxplot(data_ac_ghcn_fus$mae,data_ac_ghcn_cai$mae,data_ac_snot_fus$mae,data_ac_snot_cai$mae,names=c("fus_SNOT","CAI_SNOT","fus_GHCN","CAI_GHCN"))
|
|
| 674 |
savePlot(paste("fig12_prediction_tmax_MAE_boxplot_fus_CAI_GHCN_SNOT_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 675 |
|
|
| 676 |
boxplot(data_ac_ghcn_fus$rmse,data_ac_ghcn_cai$rmse,data_ac_snot_fus$rmse,data_ac_snot_cai$rmse,names=c("fus_SNOT","CAI_SNOT","fus_GHCN","CAI_GHCN"))
|
|
| 677 |
savePlot(paste("fig12_prediction_tmax_RMSE_boxplot_fus_CAI_GHCN_SNOT_",date_selected,out_prefix,".png", sep=""), type="png")
|
|
| 678 |
|
|
| 679 |
filename<-paste("accuracy_table_GHCN_SNOT_", date_selected,out_prefix,".RData",sep="")
|
|
| 680 |
save(ac_data_xdates,file=filename) |
|
| 681 |
|
|
| 682 |
mean(data_ac_snot_fus) |
|
| 683 |
mean(data_ac_snot_cai) |
|
| 684 |
mean(data_ac_ghcn_fus) |
|
| 685 |
mean(data_ac_ghcn_cai) |
|
| 686 |
|
|
| 687 |
### END OF CODE |
|
| 688 |
#Write a part to caculate MAE per date... |
|
| 689 |
#ac_table_metrics<-do.call(rbind,ac_tab_list) |
|
| 690 |
|
|
| 691 |
#Subset and present the average MAE and RMSE for the dataset... |
|
| 692 |
|
|
| 693 |
#calculate average per month, extract LST too...? |
|
| 694 |
|
|
| 695 |
#################################################################### |
|
| 696 |
#From this line on: code is exploratory... |
|
| 697 |
#################################################################### |
|
| 698 |
#### DO THIS FOR IMAGE STACK...DIFF and LAND COVER...#### RESIDUALS AND LAND COVER... |
|
| 699 |
# |
|
| 700 |
# dat_stack<-stack(rast_diff,rast_fus_pred,rast_cai_pred, ELEV_SRTM) |
|
| 701 |
# dat_analysis<-as(dat_stack,"SpatialGridDataFrame") |
|
| 702 |
# names(dat_analysis)<-c("diff_fc","pred_fus","pred_cai","E_SRTM")
|
|
| 703 |
# brks<-c(0,500,1000,1500,2000,2500,4000) |
|
| 704 |
# lab_brks<-1:6 |
|
| 705 |
# elev_rcstat<-cut(dat_analysis$E_SRTM,breaks=brks,labels=lab_brks,right=F) |
|
| 706 |
# dat_analysis$elev_rec<-elev_rcstat |
|
| 707 |
# |
|
| 708 |
# spplot(dat_analysis,"elev_rec") |
|
| 709 |
# spplot(dat_analysis,"diff_fc") |
|
| 710 |
# mean_diff_fc<-aggregate(diff_fc~elev_rec,data=dat_analysis,mean) |
|
| 711 |
# table(dat_analysis$elev_rec) #Number of observation per class |
|
| 712 |
# |
|
| 713 |
# diffelev_mod<-lm(diff_fc~elev_rec,data=dat_analysis) |
|
| 714 |
# summary(diffelev_mod) |
|
| 715 |
# mean_rec6_val<-0.65993+(-8.56327) |
|
| 716 |
# mean_diff_fc |
|
| 717 |
# |
|
| 718 |
# brks<-c(0,500,1000,1500,2000,2500,4000) |
|
| 719 |
# lab_brks<-1:6 |
|
| 720 |
# elev_rcstat<-cut(data_vf$ELEV_SRTM,breaks=brks,labels=lab_brks,right=F) |
|
| 721 |
# y_range<-range(c(diff_fc)) |
|
| 722 |
# x_range<-range(c(elev_rcstat)) |
|
| 723 |
# plot(elev_rcstat,diff_fc, ylab="diff_cf", xlab="ELEV_SRTM (m) ", |
|
| 724 |
# ylim=y_range, xlim=x_range) |
|
| 725 |
# text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3) |
|
| 726 |
# grid(lwd=0.5,col="black") |
|
| 727 |
# title(paste("Testing stations residuals fusion vs Elevation",date_selected,sep=" "))
|
|
| 728 |
# |
|
| 729 |
# # Combine both training and testing |
|
| 730 |
# pred_fus<-c(data_vf$pred_mod7,data_sf$pred_mod7) |
|
| 731 |
# pred_cai<-c(data_vc$pred_mod9,data_sc$pred_mod9) |
|
| 732 |
# elev_station<-c(data_vf$ELEV_SRTM,data_sf$ELEV_SRTM) |
|
| 733 |
# diff_fc<-pred_fus-pred_cai |
|
| 734 |
# |
|
| 735 |
# elev_rcstat<-cut(elev_station,breaks=brks,labels=lab_brks,right=F) |
|
| 736 |
# y_range<-range(diff_fc) |
|
| 737 |
# x_range<-range(elev_station) |
|
| 738 |
# plot(elev_station,diff_fc, ylab="diff_fc", xlab="ELEV_SRTM (m) ", |
|
| 739 |
# ylim=y_range, xlim=x_range) |
|
| 740 |
# text(elev_rcstat,diff_fc,labels=data_vf$idx,pos=3) |
|
| 741 |
# grid(lwd=0.5,col="black") |
|
| 742 |
# title(paste("Testing stations residuals fusion vs Elevation",date_selected,sep=" "))
|
|
| 743 |
# |
|
Also available in: Unified diff
Methods comp part7-task#491- SNOTEL and GHCN analyses update main script calling function