Included EOBS into knmiR

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: knmiR
 Type: Package
 Title: Import KNMI Data
-Version: 0.1.5.5
+Version: 0.1.5.6
 Authors@R: person("Martin", "Roth", email = "[email protected]",
  role = c("aut", "cre"))
 Description: Import KNMI data into R.
@@ -31,4 +31,5 @@ Roxygen: list(markdown = FALSE)
 SystemRequirements: netcdf library version 4.1 or later
 Suggests:
  testthat,
- rgeos
+ rgeos,
+ raster

NAMESPACE

@@ -14,12 +14,13 @@ export(Call)
 export(Citation)
 export(ClipQuakes)
 export(Description)
+export(EOBS)
 export(Earthquakes)
 export(HomogenPrecip)
 export(KIS)
-export(KIStemplate)
 export(KnmiData)
 export(License)
+export(LocalImportEOBS)
 export(MetaData)
 export(StudyArea)
 import(Rcpp)

NEWS.md

@@ -3,3 +3,4 @@
 - exporting ClipQuakes (fixed small bug related to object class)
 - added StudyArea
 - added examples to KIS
+- added EOBS

R/EOBS.R

@@ -39,3 +39,211 @@
 # file.remove(filename)
 # return(var)
 # }
+
+#' Imports EOBS data
+#' @param variable String from ('tg', 'tn', 'tx, ...)
+#' @param period Either numeric, timeBased or ISO-8601 style
+#' (see \code{\link[xts]{.subset.xts}})
+#' @param area Either SpatialPolygons or SpatialPolygonsDataFrame object (see
+#' package sp)
+#' @param grid String from ('0.25reg', '0.50reg', '0.25rot', '0.50rot')
+#' @param na.rm Boolean indicating if rows with NA can be deleted
+#' @param download Boolean indicating whether to download
+#' @export
+EOBS <- function(variable, period, area, grid, na.rm=TRUE,
+ download=TRUE) {
+ SanitizeInputEOBS(variable, period, area, grid)
+ url <- specifyURL(variable, grid)
+ message(paste("Loading opendapURL", url))
+ data <- GetEOBS(url, variable, area, period, na.rm)
+ return(data)
+}
+
+#' Import EOBS data from local file
+#' @note Should be merged with \code{\link{importEOBS}}
+#' @inheritParams importEOBS
+#' @param filename String containing the path to the ncdf file
+#' @export
+LocalImportEOBS <- function(variable, filename, period = NULL, area = NULL,
+ na.rm=TRUE) {
+ # Local sanitizing
+ data <- GetEOBS(filename, variable, area, period, na.rm)
+ return(data)
+}
+
+GetEOBS <- function(filename, variable, area, period, na.rm) {
+
+ result <- GetEobsBbox(filename, variable, sp::bbox(area), period)
+ result <- CreateDataTableMelt(variable, result)
+
+ if ( !is.null(area) & !is.matrix(area)) {
+ result <- removeOutsiders(result, area)
+ }
+ if ( na.rm ) result <- removeNAvalues(result)
+ result[, year := as.numeric(format(time, "%Y"))]
+ result[, month := as.numeric(format(time, "%m"))]
+ result[, day := as.numeric(format(time, "%d"))]
+ setcolorder(result, c("time", "year", "month", "day", "lat", "lon", variable, "pointID"))
+ return(result)
+}
+
+
+# Checks whether the input to importEOBS is valid or not
+# @param variable Variable name
+# @param period Period
+# @param area Area
+# @param grid Grid
+SanitizeInputEOBS <- function(variable, period, area, grid) {
+ if (variable %in% c('tg_stderr', 'tn_stderr', 'tx_stderr', 'pp_stderr',
+ 'rr_stderr')) {
+ stop("Standard error of variables not yet implemented.")
+ }
+ else if (!variable %in% c('tg', 'tn', "tx", "pp", 'rr', 'tg_stderr', 'tn_stderr',
+ 'tx_stderr', 'pp_stderr', 'rr_stderr')) {
+ stop(paste("Variable", variable, "not known."))
+ }
+ tryCatch(xts::.parseISO8601(period),
+ warning = function(e) {stop()},
+ error = function(e) {
+ stop("Period should be either Numeric, timeBased or ISO-8601 style.")
+ })
+ if (!class(area) %in% c("SpatialPolygons","SpatialPolygonsDataFrame")) {
+ stop("Area should be of class SpatialPolygons or SpatialPolygonsDataFrame.")
+ }
+ if (!grid %in% c("0.25reg", "0.50reg", "0.25rot", "0.50rot")) {
+ stop("Grid should be specified correctly.")
+ }
+}
+
+# Specifies the url based on the variableName and the grid
+# @param variableName Variable name
+# @param grid Grid
+specifyURL <- function(variableName, grid) {
+ url <- 'http:https://opendap.knmi.nl/knmi/thredds/dodsC/e-obs_'
+ if (grid=="0.50reg") {
+ url <- paste(url, '0.50regular/', sep="")
+ ending <- '_0.50deg_reg_v14.0.nc'
+ }
+ if (grid=="0.25reg") {
+ url <- paste(url, '0.25regular/', sep="")
+ ending <- '_0.25deg_reg_v14.0.nc'
+ }
+ url <- paste(url, variableName, ending, sep="")
+ return(url)
+}
+
+# Get the EOBS netcdf dimensions
+GetEobsDimensions <- function(ncdfConnection) {
+ values <- list()
+ values$lat <- ncdf4::ncvar_get(ncdfConnection, varid = 'latitude')
+ values$lon <- ncdf4::ncvar_get(ncdfConnection, varid = 'longitude')
+ values$time <- ncdf4::ncvar_get(ncdfConnection, varid = 'time')
+ return(values)
+}
+
+# Accesses the OPeNDAB server
+# @param filename String either url or local file
+# @param variableName String which variable to get
+# @param bbox Bounding box of spatial object
+# @param period Time period
+# @note This function is based on the script by Maarten Plieger
+# https://publicwiki.deltares.nl/display/OET/OPeNDAP+subsetting+with+R
+GetEobsBbox = function(filename, variableName, bbox, period){
+
+ # Open the dataset
+ dataset = ncdf4::nc_open(filename)
+
+ # Get lon and lat variables, which are the dimensions of depth.
+ values <- GetEobsDimensions(dataset)
+
+ # Determine the valid range of the dimensions based on the period and
+ # the bounding box
+ validRange <- list()
+ validRange$time <- which(findInterval(values$time,
+ periodBoundaries(values$time, period))==1)
+ validRange$lat <- which(findInterval(values$lat, bbox[2,])==1)
+ validRange$lon <- which(findInterval(values$lon, bbox[1,])==1)
+
+ # Make a selection of indices which fall in our subsetting window
+ # E.g. translate degrees to indices of arrays.
+ determineCount <- function(x) {return(c(x[1], tail(x,1) - x[1] + 1))}
+ count <- rbind(determineCount(validRange$lon),
+ determineCount(validRange$lat),
+ determineCount(validRange$time))
+
+
+ # Prepare a list with the valued values of the dimensions and the variable
+ validValues <- list()
+ validValues$lat <- values$lat[validRange$lat]
+ validValues$lon <- values$lon[validRange$lon]
+ validValues$time <- as.Date(values$time[validRange$time],
+ origin="1950-01-01")
+ validValues[[variableName]] <- ncdf4::ncvar_get(dataset, variableName,
+ start=count[, 1],
+ count=count[, 2])
+
+ # Close the data set and return data.table created from the valid values
+ ncdf4::nc_close(dataset)
+ return(validValues)
+}
+
+CreateDataTableMelt <- function(variable, validValues) {
+ if (length(validValues$time) > 1) {
+ meltedValues <- reshape2::melt(validValues[[variable]],
+ varnames=c("lon", "lat", "time"))
+ result <- as.data.table(meltedValues)
+ } else {
+ meltedValues <- reshape2::melt(validValues[[variable]],
+ varnames=c("lon", "lat"))
+ result <- as.data.table(meltedValues)
+ result[, time := 1]
+ }
+ setkey(result, lon, lat)
+ result[, pointID:=.GRP, by = key(result)]
+ setkey(result, pointID)
+ index <- result[, !all(is.na(value)), by = pointID][V1==TRUE, pointID]
+ result <- result[pointID %in% index, ]
+ result[, pointID := NULL]
+ result[, lon := validValues$lon[lon]]
+ result[, lat := validValues$lat[lat]]
+ result[, time := validValues$time[time]]
+ setnames(result, "value", paste(variable))
+ return(result)
+}
+
+# Removes points outside of the SpatialPolygons
+# Not for external use
+# @param data Data.table
+# @param area Valid area
+removeOutsiders <- function(data, area) {
+ setkey(data, lon, lat)
+ data[, pointID:=.GRP, by=key(data)]
+ coords <- data[, .(lon = unique(lon),
+ lat = unique(lat)), by = pointID][, .(lon, lat)]
+ points <- sp::SpatialPoints(coords, area@proj4string)
+ index <- data[, unique(pointID)][which(!is.na(sp::over(points,
+ as(area, 'SpatialPolygons'))))]
+ data <- data[pointID %in% index]
+ setkey(data, lon, lat)
+ return(data[, pointID:=.GRP, by=key(data)])
+}
+
+# Removes all rows with NAs
+# Not for external use
+# @param data data.table
+removeNAvalues <- function(data) {
+ # We don't check if time is NA (it should not) but date * 0 is not defined
+ data <- data[complete.cases(data[,!"time", with=FALSE]*0)]
+ setkey(data, lon, lat)
+ data[, pointID:=.GRP, by=key(data)]
+}
+
+# To define the valid range
+# @param time Time
+# @param period Period
+periodBoundaries <- function(time, period) {
+ xts <- xts::xts(time, as.Date(time, origin="1950-01-01"))
+ interval <- range(as.numeric(xts[period]))
+ interval[2] <- interval[2] + 1
+ return(interval)
+}

