## 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
variab <- "tmean"
forcing <- "SLP"
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()
library(hexbin)
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))
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.025 1.941 2.421
## 27 KARASJOK -0.017 4.071 5.599
## 23 MILAN -0.016 2.511 3.121
## 34 BRAGANCA -0.012 2.286 2.768
## 80 MADRID-BARAJAS 0.000 2.053 2.552
## 55 STORNOWAY 0.001 1.347 1.699
## 35 LISBOA-GEOFISICA 0.006 1.606 2.011
## 69 TAFJORD 0.006 2.537 3.262
## 79 PALMA-DE-MALLORCA 0.006 1.606 2.020
## 10 SODANKYLA 0.015 3.979 5.449
## 45 SAENTIS 0.025 2.990 3.745
## 39 MALAGA 0.033 1.246 1.533
## 40 NAVACERRADA 0.034 2.650 3.293
## 31 VARDOE 0.037 1.780 2.315
## 17 ZUGSPITZE 0.039 2.853 3.611
## 70 SANTIAGO-DE-COMPOSTELA 0.040 1.896 2.353
## 25 ROMA-CIAMPINO 0.048 1.732 2.182
## 54 HAPARANDA 0.050 3.395 4.780
## 51 FOKSTUA 0.054 2.841 3.584
## 44 LUGANO 0.058 1.821 2.242
## 24 CAGLIARI 0.060 1.489 1.896
## 58 MONT-AIGOUAL 0.061 2.737 3.398
## 71 JACKVIK 0.062 3.627 4.974
## 28 KJOEREMSGRENDE 0.067 3.100 4.067
## 41 SAN-SEBASTIAN-IGUELDO 0.067 2.134 2.682
## 86 SIIKAJOKI-REVONLAHTI 0.073 3.646 4.960
## 4 SONNBLICK 0.075 2.664 3.395
## 13 MARSEILLE-MARIGNANE 0.081 1.974 2.465
## 26 VERONA-VILLAFRANCA 0.081 1.963 2.472
## 30 UTSIRA 0.083 1.905 2.428
## 65 TOULOUSE-BLAGNAC 0.099 2.276 2.852
## 85 SALEN 0.100 3.403 4.756
## 47 ESKDALEMUIR 0.103 1.607 2.022
## 18 CORFU 0.104 1.478 1.865
## 72 SION 0.108 2.183 2.721
## 20 METHONI 0.131 1.246 1.582
## 76 BROCKEN 0.131 2.815 3.471
## 42 TORTOSA-OBSERVATORIO-DEL-EBRO 0.141 1.569 1.963
## 74 HVAR 0.142 1.740 2.161
## 9 JYVASKYLA-LENTOASEMA 0.146 3.596 4.739
## 56 VALLEY 0.147 1.476 1.879
## 19 LARISSA 0.163 1.836 2.327
## 2 INNSBRUCK 0.171 2.377 2.955
## 57 WADDINGTON 0.178 1.790 2.296
## 11 BOURGES 0.189 2.327 2.898
## 75 ZAVIZAN 0.193 2.500 3.211
## 64 JOKIOINEN-JOKIOISTEN 0.203 3.312 4.371
## 29 FAERDER 0.208 2.455 3.262
## 15 HOHENPEISSENBERG 0.210 2.867 3.574
## 46 ZUERISWITZERLAND 0.218 2.411 3.024
## 83 OBERSTDORF 0.220 2.421 3.021
## 12 PARIS-14E 0.222 2.204 2.783
## 8 HELSINKI-KAISANIEMI 0.223 3.015 3.956
## 43 BASEL-BINNINGEN 0.231 2.387 2.990
## 73 GOSPIC 0.235 2.323 2.972
## 78 RHEINSTETTEN 0.242 2.502 3.107
## 21 VESTERVIG 0.243 2.368 3.112
## 60 GOTEBORG A 0.248 2.727 3.570
## 50 RENNES 0.260 1.975 2.467
## 61 VISBY 0.264 2.260 2.885
## 84 REGENSBURG 0.272 2.352 2.925
## 59 SIBIU 0.275 2.232 2.869
## 62 HELGOLAND 0.277 1.960 2.503
## 36 ARAD 0.282 2.315 2.945
## 81 GIESSEN-WETTENBERG 0.285 2.357 2.965
## 77 CONSTANTA 0.289 2.106 2.730
## 1 GRAZ 0.290 2.063 2.677
## 7 ZAGREB-GRIC 0.290 2.207 2.869
## 3 SALZBURG 0.298 2.460 3.095
## 14 BREMEN 0.299 2.432 3.082
## 82 ARKONA 0.308 1.997 2.564
## 37 BUCURESTI-BANEASA 0.320 2.135 2.726
## 22 TRANEBJERG 0.331 2.131 2.709
## 16 POTSDAM 0.345 2.470 3.162
## 67 BIRZAI 0.354 2.931 3.879
## 52 LEBA 0.367 2.190 2.886
## 5 WIEN 0.368 2.149 2.798
## 33 KLAIPEDA 0.368 2.600 3.408
## 63 DRESDEN-KLOTZSCHE 0.369 2.553 3.230
## 32 KAUNAS 0.391 2.784 3.695
## 66 IASI 0.399 2.379 3.166
## 68 LAZDIJAI 0.404 2.736 3.649
## 53 SIEDLCE 0.430 2.535 3.434
## 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