R/read_vpfiles.R

Defines functions read_odim_profile_data quantity_name read_vpfiles read_vp

Documented in read_vp read_vpfiles

#' Read a vertical profile (\code{vp}) from file
#'
#' @param file A string containing the path to a vertical profile
#' generated by \link{calculate_vp}.
#'
#' @keywords internal
#'
#' @return An object of class \code{\link{summary.vp}}.
read_vp <- function(file) {
  if (!is.vpfile(file)) {
    warning(paste(file, "is not a vertical profile"))
    return(NULL)
  }
  # check input argument
  groups <- h5ls(file, recursive = FALSE)$name
  if (!("dataset1" %in% groups)) {
    stop("HDF5 file does not contain a /dataset1 group")
  }
  # extract quantities
  groups <- h5ls(file)
  groups <- groups[which(groups$name == "data"), ]$group
  quantities <- sapply(
    groups,
    function(x) {
      quantity_name(file, x)
    }
  )
  profile <- as.data.frame(lapply(
    groups,
    function(x) {
      read_odim_profile_data(file, x)
    }
  ))
  # rename "height" column to "height"
  quantities <- gsub("HGHT", "height", quantities)
  names(profile) <- quantities

  # extract attributes
  attribs.how <- h5readAttributes(file, "how")
  attribs.what <- h5readAttributes(file, "what")
  attribs.where <- h5readAttributes(file, "where")
  # add vp_filename attribute if missing
  if (is.null(attribs.how$filename_vp)) {
    attribs.how$filename_vp <- file
  }

  # convert some useful metadata
  datetime <- as.POSIXct(paste(attribs.what$date, attribs.what$time),
    format = "%Y%m%d %H%M%S", tz = "UTC"
  )
  sources <- strsplit(attribs.what$source, ",")[[1]]
  radar <- gsub("NOD:", "", sources[which(grepl("NOD:", sources))])
  if (length(radar) == 0) {
    radar <- gsub("RAD:", "", sources[which(grepl("RAD:", sources))])
    if (length(radar) == 0) {
      radar <- gsub("WMO:", "", sources[which(grepl("WMO:", sources))])
      if (length(radar) == 0) {
        radar <- "unknown"
      }
    }
  }

  # prepare output
  output <- list(
    radar = radar, datetime = datetime, data = profile,
    attributes = list(
      how = attribs.how, what = attribs.what,
      where = attribs.where
    )
  )
  class(output) <- "vp"
  output
}

#' Read a vertical profile (\code{vp}) or a list of vertical profiles
#' (\code{vp}) from files
#'
#' @param files A character vector containing the file names of
#' vertical profiles in ODIM HDF5 format generated by \link{calculate_vp}.
#'
#' @export
#'
#' @return A single \code{vp} object or a list of \code{vp} objects.
#'
#' @examples
#' # locate example profile file:
#' vpfile <- system.file("extdata", "profile.h5", package = "bioRad")
#'
#' # print the local path of the profile file:
#' vpfile
#'
#' # load the file:
#' read_vpfiles(vpfile)
#'
#' # load multiple files at once:
#' \dontrun{
#' # read_vpfiles(c("my/path/profile1.h5", "my/path/profile2.h5", ...))
#' }
#'
read_vpfiles <- function(files) {
  if (length(files) == 1) {
    return(read_vp(files))
  } else {
    vps <- lapply(files, read_vp)
    # remove nulls
    vps <- vps[!sapply(vps, is.null)]
    do.call(c.vp, vps)
  }
}

quantity_name <- function(file, group) {
  whatgroup <- h5readAttributes(file, paste(group, "/what", sep = ""))
  whatgroup$quantity
}

read_odim_profile_data <- function(file, group) {
  whatgroup <- h5readAttributes(file, sprintf("%s/what", group))
  nodata <- whatgroup$nodata
  undetect <- whatgroup$undetect
  gain <- whatgroup$gain
  offset <- whatgroup$offset
  data <- h5read(file, sprintf("%s/data", group))[1, ]
  data <- replace(data, data == nodata, NA)
  data <- replace(data, data == undetect, NaN)
  offset + gain * data
}

Try the bioRad package in your browser

Any scripts or data that you put into this service are public.

bioRad documentation built on May 10, 2022, 1:07 a.m.