Revision 1e63098f
Added by Benoit Parmentier over 10 years ago
| climate/research/oregon/interpolation/contribution_of_covariates_paper_interpolation.R | ||
|---|---|---|
| 1085 | 1085 |
file_name<-paste("table6b_paper","_",out_prefix,".txt",sep="")
|
| 1086 | 1086 |
write.table(table6b,file=file_name,sep=",") |
| 1087 | 1087 |
|
| 1088 |
##Add new figure: |
|
| 1088 |
############################### |
|
| 1089 |
##Add new figure and analyses |
|
| 1089 | 1090 |
|
| 1090 |
s.table_term_tb <- l_s_table_term_tb[[2]] #baseline 2 (lat*lon) |
|
| 1091 |
#s.table_term_tb <- l_s_table_term_tb[[2]] #baseline 2 (lat*lon)
|
|
| 1091 | 1092 |
s.table_term_tb <- l_s_table_term_tb[[1]] #baseline 1 (lat*lon+elev_s) |
| 1092 |
tt <- aggregate(s.table_term_tb$p_val_rec2~s.table_term_tb$month,FUN=mean)
|
|
| 1093 |
tt <- aggregate(s.table_term_tb$p_val_rec2~s.table_term_tb$month,FUN=mean)
|
|
| 1094 |
plot(tt)
|
|
| 1095 |
test <- subset(s.table_term_tb,mod_name=="mod4" & term_name=="s(LST)")
|
|
| 1096 |
tt<-aggregate(test$p_val_rec3~test$month,FUN=mean)
|
|
| 1097 |
plot(tt,type="l",ylim=c(0.2,1))
|
|
| 1098 |
test2 <- subset(s.table_term_tb,mod_name=="mod4" & term_name=="s(elev_s)")
|
|
| 1099 |
tt2<-aggregate(test2$p_val_rec2~test2$month,FUN=mean)
|
|
| 1100 |
plot(tt2)
|
|
| 1101 |
lines(tt2)
|
|
| 1093 |
tb_sig_p_val_rec2 <- aggregate(s.table_term_tb$p_val_rec2~s.table_term_tb$month,FUN=mean)
|
|
| 1094 |
plot(tb_sig_p_val_rec2)
|
|
| 1095 |
#Now prepare
|
|
| 1096 |
s_table_LST_mod4 <- subset(s.table_term_tb,mod_name=="mod4" & term_name=="s(LST)")
|
|
| 1097 |
tb_mod4_LST_rec3 <- aggregate(s_table_LST_mod4$p_val_rec3~s_table_term_mod4$month,FUN=mean)
|
|
| 1098 |
plot(tb_mod4_LST_rec3,type="l",ylim=c(0.2,1))
|
|
| 1099 |
s_table_elev_mod4 <- subset(s.table_term_tb,mod_name=="mod4" & term_name=="s(elev_s)")
|
|
| 1100 |
tb_mod4_elev_rec3 <- aggregate(tb_mod4_elev_rec3$p_val_rec2 ~ tb_mod4_elev_rec3$month,FUN=mean)
|
|
| 1101 |
plot(tb_mod4_elev_rec3)
|
|
| 1102 |
lines(tb_mod4_elev_rec3)
|
|
| 1102 | 1103 |
test1 <- subset(s.table_term_tb,mod_name=="mod1" & term_name=="s(elev_s)") |
| 1103 | 1104 |
tt1<-aggregate(test1$p_val_rec2~test1$month,FUN=mean) |
| 1104 | 1105 |
plot(tt1) |
| 1105 | 1106 |
#model with LST and elev |
| 1106 | 1107 |
tb1_mod4_month <- raster_prediction_obj_1$summary_month_metrics_v[[4]] #note that this is for model1 |
| 1107 | 1108 |
#model with lev only |
| 1108 |
b1_mod1_month <- raster_prediction_obj_1$summary_month_metrics_v[[1]] #note that this is for model1 |
|
| 1109 |
tb1_mod1_month <- raster_prediction_obj_1$summary_month_metrics_v[[1]] #note that this is for model1
|
|
| 1109 | 1110 |
|
| 1110 | 1111 |
tb1_mod1_month<-raster_prediction_obj_1$summary_month_metrics_v[[1]] #note that this is for model1 |
| 1111 |
plot((tb1_mod4_month$rme) - tb1_mod1_month$rmse) |
|
| 1112 |
rmse_dif <- (tb1_mod4_month$rmse) - tb1_mod1_month$rmse |
|
| 1113 |
plot(rmse_dif) |
|
| 1112 | 1114 |
|
| 1113 | 1115 |
#now get correlation with LST and elev |
| 1114 | 1116 |
list_data_s <-extract_list_from_list_obj(raster_prediction_obj_1$validation_mod_obj,"data_s") |
| ... | ... | |
| 1116 | 1118 |
|
| 1117 | 1119 |
data_v_combined <-convert_spdf_to_df_from_list(list_data_v) #long rownames |
| 1118 | 1120 |
|
| 1119 |
|
|
| 1120 |
#LSTD_bias_avgm<-aggregate(LSTD_bias~month,data=data_month_all,mean) |
|
| 1121 |
#LSTD_bias_sdm<-aggregate(LSTD_bias~month,data=data_month_all,sd) |
|
| 1122 | 1121 |
data_v_combined |
| 1123 | 1122 |
LST_avgm <- aggregate(LST~month+id,data=data_v_combined,mean) |
| 1124 | 1123 |
#test <- aggregate(LST+dailyTmax~month+id,data=data_v_combined,mean) |
| ... | ... | |
| 1142 | 1141 |
plot(df_cor$LST_tmax,ylim=c(-1,1),col="red",type="l") |
| 1143 | 1142 |
lines(df_cor$elev_tmax,ylim=c(-1,1),col="blue") |
| 1144 | 1143 |
|
| 1144 |
#### Now plot figure for paper: |
|
| 1145 |
|
|
| 1146 |
|
|
| 1147 |
res_pix<-480 |
|
| 1148 |
col_mfrow<-3 |
|
| 1149 |
row_mfrow<-1 |
|
| 1150 |
png_file_name<- paste("Figure_17_paper_","LST_elev_",out_prefix,".png", sep="")
|
|
| 1151 |
png(filename=file.path(out_dir,png_file_name), |
|
| 1152 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 1153 |
#png(filename=file.path(out_dir,png_file_name), |
|
| 1154 |
#width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
|
| 1155 |
par(mfrow=c(row_mfrow,col_mfrow)) |
|
| 1156 |
|
|
| 1157 |
plot(tb_mod4_LST_rec3,type="l",ylab="Proportion of significant LST term", |
|
| 1158 |
xlab="Month") |
|
| 1159 |
title("Proportion of significant LST term in mod4",cex=1.5)
|
|
| 1160 |
|
|
| 1161 |
plot(rmse_dif,type="l",ylab="ΔRMSE between mod1 and mod4", |
|
| 1162 |
xlab="Month") |
|
| 1163 |
abline(h=0,lty="dashed") |
|
| 1164 |
title("Monthly ΔRMSE betwen mod1 and mod4",cex=1.5)
|
|
| 1165 |
|
|
| 1166 |
plot(df_cor$LST_tmax,ylim=c(-1,1),col="red",type="l", |
|
| 1167 |
ylab="Pearson Correlation",xlab="Month") |
|
| 1168 |
legend("topright",legend=c("Elev-tmax","LST-tmax"),col=c("red","blue"),cex=1,bty="n")
|
|
| 1169 |
|
|
| 1170 |
lines(df_cor$elev_tmax,ylim=c(-1,1),col="blue") |
|
| 1171 |
title("Monthly correlation for LST-tmax and elev-tmax",cex=1.5)
|
|
| 1172 |
|
|
| 1173 |
dev.off() |
|
| 1174 |
|
|
| 1175 |
names(tb_mod4_LST_rec3) <- c("month","p_val_rec3")
|
|
| 1176 |
tb_mod4_LST_rec3$month <- c("month","p_val_rec3")
|
|
| 1177 |
|
|
| 1178 |
res_pix<-480 |
|
| 1179 |
layout_m <- c(1,3) # works if set to this?? ok set the resolution... |
|
| 1180 |
#date_selected <- c("20100101","20100901")
|
|
| 1181 |
png(paste("Figure_17_paper_","LST_elev_",out_prefix,".png", sep=""),
|
|
| 1182 |
height=res_pix*layout_m[1],width=res_pix*layout_m[2]) |
|
| 1183 |
#height=480*6,width=480*4) |
|
| 1184 |
|
|
| 1185 |
p_prop <- xyplot(p_val_rec3 ~ as.numeric(month),data=tb_mod4_LST_rec3, |
|
| 1186 |
col=c("black"),
|
|
| 1187 |
type="b", |
|
| 1188 |
ylab=list(label="\u0394RMSE between mod1 and mod4",cex=1.5), |
|
| 1189 |
xlab=list(label="Month",cex=1.5), |
|
| 1190 |
main=list(label="Proportion of significant LST term in mod4",cex=1.8)) |
|
| 1191 |
|
|
| 1192 |
p_dif <- xyplot(rmse_dif ~ 1:12, |
|
| 1193 |
col=c("black"),
|
|
| 1194 |
type="b", |
|
| 1195 |
ylab=list(label="\u0394RMSE between mod1 and mod4",cex=1.5), |
|
| 1196 |
xlab=list(label="Month",cex=1.5), |
|
| 1197 |
main=list(label="Monthly \u0394RMSE betwen mod1 and mod4",cex=1.8)) |
|
| 1198 |
p_dif <- update(p_dif, panel = function(...) {
|
|
| 1199 |
panel.abline(h = 0, lty = 2, col = "gray") |
|
| 1200 |
panel.xyplot(...) |
|
| 1201 |
}) |
|
| 1202 |
|
|
| 1203 |
# Reshape the data: |
|
| 1204 |
df_cor$id <- 1:nrow(df_cor) #add a column before "melt", solution works well for any data.frame |
|
| 1205 |
m <- melt(df_cor, id = "id") |
|
| 1206 |
|
|
| 1207 |
p_cor<- xyplot(value ~ id, data = m, groups = variable, |
|
| 1208 |
#col=c("black","red"),
|
|
| 1209 |
type="b", |
|
| 1210 |
ylim=c(-1,1), |
|
| 1211 |
ylab=list(label="Pearson Correlation",cex=1.5), |
|
| 1212 |
xlab=list(label="Month",cex=1.5), |
|
| 1213 |
main=list(label="Pearson Correlation",cex=1.8), |
|
| 1214 |
auto.key = list("topright", corner = c(0,1),# col=c("black","red"),
|
|
| 1215 |
border = FALSE, lines = TRUE,cex=1.2)) |
|
| 1216 |
p_cor <- update(p_cor, panel = function(...) {
|
|
| 1217 |
panel.abline(h = 0, lty = 2, col = "gray") |
|
| 1218 |
panel.xyplot(...) |
|
| 1219 |
}) |
|
| 1220 |
#par.settings = list(axis.text = list(font = 2, cex = 1.3),layout=layout_m, |
|
| 1221 |
#par.main.text=list(font=2,cex=2),strip.background=list(col="white")),par.strip.text=list(font=2,cex=1.5), |
|
| 1222 |
|
|
| 1223 |
#grid.arrange(p1,p2,p3,ncol=1) |
|
| 1224 |
grid.arrange(p_prop,p_dif,p_cor,ncol=3) |
|
| 1225 |
|
|
| 1226 |
dev.off() |
|
| 1227 |
|
|
| 1228 |
|
|
| 1145 | 1229 |
######################## |
| 1146 | 1230 |
### Prepare table 7: correlation matrix between covariates |
| 1147 | 1231 |
|
| ... | ... | |
| 1420 | 1504 |
png(filename=file.path(out_dir,png_file_name), |
| 1421 | 1505 |
width=col_mfrow*res_pix,height=row_mfrow*res_pix) |
| 1422 | 1506 |
par(mfrow=c(row_mfrow,col_mfrow)) |
| 1507 |
|
|
| 1423 | 1508 |
boxplot(data_v_ag$res_mod1~elev_rcstat,ylim=c(-15,15),ylab="Residuals (deg C)",main="Residuals vs elevation classes for mod1=lat*lon", |
| 1424 | 1509 |
xlab="Elevation classes (meter)",names=c("0-500","500-1000","1000-1500","1500-2000"))
|
| 1425 | 1510 |
boxplot(data_v_ag$res_mod5~elev_rcstat,ylim=c(-15,15),ylab="Residuals (deg C)",main="Residuals vs elevation classes for mod5=lat*lon+LST", |
| 1426 | 1511 |
xlab="Elevation classes (meter)",names=c("0-500","500-1000","1000-1500","1500-2000"))
|
| 1427 | 1512 |
boxplot(data_v_ag$res_mod2~elev_rcstat,ylim=c(-15,15),ylab="Residuals (deg C)",main="Residuals vs elevation classes for mod2=lat*lon+elev", |
| 1428 | 1513 |
xlab="Elevation classes (meter)",names=c("0-500","500-1000","1000-1500","1500-2000"))
|
| 1514 |
dev.off() |
|
| 1515 |
|
|
| 1516 |
layout_m <- c(1,3) # works if set to this?? ok set the resolution... |
|
| 1517 |
png(paste("Figure_16_paper_","residuals_MAE_",out_prefix,".png", sep=""),
|
|
| 1518 |
height=480*layout_m[1],width=480*layout_m[2]) |
|
| 1519 |
#height=480*6,width=480*4) |
|
| 1429 | 1520 |
|
| 1521 |
p_bw1<-bwplot(data_v_ag$res_mod1~elev_rcstat,do.out=F,ylim=c(-15,15), |
|
| 1522 |
ylab=list(label="Residuals (deg C)",cex=1.5), |
|
| 1523 |
xlab=list(label="Elevation classes (meter)",cex=1.5), |
|
| 1524 |
main=list(label="Residuals vs elevation for mod1=lat*lon",cex=1.8), |
|
| 1525 |
scales = list(x = list(at = c(1, 2, 3, 4), |
|
| 1526 |
labels = c("0-500","500-1000","1000-1500","1500-2000"))))
|
|
| 1527 |
|
|
| 1528 |
p_bw2 <- bwplot(data_v_ag$res_mod5~elev_rcstat,do.out=F,ylim=c(-15,15), |
|
| 1529 |
ylab=list(label="Residuals (deg C)",cex=1.5), |
|
| 1530 |
xlab=list(label="Elevation classes (meter)",cex=1.5), |
|
| 1531 |
main=list(label="Residuals vs elevation for mod5=lat*lon+LST",cex=1.8), |
|
| 1532 |
scales = list(x = list(at = c(1, 2, 3, 4), |
|
| 1533 |
labels = c("0-500","500-1000","1000-1500","1500-2000"))))
|
|
| 1534 |
|
|
| 1535 |
p_bw3 <- bwplot(data_v_ag$res_mod5~elev_rcstat,do.out=F,ylim=c(-15,15), |
|
| 1536 |
ylab=list(label="Residuals (deg C)",cex=1.5), |
|
| 1537 |
xlab=list(label="Elevation classes (meter)",cex=1.5), |
|
| 1538 |
main=list(label="Residuals vs elevation for mod5=lat*lon+LST",cex=1.8), |
|
| 1539 |
scales = list(x = list(at = c(1, 2, 3, 4), |
|
| 1540 |
labels = c("0-500","500-1000","1000-1500","1500-2000"))))
|
|
| 1541 |
grid.arrange(p_bw1,p_bw2,p_bw3,ncol=3) |
|
| 1430 | 1542 |
dev.off() |
| 1431 | 1543 |
|
| 1432 | 1544 |
#layout_m <- c(1,3) # works if set to this?? ok set the resolution... |
Also available in: Unified diff
contribution of covariates and methods paper: revisions- breaking down significance and accuracy of LST by month