Nothing
#' Read a polar volume (`pvol`) from file
#'
#' @param file A string containing the path to a polar volume file
#' @param sort A logical value, when `TRUE` sort scans ascending
#' by elevation.
#' @param param An atomic vector of character strings, containing the names
#' of scan parameters to read. To read all scan parameters use 'all'.
#' @param lat Latitude in decimal degrees of the radar position. If not
#' specified, value stored in file is used. If specified, value stored in
#' file is overwritten.
#' @param lon Longitude in decimal degrees of the radar position. If not
#' specified, value stored in file is used. If specified, value stored in
#' file is overwritten.
#' @param height Height of the center of the antenna in meters above sea
#' level. If not specified, value stored in file is used. If specified, value
#' stored in file is overwritten.
#' @param elev_min Minimum scan elevation to read in degrees.
#' @param elev_max Maximum scan elevation to read in degrees.
#' @param verbose A logical value, whether to print messages (`TRUE`)
#' to console.
#' @param mount (deprecated) A character string with the mount point (a directory path)
#' for the Docker container.
#' @param local_install (deprecated) String with path to local vol2bird installation,
#' to use local installation instead of Docker container
#'
#' @return An object of class [pvol][summary.pvol], which is a list
#' containing polar scans, i.e. objects of class `scan`
#'
#' @export
#'
#' @details
#' Scan parameters are named according to the OPERA data information
#' model (ODIM), see Table 16 in the
#' [ODIM specification](https://github.com/adokter/vol2bird/blob/master/doc/OPERA2014_O4_ODIM_H5-v2.2.pdf).
#' Commonly available parameters are:
#' * `DBZH`, `DBZ`: (Logged) reflectivity factor (dBZ)
#' * `TH`, `T`: (Logged) uncorrected reflectivity factor (dBZ)
#' * `VRADH`, `VRAD`: Radial velocity (m/s). Radial
#' velocities towards the radar are negative, while radial velocities away
#' from the radar are positive
#' * `RHOHV`: Correlation coefficient (unitless). Correlation
#' between vertically polarized and horizontally polarized reflectivity
#' factor
#' * `PHIDP`: Differential phase (degrees)
#' * `ZDR`: (Logged) differential reflectivity (dB)
#'
#' @examples
#' # locate example volume file:
#' pvolfile <- system.file("extdata", "volume.h5", package = "bioRad")
#'
#' # print the local path of the volume file:
#' pvolfile
#'
#' # load the file:
#' example_pvol <- read_pvolfile(pvolfile)
#'
#' # print summary info for the loaded polar volume:
#' example_pvol
#'
#' # print summary info for the scans in the polar volume:
#' example_pvol$scans
#'
#' # copy the first scan to a new object 'scan'
#' scan <- example_pvol$scans[[1]]
#'
#' # print summary info for the new object:
#' scan
read_pvolfile <- function(file, param = c(
"DBZH", "DBZ", "VRADH", "VRAD", "WRADH", "WRAD", "TH", "T", "RHOHV",
"ZDR", "PHIDP", "CELL", "BIOLOGY", "WEATHER", "BACKGROUND"
),
sort = TRUE, lat, lon, height, elev_min = 0,
elev_max = 90, verbose = TRUE,
mount = dirname(file), local_install) {
if (!missing(local_install)) warning("argument 'local_install' has been deprecated")
if (!file.exists(file)) {
stop(paste0("'", file, "' does not exist in current working directory ('", getwd(), "')."))
}
tryCatch(
read_pvolfile_body(
file, param, sort, lat, lon,
height, elev_min, elev_max,
verbose, mount, local_install
),
error = function(err) {
rhdf5::h5closeAll()
stop(err)
}
)
}
# this is the actual function read_pvolfile, without error handling that checks
# for open HDF5 files
read_pvolfile_body <- function(file, param = c(
"DBZH", "DBZ", "VRADH", "VRAD", "TH", "T", "RHOHV",
"ZDR", "PHIDP", "CELL", "BIOLOGY", "WEATHER", "BACKGROUND"
),
sort = TRUE, lat, lon, height, elev_min = 0,
elev_max = 90, verbose = TRUE,
mount = dirname(file), local_install) {
# input checks
assertthat::assert_that(assertthat::is.flag(sort),
msg = "'sort' should be logical")
if (!missing(lat)) {
lat_msg <- "'lat' should be a single numeric between -90 and 90 degrees"
assertthat::assert_that(assertthat::is.number(lat), msg = lat_msg)
assertthat::assert_that(lat > -90, msg = lat_msg)
assertthat::assert_that(lat < 90, msg = lat_msg)
}
if (!missing(lon)) {
lon_msg <- "'lon' should be a single numeric between -360 and 360 degrees"
assertthat::assert_that(assertthat::is.number(lon), msg = lon_msg)
assertthat::assert_that(lon > -360, msg = lon_msg)
assertthat::assert_that(lon < 360, msg = lon_msg)
}
if (!missing(height)) {
height_msg <- "'height' should be a single positive number of meters above sea level"
assertthat::assert_that(assertthat::is.number(height), msg = height_msg)
assertthat::assert_that(height > 0, msg = height_msg)
}
# check file type. If not ODIM HDF5, try to convert from RSL
cleanup <- FALSE
if (rhdf5::H5Fis_hdf5(file)) {
if (!is.pvolfile(file)) {
stop("Failed to read HDF5 file.")
}
} else {
pvol_tmp <- tempfile()
rlang::check_installed("vol2birdR", format_reason_vol2bird("to read `NEXRAD` files."))
config <- vol2birdR::vol2bird_config()
vol2birdR::rsl2odim(file = file, config = config, pvolfile_out = pvol_tmp, verbose = verbose)
if (!is.pvolfile(pvol_tmp)) {
file.remove(pvol_tmp)
stop("converted file contains errors")
}
cleanup <- TRUE
file <- pvol_tmp
}
# extract scan groups
scans <- rhdf5::h5ls(file, recursive = FALSE)$name
scans <- scans[grep("dataset", scans)]
# extract elevations, and make selection based on elevation
elevs <- sapply(
scans,
function(x) {
rhdf5::h5readAttributes(
file,
paste(x, "/where", sep = "")
)$elangle
}
)
scans <- scans[elevs >= elev_min & elevs <= elev_max]
# extract attributes
h5struct <- rhdf5::h5ls(file)
h5struct <- h5struct[h5struct$group == "/", ]$name
attribs.how <- attribs.what <- attribs.where <- NULL
if ("how" %in% h5struct) {
attribs.how <- rhdf5::h5readAttributes(file, "how")
}
if ("what" %in% h5struct) {
attribs.what <- rhdf5::h5readAttributes(file, "what")
}
if ("where" %in% h5struct) {
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
}
vol.lat <- c(attribs.where$lat) # need the c() to convert single element matrix to single element vector
vol.lon <- c(attribs.where$lon)
vol.height <- c(attribs.where$height)
if (is.null(vol.lat)) {
if (missing(lat)) {
if (cleanup) {
file.remove(file)
}
stop("latitude not found in file, provide 'lat' argument")
}
}
if (!missing(lat)) {
vol.lat <- lat
}
if (is.null(vol.lon)) {
if (missing(lon)) {
if (cleanup) {
file.remove(file)
}
stop("longitude not found in file, provide 'lon' argument")
}
}
if (!missing(lon)) {
vol.lon <- lon
}
if (is.null(vol.height)) {
if (missing(height)) {
if (cleanup) {
file.remove(file)
}
stop("antenna height not found in file, provide 'height' argument")
}
}
if (!missing(height)) {
vol.height <- height
}
geo <- list(lat = vol.lat, lon = vol.lon, height = vol.height)
# 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"
}
}
}
# write height, lat, lon attributes (update with potential user-defined values)
attribs.where$height <- vol.height
attribs.where$lat <- vol.lat
attribs.where$lon <- vol.lon
# assemble attributes in list
attributes <- list(
how = attribs.how, what = attribs.what,
where = attribs.where
)
# read scan groups
data <- lapply(
scans,
function(x) {
read_pvolfile_scan(file, x, param, radar, datetime, geo, attributes)
}
)
# filter out NULL output from read_pvolfile_scan
valid_scans <- which(!sapply(data, is.null))
assertthat::assert_that(length(valid_scans) > 0, msg = paste("none of the requested scan parameters found in file", file))
if(length(valid_scans) < length(scans)){
warning(paste("ignoring",length(scans)-length(valid_scans),"scan(s) in file",file,"because requested scan parameter(s) are missing."))
}
data <- data[valid_scans]
# order by elevation
if (sort) {
data <- data[order(sapply(data, get_elevation_angles))]
}
# prepare output
output <- list(
radar = radar, datetime = datetime, scans = data,
attributes = attributes, geo = geo
)
class(output) <- "pvol"
if (cleanup) {
file.remove(file)
}
output
}
read_pvolfile_scan <- function(file, scan, param, radar, datetime, geo, attributes) {
h5struct <- rhdf5::h5ls(file, all = TRUE)
groups <- h5struct[h5struct$group == paste("/", scan, sep = ""), ]$name
groups <- groups[grep("data", groups)]
dtypes <- h5struct[startsWith(h5struct$group, paste("/", scan, "/data", sep = "")), ]
dtypes <- dtypes[dtypes$name == "data", ]$dtype
h5struct <- h5struct[h5struct$group == paste("/", scan, sep = ""), ]$name
# select which scan parameters to read
if (length(param) == 1 && param == "all") {
allParam <- TRUE
} else {
allParam <- FALSE
}
if (!allParam) {
quantityNames <- sapply(
groups,
function(x) {
rhdf5::h5readAttributes(
file,
paste(scan, "/", x, "/what",
sep = ""
)
)$quantity
}
)
groups <- groups[quantityNames %in% param]
dtypes <- dtypes[quantityNames %in% param]
if (length(groups) == 0) {
return(NULL)
}
}
# read attributes
attribs.how <- attribs.what <- attribs.where <- NULL
if ("how" %in% h5struct) {
attribs.how <- rhdf5::h5readAttributes(file, paste(scan, "/how", sep = ""))
}
if ("what" %in% h5struct) {
attribs.what <- rhdf5::h5readAttributes(file, paste(scan, "/what", sep = ""))
}
if ("where" %in% h5struct) {
attribs.where <- rhdf5::h5readAttributes(file, paste(scan, "/where", sep = ""))
}
# add attributes to geo list
geo$elangle <- c(attribs.where$elangle)
geo$rscale <- c(attribs.where$rscale)
geo$ascale <- c(360 / attribs.where$nrays)
geo$astart <- attribs.how$astart
# odim stores ranges as Km in package ranges are until now in meters
geo$rstart <- attribs.where$rstart * 1000
# read scan parameters
quantities <- mapply(
function(x, y) {
read_pvolfile_quantity(
file,
paste(scan, "/", x, sep = ""),
radar, datetime, geo, y
)
},
x = groups,
y = dtypes,
SIMPLIFY = FALSE
)
quantityNames <- sapply(quantities, "[[", "quantityName")
quantities <- lapply(quantities, "[[", "quantity")
names(quantities) <- quantityNames
# if wavelength is attribute is missing at the scan level, copy it from the pvol level
if (is.null(attribs.how$wavelength)) attribs.how$wavelength <- attributes$how$wavelength
output <- list(
radar = radar, datetime = datetime, params = quantities,
attributes = list(
how = attribs.how, what = attribs.what,
where = attribs.where
), geo = geo
)
class(output) <- "scan"
output
}
read_pvolfile_quantity <- function(file, quantity, radar, datetime, geo, dtype) {
data <- rhdf5::h5read(file, quantity)$data
# convert storage mode from raw to numeric:
storage.mode(data) <- "numeric"
attr <- rhdf5::h5readAttributes(file, paste(quantity, "/what", sep = ""))
data <- replace(data, data == as.numeric(attr$nodata), NA)
data <- replace(data, data == as.numeric(attr$undetect), NaN)
data <- as.numeric(attr$offset) + as.numeric(attr$gain) * data
if (attr$quantity == "RHOHV") {
data <- replace(data, data > 10, NaN)
}
conversion <- list(
gain = as.numeric(attr$gain), offset = as.numeric(attr$offset),
nodata = as.numeric(attr$nodata), undetect = as.numeric(attr$undetect),
dtype = dtype
)
class(data) <- c("param", class(data))
attributes(data)$radar <- radar
attributes(data)$datetime <- datetime
attributes(data)$geo <- geo
attributes(data)$param <- as.character(attr$quantity)
attributes(data)$conversion <- conversion
list(quantityName = attr$quantity, quantity = data)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.