## simple validation for value.DS
# Sys.setenv(LANGUAGE='en', LANG="en_EN.UTF-8") # switching communicates to English
library(esd)
## Loading required package: ncdf
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## Następujące obiekty zostały zakryte z 'package:base':
##
## as.Date, as.Date.numeric
##
##
## Attaching package: 'esd'
##
## Następujący obiekt został zakryty z 'package:base':
##
## subset.matrix
library(plotrix) # required for plotting taylor diagrams only
library(hexbin)
variab <- "tmean"
forcing <- "T2M"
results <- readRDS(file=paste0("value.DS.ERA.",forcing,variab,".rda"))
ecad <- as.monthly(subset(readRDS(paste0("ecad.",variab,".rda")),it=1979:2008))
anomaly.ecad <-anomaly(ecad)
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in max(x, na.rm = TRUE): brak argumentów w max; zwracanie wartości
## -Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in max(x, na.rm = TRUE): brak argumentów w max; zwracanie wartości
## -Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in max(x, na.rm = TRUE): brak argumentów w max; zwracanie wartości
## -Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in max(x, na.rm = TRUE): brak argumentów w max; zwracanie wartości
## -Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in max(x, na.rm = TRUE): brak argumentów w max; zwracanie wartości
## -Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in max(x, na.rm = TRUE): brak argumentów w max; zwracanie wartości
## -Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
anomaly.results <-anomaly(results)
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in max(x, na.rm = TRUE): brak argumentów w max; zwracanie wartości
## -Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in max(x, na.rm = TRUE): brak argumentów w max; zwracanie wartości
## -Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in max(x, na.rm = TRUE): brak argumentów w max; zwracanie wartości
## -Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in max(x, na.rm = TRUE): brak argumentów w max; zwracanie wartości
## -Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in max(x, na.rm = TRUE): brak argumentów w max; zwracanie wartości
## -Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
## Warning in max(x, na.rm = TRUE): brak argumentów w max; zwracanie wartości
## -Inf
## Warning in min(x, na.rm = TRUE): brak argumentów w min; zwracanie wartości
## Inf
ecad.df <- as.data.frame(anomaly.ecad)
results.df <- as.data.frame(anomaly.results)
#corr <- round(sapply(1:ncol(ecad.df), function(i) cor(ecad.df[,i], results.df[,i],use = "pairwise.complete.obs",method = "spearman")),2)
#dev.off()
#hist(apply(ecad.df, 2, mean)-apply(results.df, 2, mean))
plot(ecad.df[,1],results.df[,1], xlim=c(-15,15), ylim=c(-15,15), col="#FF505050", main=paste0(variab,": ","1979-2008", "\nforcing: ",forcing), xlab="ECAD", ylab="DOWNSCALED")
abline(0,1); grid(lwd=1, col="black")
for (i in (1:ncol(ecad.df))) points(ecad.df[,i], results.df[,i], col="#0000FF20", cex=0.7, pch=19)
#dev.off()
hex <- hexbin(as.numeric(coredata(ecad)),as.numeric(coredata(results)))
plot(hex)
hist(as.numeric(unlist(ecad.df-results.df)), breaks=c(-25,-10:10,25), ylim=c(0,0.15))
abline(h=0.1)
# plotting taylor diagrams
# comparison for the whole period (1979-2008):
for (i in 1:dim(ecad.df)[2]){
mod <- results.df[,i]
obs <- ecad.df[,i]
correlation <- round(cor(mod,obs,use = "pairwise.complete.obs"),3)
bias <- round(mean(abs(mod-obs),na.rm=T),3) # mean absolute error
rmse <- round(sqrt(mean((mod - obs)^2, na.rm=T)),3)
station <- names(ecad.df)[i]
if(i==1) staty <- NULL # creating empty object for storing results
staty <- rbind(staty,cbind.data.frame(station,correlation,bias,rmse))
if(i==1) taylor.diagram(obs,mod,sd.arcs = T, col="#FF005050", ref.sd = T,normalize = T, main=paste0(variab,": ","1979-2008", "\nforcing: ",forcing))
if(i>1) taylor.diagram(obs,mod,sd.arcs = T, col="#FF005050", add=T, normalize=T)
# print(i) # checking progress
}
staty <- staty[with(staty, order(correlation)), ]
print(staty)
## station correlation bias rmse
## 38 BADAJOZ-TALAVERA LA REAL -0.008 2.053 2.559
## 19 LARISSA 0.010 2.329 2.908
## 70 SANTIAGO-DE-COMPOSTELA 0.022 2.079 2.584
## 80 MADRID-BARAJAS 0.022 2.221 2.752
## 18 CORFU 0.033 1.890 2.341
## 20 METHONI 0.036 1.533 1.943
## 34 BRAGANCA 0.039 2.421 2.952
## 42 TORTOSA-OBSERVATORIO-DEL-EBRO 0.045 1.897 2.327
## 35 LISBOA-GEOFISICA 0.058 1.703 2.094
## 40 NAVACERRADA 0.059 2.769 3.538
## 39 MALAGA 0.071 1.314 1.650
## 79 PALMA-DE-MALLORCA 0.078 1.851 2.249
## 41 SAN-SEBASTIAN-IGUELDO 0.081 2.349 2.928
## 65 TOULOUSE-BLAGNAC 0.105 2.598 3.219
## 24 CAGLIARI 0.107 1.778 2.209
## 23 MILAN 0.127 2.805 3.442
## 25 ROMA-CIAMPINO 0.127 2.099 2.560
## 13 MARSEILLE-MARIGNANE 0.132 2.335 2.918
## 59 SIBIU 0.144 2.905 3.762
## 74 HVAR 0.157 2.097 2.616
## 77 CONSTANTA 0.159 2.574 3.405
## 58 MONT-AIGOUAL 0.162 2.979 3.652
## 72 SION 0.176 2.498 3.124
## 55 STORNOWAY 0.183 1.462 1.832
## 26 VERONA-VILLAFRANCA 0.186 2.290 2.796
## 47 ESKDALEMUIR 0.193 1.860 2.313
## 56 VALLEY 0.197 1.671 2.083
## 37 BUCURESTI-BANEASA 0.198 2.686 3.491
## 50 RENNES 0.199 2.398 2.946
## 44 LUGANO 0.202 1.957 2.455
## 45 SAENTIS 0.212 3.164 3.856
## 73 GOSPIC 0.213 2.859 3.639
## 57 WADDINGTON 0.218 2.073 2.562
## 17 ZUGSPITZE 0.229 3.050 3.746
## 36 ARAD 0.233 2.882 3.710
## 11 BOURGES 0.234 2.647 3.273
## 75 ZAVIZAN 0.234 3.020 3.743
## 4 SONNBLICK 0.235 2.957 3.637
## 12 PARIS-14E 0.263 2.541 3.133
## 2 INNSBRUCK 0.273 2.679 3.320
## 46 ZUERISWITZERLAND 0.282 2.777 3.454
## 83 OBERSTDORF 0.285 2.798 3.443
## 43 BASEL-BINNINGEN 0.286 2.705 3.408
## 7 ZAGREB-GRIC 0.300 2.729 3.494
## 31 VARDOE 0.300 1.953 2.488
## 78 RHEINSTETTEN 0.305 2.745 3.505
## 15 HOHENPEISSENBERG 0.310 3.163 3.952
## 71 JACKVIK 0.320 3.816 5.378
## 66 IASI 0.340 3.006 3.894
## 3 SALZBURG 0.345 2.802 3.564
## 27 KARASJOK 0.356 4.261 5.578
## 1 GRAZ 0.359 2.425 3.123
## 84 REGENSBURG 0.363 2.588 3.300
## 69 TAFJORD 0.371 2.420 3.120
## 81 GIESSEN-WETTENBERG 0.372 2.542 3.255
## 30 UTSIRA 0.391 1.900 2.391
## 76 BROCKEN 0.393 2.840 3.473
## 51 FOKSTUA 0.408 2.752 3.504
## 28 KJOEREMSGRENDE 0.411 3.066 3.986
## 5 WIEN 0.418 2.527 3.257
## 10 SODANKYLA 0.435 3.964 5.289
## 14 BREMEN 0.436 2.559 3.275
## 62 HELGOLAND 0.446 2.038 2.603
## 54 HAPARANDA 0.456 3.373 4.601
## 63 DRESDEN-KLOTZSCHE 0.472 2.714 3.514
## 16 POTSDAM 0.486 2.580 3.358
## 85 SALEN 0.487 3.183 4.414
## 21 VESTERVIG 0.508 2.378 3.038
## 82 ARKONA 0.509 2.020 2.614
## 29 FAERDER 0.515 2.370 3.086
## 86 SIIKAJOKI-REVONLAHTI 0.528 3.310 4.537
## 60 GOTEBORG A 0.550 2.642 3.405
## 53 SIEDLCE 0.569 2.619 3.604
## 52 LEBA 0.595 2.105 2.838
## 9 JYVASKYLA-LENTOASEMA 0.619 3.014 4.135
## 33 KLAIPEDA 0.625 2.426 3.273
## 61 VISBY 0.625 2.027 2.647
## 22 TRANEBJERG 0.626 2.166 2.592
## 68 LAZDIJAI 0.639 2.597 3.554
## 32 KAUNAS 0.642 2.631 3.586
## 64 JOKIOINEN-JOKIOISTEN 0.642 2.798 3.829
## 8 HELSINKI-KAISANIEMI 0.651 2.579 3.484
## 67 BIRZAI 0.655 2.597 3.605
## 6 UCCLE NA NaN NaN
## 48 OXFORD NA NaN NaN
## 49 WICK NA NaN NaN
# simple boxplots:
boxplot(results.df, col = "#FF00FF50", main=paste0(variab,": ","1979-2008", "\nforcing: ",forcing))
boxplot(ecad.df,add=T, col="#00FF0050")
legend("topleft",legend = c("model","ecad"),col = c("#FF00FF20","#00FF0020"), lty=1,lwd=3, cex=0.8)
plot(subset(ecad, is=52), col="#0000FF95", lwd=2)
lines(subset(results, is=52), col="#FF0000", lwd=2)
dev.off()
## null device
## 1