R/KIS.R

@@ -44,17 +44,6 @@ SetRightColumnClass <- function(dt) {
  dt[, c(3) := as.numeric(eval(as.name(colName)))]
 }
 
-#' Python test template
-#' @param var variable
-#' @param geoIdentifier location
-#' @param period period
-#' @export
-KIStemplate <- function(var, geoIdentifier, period) {
- data.table(date = as.Date(0 : 19, origin = as.Date("2015-01-01")),
- loc = geoIdentifier,
- var = rnorm(20))
-}
-
 WriteKISRecipe <- function(var, locationID, period) {
  # period is not yet used in the recipe
  # max results does not seem to have any effect

R/knmiR.R

@@ -0,0 +1,17 @@
+#' knmiR: A package to import and KNMI data
+#'
+#' The package provides an R interface to some data sources of KNMI.
+#'
+#' Available data sources:
+#'
+#' Earthquake catalog via the function Earthquakes()
+#'
+#' KIS data base via the function KIS() (internally only)
+#'
+#'
+# #' @section eobsR functions:
+# #' Do the following
+#'
+#' @docType package
+#' @name knmiR
+NULL

tests/testthat/test-EOBS.R

@@ -0,0 +1,19 @@
+context("EOBS")
+
+library(knmiR)
+adm0 <- raster::getData('GADM', country='NL', level=0)
+context("Input testing:")
+
+expect_error(EOBS("foo"), "Variable foo not known.")
+expect_error(EOBS('tg', 'A', "foo"), "Period should be either Numeric, timeBased or ISO-8601 style.")
+expect_error(EOBS('tg', '2014', "foo"), "Area should be of class SpatialPolygons or SpatialPolygonsDataFrame.")
+expect_error(EOBS('tg', '2014', adm0, "foo"), "Grid should be specified correctly.")
+
+context("Output testing:")
+expect_equal_to_reference(EOBS('tg', '2014', adm0, '0.50reg'), file="EOBSreference/output.rds")
+expect_equal_to_reference(EOBS('rr', '2014', adm0, '0.50reg'), file="EOBSreference/output_rr.rds")
+expect_equal_to_reference(EOBS("tg", "2015-06-01", adm0, grid = "0.50reg"),
+ file = "EOBSreference/output_one_timestep.rds")
+expect_equal_to_reference(LocalImportEOBS("tg", "tg_0.50deg_reg_v12.0_plus_2015_ANN_avg.nc",
+ "2000/2015", adm0),
+ file = "EOBSreference/output_local.rds")