library(ncdf4)
library(gtools)
library(verification)
## Loading required package: fields
## Loading required package: spam
## Loading required package: grid
## Spam version 1.0-1 (2014-09-09) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
##
## Następujące obiekty zostały zakryte z 'package:base':
##
## backsolve, forwardsolve
##
## Loading required package: maps
## Loading required package: boot
##
## Attaching package: 'boot'
##
## Następujące obiekty zostały zakryte z 'package:gtools':
##
## inv.logit, logit
##
## Loading required package: CircStats
## Loading required package: MASS
## Loading required package: dtw
## Loading required package: proxy
##
## Attaching package: 'proxy'
##
## Następujący obiekt został zakryty z 'package:spam':
##
## as.matrix
##
## Następujące obiekty zostały zakryte z 'package:stats':
##
## as.dist, dist
##
## Następujący obiekt został zakryty z 'package:base':
##
## as.matrix
##
## Loaded dtw v1.17-1. See ?dtw for help, citation("dtw") for use in publication.
library(plyr)
##
## Attaching package: 'plyr'
##
## Następujący obiekt został zakryty z 'package:fields':
##
## ozone
library(plotrix)
##
## Attaching package: 'plotrix'
##
## Następujący obiekt został zakryty z 'package:fields':
##
## color.scale
library(hydroGOF)
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## Następujące obiekty zostały zakryte z 'package:base':
##
## as.Date, as.Date.numeric
library(ggplot2)
library(magrittr)
library(dplyr)
##
## Attaching package: 'dplyr'
##
## Następujące obiekty zostały zakryte z 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
##
## Następujący obiekt został zakryty z 'package:MASS':
##
## select
##
## Następujące obiekty zostały zakryte z 'package:stats':
##
## filter, lag
##
## Następujące obiekty zostały zakryte z 'package:base':
##
## intersect, setdiff, setequal, union
# wczesniejsze ustalenia koordynat
geo_hopen <- c(25.04,76.30)
geo_hornsund <- c(15.30,77.00)
geo_bjornoya <- c(19.01,74.31)
#
# bjornoya <- NULL
# hopen <- NULL
# hornsund <- NULL
#
# katalogi <- dir("~/Pulpit/svalbard/156.17.59.39/wrf_svalbard/")
# domena <- "d03"
#
# for (i in 1:length(katalogi)){
# setwd(paste0("~/Pulpit/svalbard/156.17.59.39/wrf_svalbard/",katalogi[i]))
#
# plik=dir(pattern="d03_2008-06")
#
# for (pliki in 1:30){
# cat(paste(katalogi[i],domena),plik[pliki],"\n")
#
# nc<-nc_open(plik[pliki], write=F,readunlim=F,verbose=F) # otwieranie pliku
# lon =ncvar_get(nc,varid="XLONG") ; lon <- lon[,,1]
# lat =ncvar_get(nc,varid="XLAT") ; lat <- lat[,,1]
# #
# ind_hopen <- which( (abs(abs(lon-geo_hopen[1]) + (abs(lat-geo_hopen[2])))) == min(abs(abs(lon-geo_hopen[1]) + (abs(lat-geo_hopen[2])))),arr.ind=T) #109,63
# ind_hornsund <- which( (abs(abs(lon-geo_hornsund[1]) + (abs(lat-geo_hornsund[2])))) == min(abs(abs(lon-geo_hornsund[1]) + (abs(lat-geo_hornsund[2])))),arr.ind=T) #58,74
# ind_bjornoya <- which( (abs(abs(lon-geo_bjornoya[1]) + (abs(lat-geo_bjornoya[2])))) == min(abs(abs(lon-geo_bjornoya[1]) + (abs(lat-geo_bjornoya[2])))),arr.ind=T )# 80,15
#
#
# U10 =ncvar_get(nc,varid="U10")
# V10 =ncvar_get(nc,varid="V10")
# PSFC =ncvar_get(nc,varid="PSFC")
# T2 =ncvar_get(nc,varid="T2")
# WS <- sqrt(U10^2 +V10^2)
#
# time =ncvar_get(nc,varid="Times")
# time=strptime(time,format="%Y-%m-%d_%H:%M:%S", tz="GMT")
#
# bjornoya_tmp <- cbind.data.frame(time,WS[ind_bjornoya[1],ind_bjornoya[2],], T2[ind_bjornoya[1],ind_bjornoya[2],]-273.15,PSFC[ind_bjornoya[1],ind_bjornoya[2],]/100, katalogi[i], domena)
# bjornoya_tmp[,2:4] <- round(bjornoya_tmp[,2:4],1)
# names(bjornoya_tmp)[2:6] <- c("WS","T2","PSFC","katalog","domena")
#
# hornsund_tmp <- cbind.data.frame(time,WS[ind_hornsund[1],ind_hornsund[2],], T2[ind_hornsund[1],ind_hornsund[2],]-273.15,PSFC[ind_hornsund[1],ind_hornsund[2],]/100, katalogi[i], domena)
# hornsund_tmp[,2:4] <- round(hornsund_tmp[,2:4],1)
# names(hornsund_tmp)[2:6] <- c("WS","T2","PSFC","katalog","domena")
#
# hopen_tmp <- cbind.data.frame(time,WS[ind_hopen[1],ind_hopen[2],], T2[ind_hopen[1],ind_hopen[2],]-273.15,PSFC[ind_hopen[1],ind_hopen[2],]/100, katalogi[i], domena)
# hopen_tmp[,2:4] <- round(hopen_tmp[,2:4],1)
# names(hopen_tmp)[2:6] <- c("WS","T2","PSFC","katalog","domena")
#
# bjornoya <- rbind.data.frame(bjornoya,bjornoya_tmp)
# hopen <- rbind.data.frame(hopen,hopen_tmp)
# hornsund <- rbind.data.frame(hornsund,hornsund_tmp)
#
# nc_close(nc)
#
# }
#
# }
#
#rm(list=(c("bjornoya_tmp","hopen_tmp","hornsund_tmp","pliki","U10","V10","PSFC","lon","lat", "T2","WS","time","nc","i")))
#save(hopen,hornsund,bjornoya,file = '/home/bartosz/Dokumenty/publikacje własne/2015/svalbard_springer/wrf.rdata')
load(file = '/home/bartosz/Dokumenty/publikacje własne/2015/svalbard_springer/wrf.rdata')
wczyt_danych <- function(wyjscie, sciezka){
tmp <- read.csv(sciezka,sep='\t',na.strings = "----")
tmp$DATA <- as.character(tmp$DATA)
time=strptime(tmp$DATA,format="%m/%d/%Y", tz="GMT")
time <- time + tmp$HOUR*3600
tmp$DATA <- time
tmp$WS <- tmp$WS/3.6
assign(wyjscie, tmp, envir = .GlobalEnv)
}
sciezka <- "/home/bartosz/Dropbox/Dokumenty/publikacje własne/2015/svalbard_springer/obs/bjornoya.tsv"
wczyt_danych("bjornoya_obs",sciezka)
wczyt_danych("hornsund_obs","/home/bartosz/Dropbox/Dokumenty/publikacje własne/2015/svalbard_springer/obs/hornsund.tsv")
wczyt_danych("hopen_obs","/home/bartosz/Dropbox/Dokumenty/publikacje własne/2015/svalbard_springer/obs/hopen.tsv")
names(bjornoya_obs)[1] <- "time"
names(hornsund_obs)[1] <- "time"
names(hopen_obs)[1] <- "time"
BJORN <- join(bjornoya_obs,bjornoya, by="time",type="inner")
HORN <- join(hornsund_obs,hornsund, by="time",type="inner")
HOPEN <- join(hopen_obs,hopen, by="time",type="inner")
BJORN <- cbind.data.frame("BJORNOYA",BJORN)
names(BJORN) <- c("stacja","time","hh","T2obs","DIRobs","WSobs","GUSTobs","SLPobs","SLP0obs","WS", "T2","PSFC","katalog","domena")
BJORN <- (select(BJORN,-hh,-GUSTobs))
HORN <- cbind.data.frame("HORN",HORN)
names(HORN) <- c("stacja","time","hh","T2obs","DIRobs","WSobs","GUSTobs","SLPobs","SLP0obs","WS", "T2","PSFC","katalog","domena")
HORN <- (select(HORN,-hh,-GUSTobs))
HORN$SLP0obs <- HORN$SLP0obs -1.4
HOPEN <- cbind.data.frame("HOPEN",HOPEN)
names(HOPEN) <- c("stacja","time","hh","T2obs","DIRobs","WSobs","GUSTobs","SLPobs","SLP0obs","WS", "T2","PSFC","katalog","domena")
HOPEN <- (select(HOPEN,-hh,-GUSTobs))
DANE <- (rbind.data.frame(BJORN,HORN,HOPEN))
DANE <- (filter(DANE, !grepl('shave_sim', katalog))) # usuwamy "niewypaly"
DANE <- DANE %>% select(-DIRobs)
DANE %>% filter(stacja=="BJORNOYA" & domena=="d03" ) %>%
filter(katalog=="sim1b") %>%
select(-time,-katalog,-domena,-stacja,-SLPobs) -> tmp
# tu fragment kodu z obliczeniem statystyk do tabeli w artykule
# rbind(round(c(cor(tmp[,6],tmp[,3],use = "pairwise.complete.obs"), me(tmp[,6],tmp[,3]), mae(tmp[,6],tmp[,3]), d(tmp[,6],tmp[,3])),3)) # SLP
# rbind(c(round(cor(tmp[,5],tmp[,1]),2), round(me(tmp[,5],tmp[,1]),2), round(mae(tmp[,5],tmp[,1]),2), round(d(tmp[,5],tmp[,1]),2))) #T2
# rbind(c(round(cor(tmp[,4],tmp[,2]),2), round(me(tmp[,4],tmp[,2]),2), round(mae(tmp[,4],tmp[,2]),2), round(d(tmp[,4],tmp[,2]),2))) # WS
###
# przykladowe ploty do sprawdzenia poprawnosci obliczen
par(oma=c(0,0,0,0))
par(mar=c(0,0,0,0))
#plot(verify((subset(DANE,domena=="d03"&katalog=="sim1b")[,5]),(subset(DANE,domena=="d03"&katalog=="sim1b")[,9]), baseline = 20^20, frcst.type = "cont", obs.type = "cont"))
## STARY KOD DLA CELOW WERYFIKACJI:
# taylor.diagram((subset(HOPEN,domena=="d03"&katalog=="sim1b")[,5]),(subset(HOPEN,domena=="d03"&katalog=="sim1b")[,9]),col = "red",sd.arcs = T, main="HOPEN")
# taylor.diagram((subset(HOPEN,domena=="d02"&katalog=="sim1b")[,5]),(subset(HOPEN,domena=="d02"&katalog=="sim1b")[,9]),col="orange",add=T)
# taylor.diagram((subset(HOPEN,domena=="d01"&katalog=="sim1b")[,5]),(subset(HOPEN,domena=="d01"&katalog=="sim1b")[,9]),col="yellow",add=T)
#
library(ggplot2)
tmp2 <- subset(HOPEN,domena=="d03"&katalog=="sim1b")
ggplot(tmp2,aes(x=T2,y=T2obs)) + geom_point()
library(magrittr)
library(dplyr)
subset(HORN,domena=="d03"&katalog=="sim3b") %>% # select(T2,T2.1) %>% cor()
qplot(T2, T2obs, data = ., geom = c("smooth","point")) %>% print
## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.
