R/read_vpfiles.R

In adokter/bioRad: Biological Analysis and Visualization of Weather Radar Data

#' Read a vertical profile (`vp`) from file
#'
#' @param file A string containing the path to a vertical profile
#' generated by [calculate_vp].
#'
#' @keywords internal
#'
#' @return An object of class [summary.vp()].
read_vp <- function(file) {
  if (!is.vpfile(file)) {
    warning(paste(file, "is not a vertical profile"))
    return(NULL)
  }
  # check input argument
  groups <- rhdf5::h5ls(file, recursive = FALSE)$name
  if (!("dataset1" %in% groups)) {
    stop("HDF5 file does not contain a /dataset1 group")
  }
  # extract quantities
  groups <- rhdf5::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 <- rhdf5::h5readAttributes(file, "how")
  attribs.what <- rhdf5::h5readAttributes(file, "what")
  attribs.where <- rhdf5::h5readAttributes(file, "where")

  # construct wavelength attribute from frequency attribute if possible:
  if(is.null(attribs.how$wavelength) & !is.null(attribs.how$frequency)){
    speed_of_light = 299792458
    attribs.how$wavelength = 100*speed_of_light/attribs.how$frequency
  }

  # 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"
  )
  if(is.null(attribs.what$source)) attribs.what$source=""
  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 (`vp`) or a list of vertical profiles
#' (`vp`) from files
#'
#' @param files A character vector containing the file names of
#' vertical profiles in ODIM HDF5 format generated by [calculate_vp].
#'
#' @export
#'
#' @return A single `vp` object or a list of `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)
#'
#' # to load multiple files at once:
#' read_vpfiles(c(vpfile, vpfile))
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 <- rhdf5::h5readAttributes(file, paste(group, "/what", sep = ""))
  whatgroup$quantity
}

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