1
|
#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
|
}
|