1 |
6c21df27
|
Benoit Parmentier
|
#for(i in 1:length(dates)){
|
2 |
|
|
accuracy_comp_CAI_fus_function <- function(i){
|
3 |
|
|
date_selected<-dates[i]
|
4 |
|
|
|
5 |
|
|
## Get the relevant raster layers with prediction for fusion and CAI
|
6 |
|
|
oldpath<-getwd()
|
7 |
|
|
setwd(path_data_cai)
|
8 |
|
|
file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_CAI2_const_all_10312012.rst",sep="")) #Search for files in relation to fusion
|
9 |
|
|
lf_cai2c<-list.files(pattern=file_pat) #Search for files in relation to fusion
|
10 |
|
|
rast_cai2c<-stack(lf_cai2c) #lf_cai2c CAI results with constant sampling over 365 dates
|
11 |
|
|
rast_cai2c<-mask(rast_cai2c,mask_ELEV_SRTM)
|
12 |
|
|
|
13 |
|
|
oldpath<-getwd()
|
14 |
|
|
setwd(path_data_fus)
|
15 |
|
|
file_pat<-glob2rx(paste("*tmax_predicted*",date_selected,"*_365d_GAM_fusion_const_all_lstd_11022012.rst",sep="")) #Search for files in relation to fusion
|
16 |
|
|
lf_fus1c<-list.files(pattern=file_pat) #Search for files in relation to fusion
|
17 |
|
|
rast_fus1c<-stack(lf_fus1c)
|
18 |
|
|
rast_fus1c<-mask(rast_fus1c,mask_ELEV_SRTM)
|
19 |
|
|
|
20 |
|
|
#PLOT ALL MODELS
|
21 |
|
|
#Prepare for plotting
|
22 |
|
|
|
23 |
|
|
setwd(path) #set path to the output path
|
24 |
|
|
|
25 |
|
|
rast_fus_pred<-raster(rast_fus1c,1) # Select the first model from the stack i.e fusion with kriging for both steps
|
26 |
|
|
rast_cai_pred<-raster(rast_cai2c,1)
|
27 |
|
|
layerNames(rast_cai_pred)<-paste("cai",date_selected,sep="_")
|
28 |
|
|
layerNames(rast_fus_pred)<-paste("fus",date_selected,sep="_")
|
29 |
|
|
rast_pred2<-stack(rast_fus_pred,rast_cai_pred)
|
30 |
|
|
#function to extract training and test from object from object models created earlier during interpolation...
|
31 |
|
|
|
32 |
|
|
#load training and testing date for the specified date for fusion and CAI
|
33 |
|
|
data_vf<-station_data_interp(date_selected,file.path(path_data_fus,obj_mod_fus_name),training=FALSE,testing=TRUE)
|
34 |
|
|
#data_sf<-station_data_interp(date_selected,file.path(path_data_fus,obj_mod_fus_name),training=TRUE,testing=FALSE)
|
35 |
|
|
data_vc<-station_data_interp(date_selected,file.path(path_data_cai,obj_mod_cai_name),training=FALSE,testing=TRUE)
|
36 |
|
|
#data_sc<-station_data_interp(date_selected,file.path(path_data_cai,obj_mod_cai_name),training=TRUE,testing=FALSE)
|
37 |
|
|
|
38 |
|
|
date_selected_snot<-strptime(date_selected,"%Y%m%d")
|
39 |
|
|
snot_selected <-snot_OR_2010_sp[snot_OR_2010_sp$date_formatted==date_selected_snot,]
|
40 |
|
|
#snot_selected<-na.omit(as.data.frame(snot_OR_2010_sp[snot_OR_2010_sp$date==90110,]))
|
41 |
|
|
rast_diff_fc<-rast_fus_pred-rast_cai_pred
|
42 |
|
|
LC_stack<-stack(LC1,LC2,LC3,LC4,LC6,LC7)
|
43 |
|
|
rast_pred3<-stack(rast_diff_fc,rast_pred2,ELEV_SRTM,LC_stack)
|
44 |
|
|
layerNames(rast_pred3)<-c("diff_fc","fus","CAI","ELEV_SRTM","LC1","LC2","LC3","LC4","LC6","LC7") #extract amount of veg...
|
45 |
|
|
|
46 |
|
|
#extract predicted tmax corresponding to
|
47 |
|
|
extract_snot<-extract(rast_pred3,snot_selected) #return value from extract is a matrix (with input SPDF)
|
48 |
|
|
snot_data_selected<-cbind(as.data.frame(snot_selected),extract_snot) #bind data together
|
49 |
|
|
snot_data_selected$res_f<-snot_data_selected$fus-snot_data_selected$tmax #calculate the residuals for Fusion
|
50 |
|
|
snot_data_selected$res_c<-snot_data_selected$CAI-snot_data_selected$tmax #calculate the residuals for CAI
|
51 |
|
|
#snot_data_selected<-(na.omit(as.data.frame(snot_data_selected))) #remove rows containing NA, this may need to be modified later.
|
52 |
|
|
|
53 |
|
|
###fig3: Plot predicted vs observed tmax
|
54 |
|
|
#fig3a: FUS
|
55 |
|
|
png(paste("fig3_testing_scatterplot_pred_fus_CAI_observed_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""))
|
56 |
|
|
par(mfrow=c(1,2))
|
57 |
|
|
x_range<-range(c(data_vf$pred_mod7,snot_data_selected$fus,data_vc$pred_mod9,snot_data_selected$CAI),na.rm=T)
|
58 |
|
|
y_range<-range(c(data_vf$dailyTmax,snot_data_selected$tmax,data_vc$dailyTmax,snot_data_selected$tmax),na.rm=T)
|
59 |
|
|
plot(data_vf$pred_mod7,data_vf$dailyTmax, ylab="Observed daily tmax (C)", xlab="Fusion predicted daily tmax (C)",
|
60 |
|
|
ylim=y_range,xlim=x_range)
|
61 |
|
|
#text(data_vf$pred_mod7,data_vf$dailyTmax,labels=data_vf$idx,pos=3)
|
62 |
|
|
abline(0,1) #takes intercept at 0 and slope as 1 so display 1:1 ine
|
63 |
|
|
grid(lwd=0.5,col="black")
|
64 |
|
|
points(snot_data_selected$fus,snot_data_selected$tmax,pch=2,co="red")
|
65 |
|
|
title(paste("Testing stations tmax fusion vs daily tmax",date_selected,sep=" "))
|
66 |
|
|
legend("topleft",legend=c("GHCN", "SNOT"),
|
67 |
|
|
cex=1.2, col=c("black","red"),
|
68 |
|
|
pch=c(1,2))
|
69 |
|
|
#fig 3b: CAI
|
70 |
|
|
#x_range<-range(c(data_vc$pred_mod9,snot_data_selected$CAI))
|
71 |
|
|
#y_range<-range(c(data_vc$dailyTmax,snot_data_selected$tmax))
|
72 |
|
|
plot(data_vc$pred_mod9,data_vc$dailyTmax, ylab="Observed daily tmax (C)", xlab="CAI predicted daily tmax (C)",
|
73 |
|
|
ylim=y_range,xlim=x_range)
|
74 |
|
|
#text(data_vc$pred_mod9,data_vc$dailyTmax,labels=data_vf$idx,pos=3)
|
75 |
|
|
abline(0,1) #takes intercept at 0 and slope as 1 so display 1:1 ine
|
76 |
|
|
grid(lwd=0.5,col="black")
|
77 |
|
|
points(snot_data_selected$CAI,snot_data_selected$tmax,pch=2,co="red")
|
78 |
|
|
#text(snot_data_selected$CAI,snot_data_selected$tmax,labels=1:nrow(snot_data_selected),pos=3)
|
79 |
|
|
#title(paste("Testing stations tmax CAI vs daily tmax",date_selected,sep=" "))
|
80 |
|
|
legend("topleft",legend=c("GHCN", "SNOT"),
|
81 |
|
|
cex=1.2, col=c("black","red"),
|
82 |
|
|
pch=c(1,2))
|
83 |
|
|
#savePlot(paste("fig3_testing_scatterplot_pred_fus_CAI_observed_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""), type="png")
|
84 |
|
|
dev.off()
|
85 |
|
|
|
86 |
|
|
##### Fig4a: ELEV-CAI
|
87 |
|
|
png(paste("fig4_testing_scatterplot_pred_fus_CIA_elev_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""))
|
88 |
|
|
par(mfrow=c(1,2))
|
89 |
|
|
y_range<-range(c(data_vc$pred_mod9,snot_data_selected$CAI),na.rm=T)
|
90 |
|
|
#y_range<-range(c(data_vc$pred_mod9,snot_data_selected$CAI),na.rm=T)
|
91 |
|
|
x_range<-range(c(data_vc$ELEV_SRTM,snot_data_selected$ELEV_SRTM),na.rm=T)
|
92 |
|
|
lm_mod1<-lm(data_vc$pred_mod9~data_vc$ELEV_SRTM)
|
93 |
|
|
lm_mod2<-lm(snot_data_selected$CAI~snot_data_selected$ELEV_SRTM)
|
94 |
|
|
plot(data_vc$ELEV_SRTM,data_vc$pred_mod9,ylab="Observed daily tmax (C)", xlab="Elevation (m)",
|
95 |
|
|
ylim=y_range,xlim=x_range)
|
96 |
|
|
#text(data_vc$ELEV_SRTM,data_vc$pred_mod9,labels=data_vc$idx,pos=3)
|
97 |
|
|
abline(lm_mod1) #takes intercept at 0 and slope as 1 so display 1:1 ine
|
98 |
|
|
abline(lm_mod2,col="red") #takes intercept at 0 and slope as 1 so display 1:1 ine
|
99 |
|
|
grid(lwd=0.5, col="black")
|
100 |
|
|
points(snot_data_selected$ELEV_SRTM,snot_data_selected$CAI,pch=2,co="red")
|
101 |
|
|
title(paste("Testing stations tmax CAI vs elevation",date_selected,sep=" "))
|
102 |
|
|
legend("topleft",legend=c("GHCN", "SNOT"),
|
103 |
|
|
cex=1.2, col=c("black","red"),
|
104 |
|
|
pch=c(1,2))
|
105 |
|
|
|
106 |
|
|
#Fig4bELEV-FUS
|
107 |
|
|
y_range<-range(c(data_vf$pred_mod7,snot_data_selected$fus),na.rm=T)
|
108 |
|
|
x_range<-range(c(data_vf$ELEV_SRTM,snot_data_selected$ELEV_SRTM),na.rm=T)
|
109 |
|
|
lm_mod1<-lm(data_vf$pred_mod7~data_vf$ELEV_SRTM)
|
110 |
|
|
lm_mod2<-lm(snot_data_selected$fus~snot_data_selected$ELEV_SRTM)
|
111 |
|
|
plot(data_vf$ELEV_SRTM,data_vf$pred_mod7,ylab="Observed daily tmax (C)", xlab="Elevation (m)",
|
112 |
|
|
ylim=y_range,xlim=x_range)
|
113 |
|
|
#text(data_vc$ELEV_SRTM,data_vc$pred_mod9,labels=data_vc$idx,pos=3)
|
114 |
|
|
abline(lm_mod1) #takes intercept at 0 and slope as 1 so display 1:1 ine
|
115 |
|
|
abline(lm_mod2,col="red") #takes intercept at 0 and slope as 1 so display 1:1 ine
|
116 |
|
|
grid(lwd=0.5, col="black")
|
117 |
|
|
points(snot_data_selected$ELEV_SRTM,snot_data_selected$fus,pch=2,co="red")
|
118 |
|
|
title(paste("Testing stations tmax vs elevation",date_selected,sep=" "))
|
119 |
|
|
legend("topleft",legend=c("GHCN", "SNOT"),
|
120 |
|
|
cex=1.2, col=c("black","red"),
|
121 |
|
|
pch=c(1,2))
|
122 |
|
|
#savePlot(paste("fig4_testing_scatterplot_pred_fus_CIA_elev_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""), type="png")
|
123 |
|
|
dev.off()
|
124 |
|
|
############ ACCURACY METRICS AND RESIDUALS #############
|
125 |
|
|
|
126 |
|
|
#START FIG 5
|
127 |
|
|
#####Fig5a: CAI vs FUSION: difference by plotting on in terms of the other
|
128 |
|
|
png(paste("fig5_testing_scatterplot_pred_fus_CAI_observed_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""))
|
129 |
|
|
par(mfrow=c(1,2))
|
130 |
|
|
lm_mod<-lm(snot_data_selected$CAI~snot_data_selected$fus)
|
131 |
|
|
y_range<-range(c(data_vc$pred_mod9,snot_data_selected$CAI),na.rm=T)
|
132 |
|
|
x_range<-range(c(data_vf$pred_mod7,snot_data_selected$fus),na.rm=T)
|
133 |
|
|
|
134 |
|
|
plot(data_vf$pred_mod7,data_vc$pred_mod9,ylab="Predicted CAI daily tmax (C)", xlab="Predicted fusion daily tmax (C)",
|
135 |
|
|
ylim=y_range,xlim=x_range)
|
136 |
|
|
#text(data_vc$ELEV_SRTM,data_vc$dailyTmax,labels=data_vc$idx,pos=3)
|
137 |
|
|
abline(0,1) #takes intercept at 0 and slope as 1 so display 1:1 ine
|
138 |
|
|
abline(lm_mod,col="red")
|
139 |
|
|
grid(lwd=0.5, col="black")
|
140 |
|
|
points(snot_data_selected$fus,snot_data_selected$CAI,pch=2,co="red")
|
141 |
|
|
title(paste("Testing stations predicted tmax fusion vs CAI tmax",date_selected,sep=" "))
|
142 |
|
|
legend("topleft",legend=c("GHCN", "SNOT"),
|
143 |
|
|
cex=1.2, col=c("black","red"),
|
144 |
|
|
pch=c(1,2))
|
145 |
|
|
####Fig5b: diff vs elev: difference by plotting on in terms of elev
|
146 |
|
|
diff_fc<-data_vf$pred_mod7-data_vc$pred_mod9
|
147 |
|
|
plot(snot_data_selected$ELEV_SRTM,snot_data_selected$diff_fc,pch=2,col="red")
|
148 |
|
|
lm_mod<-lm(snot_data_selected$diff_fc~snot_data_selected$ELEV_SRTM)
|
149 |
|
|
abline(lm_mod,col="red")
|
150 |
|
|
points(data_vf$ELEV_SRTM,diff_fc)
|
151 |
|
|
lm_mod<-lm(diff_fc~data_vf$ELEV_SRTM)
|
152 |
|
|
abline(lm_mod)
|
153 |
|
|
legend("topleft",legend=c("GHCN", "SNOT"),
|
154 |
|
|
cex=1.2, col=c("black","red"),
|
155 |
|
|
pch=c(1,2))
|
156 |
|
|
title(paste("Prediction tmax difference and elevation ",sep=""))
|
157 |
|
|
dev.off()
|
158 |
|
|
#savePlot(paste("fig5_testing_scatterplot_pred_fus_CAI_observed_SNOT_GHCN_",date_selected,out_prefix,".png", sep=""), type="png")
|
159 |
|
|
|
160 |
|
|
#DO diff IN TERM OF ELEVATION CLASSES as well as diff..
|
161 |
|
|
|
162 |
|
|
#### START FIG 6: difference fc vs elev
|
163 |
|
|
#fig6a
|
164 |
|
|
png(paste("fig6_elevation_classes_diff_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""))
|
165 |
|
|
par(mfrow=c(1,2))
|
166 |
|
|
brks<-c(0,500,1000,1500,2000,2500,4000)
|
167 |
|
|
lab_brks<-1:6
|
168 |
|
|
elev_rcstat<-cut(snot_data_selected$ELEV_SRTM,breaks=brks,labels=lab_brks,right=F)
|
169 |
|
|
snot_data_selected$elev_rec<-elev_rcstat
|
170 |
|
|
y_range<-range(c(snot_data_selected$diff_fc),na.rm=T)
|
171 |
|
|
x_range<-range(c(elev_rcstat),na.rm=T)
|
172 |
|
|
plot(elev_rcstat,snot_data_selected$diff_fc, ylab="diff_fc", xlab="ELEV_SRTM (m) ",
|
173 |
|
|
ylim=y_range, xlim=x_range)
|
174 |
|
|
#text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3)
|
175 |
|
|
grid(lwd=0.5,col="black")
|
176 |
|
|
title(paste("SNOT stations diff f vs Elevation",date_selected,sep=" "))
|
177 |
|
|
|
178 |
|
|
###With fewer classes...fig6b
|
179 |
|
|
brks<-c(0,1000,2000,3000,4000)
|
180 |
|
|
lab_brks<-1:4
|
181 |
|
|
elev_rcstat<-cut(snot_data_selected$ELEV_SRTM,breaks=brks,labels=lab_brks,right=F)
|
182 |
|
|
snot_data_selected$elev_rec<-elev_rcstat
|
183 |
|
|
y_range<-range(c(snot_data_selected$diff_fc),na.rm=T)
|
184 |
|
|
x_range<-range(c(elev_rcstat),na.rm=T)
|
185 |
|
|
plot(elev_rcstat,snot_data_selected$diff_fc, ylab="diff_fc", xlab="ELEV_SRTM (m) ",
|
186 |
|
|
ylim=y_range, xlim=x_range)
|
187 |
|
|
#text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3)
|
188 |
|
|
grid(lwd=0.5,col="black")
|
189 |
|
|
title(paste("SNOT stations diff f vs Elevation",date_selected,sep=" "))
|
190 |
|
|
#savePlot(paste("fig6_elevation_classes_diff_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""), type="png")
|
191 |
|
|
dev.off()
|
192 |
|
|
#START FIG 7 with residuals
|
193 |
|
|
#fig 7a
|
194 |
|
|
png(paste("fig7_elevation_classes_residuals_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""))
|
195 |
|
|
par(mfrow=c(1,2))
|
196 |
|
|
brks<-c(0,1000,2000,3000,4000)
|
197 |
|
|
lab_brks<-1:4
|
198 |
|
|
elev_rcstat<-cut(snot_data_selected$ELEV_SRTM,breaks=brks,labels=lab_brks,right=F)
|
199 |
|
|
snot_data_selected$elev_rec<-elev_rcstat
|
200 |
|
|
y_range<-range(c(snot_data_selected$res_f,snot_data_selected$res_c),na.rm=T)
|
201 |
|
|
x_range<-range(c(elev_rcstat),na.rm=T)
|
202 |
|
|
plot(elev_rcstat,snot_data_selected$res_f, ylab="res_f", xlab="ELEV_SRTM (m) ",
|
203 |
|
|
ylim=y_range, xlim=x_range)
|
204 |
|
|
#text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3)
|
205 |
|
|
grid(lwd=0.5,col="black")
|
206 |
|
|
title(paste("SNOT stations residuals fusion vs Elevation",date_selected,sep=" "))
|
207 |
|
|
#fig 7b
|
208 |
|
|
elev_rcstat<-cut(snot_data_selected$ELEV_SRTM,breaks=brks,labels=lab_brks,right=F)
|
209 |
|
|
y_range<-range(c(snot_data_selected$res_c,snot_data_selected$res_f),na.rm=T)
|
210 |
|
|
x_range<-range(c(elev_rcstat))
|
211 |
|
|
plot(elev_rcstat,snot_data_selected$res_c, ylab="res_c", xlab="ELEV_SRTM (m) ",
|
212 |
|
|
ylim=y_range, xlim=x_range)
|
213 |
|
|
#text(elev_rcstat,diff_cf,labels=data_vf$idx,pos=3)
|
214 |
|
|
grid(lwd=0.5,col="black")
|
215 |
|
|
title(paste("SNOT stations residuals CAI vs Elevation",date_selected,sep=" "))
|
216 |
|
|
#savePlot(paste("fig7_elevation_classes_residuals_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""), type="png")
|
217 |
|
|
dev.off()
|
218 |
|
|
####### COMPARE CAI FUSION USING SNOTEL DATA WITH ACCURACY METRICS###############
|
219 |
|
|
################ RESIDUALS and MAE etc. #####################
|
220 |
|
|
browser()
|
221 |
|
|
### Run for full list of date? --365
|
222 |
|
|
ac_tab_snot_fus<-calc_accuracy_metrics(snot_data_selected$tmax,snot_data_selected$fus)
|
223 |
|
|
ac_tab_snot_cai<-calc_accuracy_metrics(snot_data_selected$tmax,snot_data_selected$CAI)
|
224 |
|
|
ac_tab_ghcn_fus<-calc_accuracy_metrics(data_vf$dailyTmax,data_vf$pred_mod7)
|
225 |
|
|
ac_tab_ghcn_cai<-calc_accuracy_metrics(data_vc$dailyTmax,data_vc$pred_mod9)
|
226 |
|
|
|
227 |
|
|
ac_tab<-do.call(rbind,list(ac_tab_snot_fus,ac_tab_snot_cai,ac_tab_ghcn_fus,ac_tab_ghcn_cai))
|
228 |
|
|
ac_tab$mod_id<-c("snot_fus","snot_cai","ghcn_fus","ghcn_cai")
|
229 |
|
|
ac_tab$date<-date_selected
|
230 |
|
|
|
231 |
|
|
list_ac_tab[[i]]<-ac_tab #storing the accuracy metric data.frame in a list...
|
232 |
|
|
#save(list_ac_tab,)
|
233 |
|
|
save(list_ac_tab,file= paste("list_ac_tab_", date_selected,out_prefix,".RData",sep=""))
|
234 |
|
|
|
235 |
|
|
#FIG8: boxplot of residuals for methods (fus, cai) using SNOT and GHCN data
|
236 |
|
|
#fig8a
|
237 |
|
|
png(paste("fig8_residuals_boxplot_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""))
|
238 |
|
|
par(mfrow=c(1,2))
|
239 |
|
|
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)
|
240 |
|
|
boxplot(snot_data_selected$res_f,snot_data_selected$res_c,names=c("FUS","CAI"),ylim=y_range,ylab="Residuals tmax degree C")
|
241 |
|
|
title(paste("Residuals for fusion and CAI methods for SNOT data ",date_selected,sep=" "))
|
242 |
|
|
#fig8b
|
243 |
|
|
boxplot(data_vf$res_mod7,data_vc$res_mod9,names=c("FUS","CAI"),ylim=y_range,ylab="Residuals tmax degree C")
|
244 |
|
|
title(paste("Residuals for fusion and CAI methods for GHCN data ",date_selected,sep=" "))
|
245 |
|
|
#savePlot(paste("fig8_residuals_boxplot_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""), type="png")
|
246 |
|
|
dev.off()
|
247 |
|
|
mae_fun<-function(residuals){
|
248 |
|
|
mean(abs(residuals),na.rm=T)
|
249 |
|
|
}
|
250 |
|
|
|
251 |
|
|
mean_diff_fc<-aggregate(diff_fc~elev_rec,data=snot_data_selected,mean)
|
252 |
|
|
mean_mae_c<-aggregate(res_c~elev_rec,data=snot_data_selected,mae_fun)
|
253 |
|
|
mean_mae_f<-aggregate(res_f~elev_rec,data=snot_data_selected,mae_fun)
|
254 |
|
|
|
255 |
|
|
####FIG 9: plot MAE for fusion and CAI as well as boxplots of both thechnique
|
256 |
|
|
#fig 9a: boxplot of residuals for MAE and CAI
|
257 |
|
|
png(paste("fig9_residuals_boxplot_MAE_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""))
|
258 |
|
|
par(mfrow=c(1,2))
|
259 |
|
|
height<-cbind(snot_data_selected$res_f,snot_data_selected$res_c)
|
260 |
|
|
boxplot(height,names=c("FUS","CAI"),ylab="Residuals tmax degree C")
|
261 |
|
|
title(paste("Residuals for fusion and CAI methods for SNOT data ",date_selected,sep=" "))
|
262 |
|
|
#par(new=TRUE)
|
263 |
|
|
#abline(h=ac_tab[1,1],col="red")
|
264 |
|
|
points(1,ac_tab[1,1],pch=5,col="red")
|
265 |
|
|
points(2,ac_tab[2,1],pch=5,col="black")
|
266 |
|
|
legend("bottom",legend=c("FUS_MAE", "CAI_MAE"),
|
267 |
|
|
cex=0.8, col=c("red","black"),
|
268 |
|
|
pch=c(2,1))
|
269 |
|
|
#fig 9b: MAE per 3 elevation classes:0-1000,1000-2000,2000-3000,3000-4000
|
270 |
|
|
y_range<-c(0,max(c(mean_mae_c[,2],mean_mae_f[,2]),na.rm=T))
|
271 |
|
|
plot(1:3,mean_mae_c[,2],ylim=y_range,type="n",ylab="MAE in degree C",xlab="elevation classes")
|
272 |
|
|
points(mean_mae_c,ylim=y_range)
|
273 |
|
|
lines(1:3,mean_mae_c[,2],col="black")
|
274 |
|
|
par(new=TRUE) # key: ask for new plot without erasing old
|
275 |
|
|
points(mean_mae_f,ylim=y_range)
|
276 |
|
|
lines(1:3,mean_mae_f[,2],col="red")
|
277 |
|
|
legend("bottom",legend=c("FUS_MAE", "CAI_MAE"),
|
278 |
|
|
cex=0.8, col=c("red","black"),
|
279 |
|
|
pch=c(2,1))
|
280 |
|
|
title(paste("MAE per elevation classes for SNOT data ",date_selected,sep=" "))
|
281 |
|
|
#savePlot(paste("fig9_residuals_boxplot_MAE_SNOT_GHCN_network",date_selected,out_prefix,".png", sep=""), type="png")
|
282 |
|
|
dev.off()
|
283 |
|
|
### LM MODELS for difference and elevation categories
|
284 |
|
|
## Are the differences plotted on fig 9 significant??
|
285 |
|
|
diffelev_mod<-lm(diff_fc~elev_rec,data=snot_data_selected)
|
286 |
|
|
summary(diffelev_mod)
|
287 |
|
|
##LM MODEL MAE PER ELEVATION CLASS: residuals for CAI
|
288 |
|
|
diffelev_mod<-lm(res_c~elev_rec,data=snot_data_selected)
|
289 |
|
|
summary(diffelev_mod)
|
290 |
|
|
##LM MODEL MAE PER ELEVATION CLASS: residuals for Fusions
|
291 |
|
|
diffelev_mod<-lm(res_f~elev_rec,data=snot_data_selected)
|
292 |
|
|
summary(diffelev_mod)
|
293 |
|
|
|
294 |
|
|
### LM MODELS for RESIDUALS BETWEEN CAI AND FUSION
|
295 |
|
|
## Are the differences plotted on fig 9 significant??
|
296 |
|
|
## STORE THE p values...?? overall and per cat?
|
297 |
|
|
|
298 |
|
|
#diffelev_mod<-lm(res_f~elev_rec,data=snot_data_selected)
|
299 |
|
|
#table(snot_data_selected$elev_rec) #Number of observation per class
|
300 |
|
|
#max(snot_data_selected$E_STRM)
|
301 |
|
|
|
302 |
|
|
#res
|
303 |
|
|
|
304 |
|
|
#############################################
|
305 |
|
|
#USING BOTH validation and training
|
306 |
|
|
#This part is exploratory....
|
307 |
|
|
################## EXAMINING RESIDUALS AND DIFFERENCES IN LAND COVER......############
|
308 |
|
|
######
|
309 |
|
|
|
310 |
|
|
#LC_names<-c("LC1_rec","LC2_rec","LC3_rec","LC4_rec","LC6_rec")
|
311 |
|
|
suf_name<-c("rec1")
|
312 |
|
|
sum_var<-c("diff_fc")
|
313 |
|
|
LC_names<-c("LC1","LC2","LC3","LC4","LC6")
|
314 |
|
|
brks<-c(-1,20,40,60,80,101)
|
315 |
|
|
lab_brks<-seq(1,5,1)
|
316 |
|
|
#var_name<-LC_names; suffix<-"rec1"; s_function<-"mean";df<-snot_data_selected;summary_var<-"diff_fc"
|
317 |
|
|
#reclassify_df(snot_data_selected,LC_names,var_name,brks,lab_brks,suffix,summary_var)
|
318 |
|
|
|
319 |
|
|
#Calculate mean per land cover percentage
|
320 |
|
|
data_agg<-reclassify_df(snot_data_selected,LC_names,brks,lab_brks,suf_name,sum_var)
|
321 |
|
|
data_lc<-data_agg[[1]]
|
322 |
|
|
snot_data_selected<-data_agg[[2]]
|
323 |
|
|
|
324 |
|
|
by_name<-"rec1"
|
325 |
|
|
df_lc_diff_fc<-merge_multiple_df(data_lc,by_name)
|
326 |
|
|
|
327 |
|
|
###### FIG10: PLOT LAND COVER
|
328 |
|
|
png(paste("fig10_diff_prediction_tmax_diff_res_f_land cover",date_selected,out_prefix,".png", sep=""))
|
329 |
|
|
par(mfrow=c(1,2))
|
330 |
|
|
zones_stat<-df_lc_diff_fc #first land cover
|
331 |
|
|
#names(zones_stat)<-c("lab_brks","LC")
|
332 |
|
|
y_range<-range(as.vector(t(zones_stat[,-1])),na.rm=T)
|
333 |
|
|
lab_brks_mid<-c(10,30,50,70,90)
|
334 |
|
|
plot(lab_brks_mid,zones_stat[,2],type="b",ylim=y_range,col="black", lwd=2,
|
335 |
|
|
ylab="difference between fusion and CAI",xlab="land cover percent classes")
|
336 |
|
|
lines(lab_brks_mid,zones_stat[,3],col="red",type="b")
|
337 |
|
|
lines(lab_brks_mid,zones_stat[,4],col="blue",type="b")
|
338 |
|
|
lines(lab_brks_mid,zones_stat[,5],col="darkgreen",type="b")
|
339 |
|
|
lines(lab_brks_mid,zones_stat[,6],col="purple",type="b")
|
340 |
|
|
legend("topleft",legend=c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban"),
|
341 |
|
|
cex=1.2, col=c("black","red","blue","darkgreen","purple"),
|
342 |
|
|
lty=1,lwd=1.8)
|
343 |
|
|
title(paste("Prediction tmax difference and land cover ",date_selected,sep=""))
|
344 |
|
|
|
345 |
|
|
###NOW USE RESIDUALS FOR FUSION
|
346 |
|
|
sum_var<-"res_f"
|
347 |
|
|
suf_name<-"rec2"
|
348 |
|
|
data_agg2<-reclassify_df(snot_data_selected,LC_names,brks,lab_brks,suf_name,sum_var)
|
349 |
|
|
data_resf_lc<-data_agg2[[1]]
|
350 |
|
|
#snot_data_selected<-data_agg[[2]]
|
351 |
|
|
|
352 |
|
|
by_name<-"rec2"
|
353 |
|
|
df_lc_resf<-merge_multiple_df(data_resf_lc,by_name)
|
354 |
|
|
|
355 |
|
|
zones_stat<-df_lc_resf #first land cover
|
356 |
|
|
#names(zones_stat)<-c("lab_brks","LC")
|
357 |
|
|
lab_brks_mid<-c(10,30,50,70,90)
|
358 |
|
|
plot(lab_brks_mid,zones_stat[,2],type="b",ylim=y_range,col="black",lwd=2,
|
359 |
|
|
ylab="tmax residuals fusion ",xlab="land cover percent classes")
|
360 |
|
|
lines(lab_brks_mid,zones_stat[,3],col="red",type="b")
|
361 |
|
|
lines(lab_brks_mid,zones_stat[,4],col="blue",type="b")
|
362 |
|
|
lines(lab_brks_mid,zones_stat[,5],col="darkgreen",type="b")
|
363 |
|
|
lines(lab_brks_mid,zones_stat[,6],col="purple",type="b")
|
364 |
|
|
legend("topleft",legend=c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban"),
|
365 |
|
|
cex=1.2, col=c("black","red","blue","darkgreen","purple"),
|
366 |
|
|
lty=1,lwd=1.2)
|
367 |
|
|
title(paste("Prediction tmax residuals and land cover ",date_selected,sep=""))
|
368 |
|
|
#savePlot(paste("fig10_diff_prediction_tmax_diff_res_f_land cover",date_selected,out_prefix,".png", sep=""), type="png")
|
369 |
|
|
dev.off()
|
370 |
|
|
#### FIGURE11: res_f and res_c per land cover
|
371 |
|
|
|
372 |
|
|
sum_var<-"res_c"
|
373 |
|
|
suf_name<-"rec3"
|
374 |
|
|
data_agg3<-reclassify_df(snot_data_selected,LC_names,brks,lab_brks,suf_name,sum_var)
|
375 |
|
|
data_resc_lc<-data_agg3[[1]]
|
376 |
|
|
snot_data_selected<-data_agg3[[2]]
|
377 |
|
|
|
378 |
|
|
by_name<-"rec3"
|
379 |
|
|
df_lc_resc<-merge_multiple_df(data_resc_lc,by_name)
|
380 |
|
|
|
381 |
|
|
zones_stat<-df_lc_resc #first land cover
|
382 |
|
|
#names(zones_stat)<-c("lab_brks","LC")
|
383 |
|
|
png(paste("fig11_prediction_tmax_res_f_res_c_land cover",date_selected,out_prefix,".png", sep=""))
|
384 |
|
|
par(mfrow=c(1,2))
|
385 |
|
|
y_range<-range(as.vector(t(zones_stat[,-1])),na.rm=T)
|
386 |
|
|
lab_brks_mid<-c(10,30,50,70,90)
|
387 |
|
|
plot(lab_brks_mid,zones_stat[,2],type="b",ylim=y_range,col="black",lwd=2,
|
388 |
|
|
ylab="tmax residuals CAI",xlab="land cover percent classes")
|
389 |
|
|
lines(lab_brks_mid,zones_stat[,3],col="red",type="b")
|
390 |
|
|
lines(lab_brks_mid,zones_stat[,4],col="blue",type="b")
|
391 |
|
|
lines(lab_brks_mid,zones_stat[,5],col="darkgreen",type="b")
|
392 |
|
|
lines(lab_brks_mid,zones_stat[,6],col="purple",type="b")
|
393 |
|
|
legend("topleft",legend=c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban"),
|
394 |
|
|
cex=1.2, col=c("black","red","blue","darkgreen","purple"),
|
395 |
|
|
lty=1,lwd=1.2)
|
396 |
|
|
title(paste("Prediction tmax residuals CAI and land cover ",date_selected,sep=""))
|
397 |
|
|
|
398 |
|
|
#fig11b
|
399 |
|
|
zones_stat<-df_lc_resf #first land cover
|
400 |
|
|
#names(zones_stat)<-c("lab_brks","LC")
|
401 |
|
|
y_range<-range(as.vector(t(zones_stat[,-1])),na.rm=T)
|
402 |
|
|
lab_brks_mid<-c(10,30,50,70,90)
|
403 |
|
|
plot(lab_brks_mid,zones_stat[,2],type="b",ylim=y_range,col="black",lwd=2,
|
404 |
|
|
ylab="tmax residuals fusion ",xlab="land cover percent classes")
|
405 |
|
|
lines(lab_brks_mid,zones_stat[,3],col="red",type="b")
|
406 |
|
|
lines(lab_brks_mid,zones_stat[,4],col="blue",type="b")
|
407 |
|
|
lines(lab_brks_mid,zones_stat[,5],col="darkgreen",type="b")
|
408 |
|
|
lines(lab_brks_mid,zones_stat[,6],col="purple",type="b")
|
409 |
|
|
legend("topleft",legend=c("LC1_forest", "LC2_shrub", "LC3_grass", "LC4_crop", "LC6_urban"),
|
410 |
|
|
cex=1.2, col=c("black","red","blue","darkgreen","purple"),
|
411 |
|
|
lty=1,lwd=1.2)
|
412 |
|
|
title(paste("Prediction tmax residuals and land cover ",date_selected,sep=""))
|
413 |
|
|
#savePlot(paste("fig10_diff_prediction_tmax_diff_res_f_land cover",date_selected,out_prefix,".png", sep=""), type="png")
|
414 |
|
|
#savePlot(paste("fig11_prediction_tmax_res_f_res_c_land cover",date_selected,out_prefix,".png", sep=""), type="png")
|
415 |
|
|
dev.off()
|
416 |
|
|
ac_data_obj<-list(snot_data_selected,data_vf,data_vc,ac_tab)
|
417 |
|
|
names(ac_data_obj)<-c("snot_data_selected","data_vf", "data_vc","ac_tab")
|
418 |
|
|
return(ac_data_obj)
|
419 |
|
|
}
|