In MartinRoth/knmipluim: Return Levels for KNMI pluim

knitr::opts_chunk$set(echo = TRUE)

KNMI pluim

library(foreach)
library(iterators)
library(data.table)
library(ggplot2)
library(knmipluim)

The calculation for the current version of KNMI pluim climate information was done using the knmipluim package (version r packageVersion("knmipluim")).

windowSize <- 20L
inputPath <- system.file("OfficialData", package="knmipluim")
outputPath <- "./inst/OfficialOutput"
txFiles <- c("KNMI14____ref_tx___19810101-20101231_v3.2.txt",
             "KNMI14_GL_2050_tx___19810101-20101231_v3.2.txt",
             "KNMI14_WH_2050_tx___19810101-20101231_v3.2.txt")
tnFiles <- c("KNMI14____ref_tn___19810101-20101231_v3.2.txt",
             "KNMI14_GL_2050_tn___19810101-20101231_v3.2.txt",
             "KNMI14_WH_2050_tn___19810101-20101231_v3.2.txt")

result <- foreach(f = iter(txFiles)) %do% {
  tmp <- fread(paste0(inputPath, "/", f))
  setnames(tmp, c("date", unlist(tmp[1, -1, with=FALSE])))
  tmp <- tmp[date != 0]
  tmp[, date := as.Date(paste(date), format = "%Y%m%d")]
  baseName <- strsplit(f, ".txt")[[1]]
  tmp <- melt(tmp, id.vars = "date", variable.name = "station", value.name = "tx")
  tmp[, station := as.character(station)]

  tmpStation <- foreach(s = iter(unique(tmp[, station]))) %do% {
    maxDat <- tmp[station == s]
    retClim <- ReturnLevelClimatology(maxDat, "tx", "max",
                                windowSize = windowSize,
                                kLoc = 15, kScale = 8, kShape = 4)
    if (s == "260") {
      p <- PlotReturnLevels(retClim)
      print(p + ggtitle(paste(baseName, "at De Bilt")))
    }
    retClim <- retClim$returnLevels
    retClim[, baseName := baseName]
    retClim[, station := s]
    write.csv(retClim, file = paste0(outputPath, "/tx_retLevel", baseName,"_", s, ".csv"), row.names = FALSE)
  } 
  return(tmpStation)
}

result <- foreach(f = iter(tnFiles)) %do% {
  tmp <- fread(paste0(inputPath, "/", f))
  setnames(tmp, c("date", unlist(tmp[1, -1, with=FALSE])))
  tmp <- tmp[date != 0]
  tmp[, date := as.Date(paste(date), format = "%Y%m%d")]
  baseName <- strsplit(f, ".txt")[[1]]
  tmp <- melt(tmp, id.vars = "date", variable.name = "station", value.name = "tn")
  tmp[, station := as.character(station)]

  tmpStation <- foreach(s = iter(unique(tmp[, station]))) %do% {
    minDat <- tmp[station == s]
    retClim <- ReturnLevelClimatology(minDat, "tn", "min",
                                windowSize = windowSize,
                                kLoc = 15, kScale = 8, kShape = 4)
    if (s == "260") {
      p <- PlotReturnLevels(retClim)
      print(p + ggtitle(paste(baseName, "at De Bilt")))
    }
    retClim <- retClim$returnLevels
    retClim[, baseName := baseName]
    retClim[, station := s]
    write.csv(retClim, file = paste0(outputPath, "/tn_retLevel", baseName,"_", s, ".csv"), row.names = FALSE)
  } 
  return(tmpStation)
}