In MartinRoth/knmitransformer: Transformation program

library(data.table)
library(knmitransformer)

The knmitransformer pacakge can be used with both standard stations, as the ones used in the published Climate Scenarions by KNMI, and user definied stations. This document provides examples for both situations and gives step by step examples for calculating transformed series of precipitation, temperature, radiation and Makkink evaporation.

Standard Stations

Standard Stations refer to the stations already used in the published climate scenarios with reference period 1981-2010. These reference data can be also found in the KNMI'14-klimaatscenario's webpage

Precipitation

To calculate transformed precipitation series using the reference data from the period 1981-2010, we have to specify the scenario, horizon, and subscenario (please use ?TransformPrecip for more information)

subscenario <- "centr"
input       <- KnmiRefFile("KNMI14____ref_rrcentr___19810101-20101231_v3.2.txt")
ofile       <- "tmp.txt" # output file - used only temporary
scenario    <- "GL"
horizon     <- 2030

futPrecip <- TransformPrecip(input=input, ofile=NA,
                             scenario=scenario, horizon=horizon,
                             subscenario=subscenario)

head(futPrecip[, 1 : 5])

The first five lines of columns V2 - V5 of the transformed input represent the station number, rd northing and rd easting, lat and lon of each station.

For Temperature, Radiation and Makkink evaporation similar steps are needed:

Temperature

var      <- "tg" #for daily average temperature; other options are daily maximum (tx), minimum (tn) temperature
input    <- KnmiRefFile("KNMI14____ref_tg___19810101-20101231_v3.2.txt")
regions  <-  MatchRegionsOnStationId(ReadInput(var, input)$header[1, -1])
scenario <- "GL"
horizon  <- 2030

futTemp <- TransformTemp(input=input, ofile=NA,
                         scenario=scenario, horizon=horizon,
                         var=var, regions=regions)

head(futTemp[, 1 : 5])

Radiation

input    <- KnmiRefFile("KNMI14____ref_rsds___19810101-20101231_v3.2.txt")
scenario <- "GL"
horizon  <- 2030

futRad <- TransformRadiation(input=input, ofile=NA,
                         scenario=scenario, horizon=horizon)

head(futRad[, 1 : 5])

Makkink evaporation

inputTemp <- KnmiRefFile("KNMI14____ref_tg___19810101-20101231_v3.2.txt")
inputRad  <- KnmiRefFile("KNMI14____ref_rsds___19810101-20101231_v3.2.txt")
scenario  <- "GL"
horizon   <- 2030
regions   <-  MatchRegionsOnStationId(ReadInput("tg", inputTemp)$header[1, -1])

futEvap <- TransformEvap(inputTemp=inputTemp, inputRad=inputRad,
                         scenario=scenario, horizon=horizon,
                         regions = regions)

head(futEvap[, 1 : 5])

User definied stations

User definied stations, such as the data set in this example for daily temperature data set, which was downloaded from the Climate Explorer will need some actions beforehand. In this example data sets from De Bilt station are used.

As you can see the data set has a slightly different format:

filename <- system.file("exampledata",
                         "tgtg260.dat",
                         package="knmitransformer")
tgData <- fread(filename, skip = 6)
head(tgData, n = 10)

For such data sets we will need to define an extra date column, provide latitude and longitude details, and select an appropriate base period of 30 years. A 30 year period is important for statistical stability; keep in mind that the transformation programme was written with the base period 1981 - 2010 in mind.

First, we append one single date column, which is necessary:

tgData[, date := as.integer(format(as.Date(paste(V1, V2, V3, sep = "-"),
                        format = "%Y-%m-%d"), format = "%Y%m%d"))]
head(tgData, n = 10)

Then we subset the correct period and ensure that the column names are date and values. Latitude and longitude are provided as well.

tgData <- tgData[V1 %in% 1981 : 2010, .(date, values = V4)]
head(tgData, n = 10)
lat <- 52.1
lon <- 5.18

Now we can create the input object and transform it in the next step (use ?CreateKnmiTFInput for more information)

input <- CreateKnmiTFInput(tgData, 1, lat, lon)
inputTrans <- TransformTemp(input, scenario = "GL", var = "tg", horizon = 2030,
                            regions = "MON") 
head(inputTrans, n = 15)

The first five lines of the second column of the transformed input represent the station number, rd northing and rd easting, lat and lon.

For precipitation we can perform a similar exercise ...

filename <- system.file("exampledata",
                        "hom1951550.dat",
                        package="knmitransformer")
precipData <- fread(filename, skip = 6)

setnames(precipData, c("date", "values"))

precipData <- precipData[date %/% 10000 %in% 1981 : 2010, ]
head(precipData, n = 10)

input <- CreateKnmiTFInput(precipData, 550, lat, lon)

inputTrans <- TransformPrecip(input, scenario = "GL", horizon = 2030)
head(inputTrans, n = 15)

And for radiation also...

filename <- system.file("exampledata",
                        "qqqq260.dat",
                        package="knmitransformer")
radiationData <- fread(filename, skip = 6)

head(radiationData, n = 10)

radiationData[, date := as.integer(format(as.Date(paste(V1, V2, V3, sep = "-"),
                        format = "%Y-%m-%d"), format = "%Y%m%d"))]
head(radiationData, n = 10)

radiationData <- radiationData[V1 %in% 1981 : 2010, .(date, values = V4)]

head(radiationData, n = 10)

The unit here is w/m2 instead of kJ/m2. Therefore, we have to adapt the data first

radiationData[, values := values * 86.4]
head(radiationData, n = 10)

input <- CreateKnmiTFInput(radiationData, 260, lat, lon)

inputTrans <- TransformRadiation(input, scenario = "GL", horizon = 2030) 
head(inputTrans, n = 15)

The final example is for evaporation:

inputTemp <- CreateKnmiTFInput(tgData, 260, lat, lon)
inputRad  <- CreateKnmiTFInput(radiationData, 260, lat, lon)

inputTrans <- TransformEvap(inputTemp = inputTemp, inputRad = inputRad,
                            scenario = "GL", horizon = 2030,
                            regions = "MON")
head(inputTrans, n = 15)