R/read_netcdf.R

In andrew-MET/mepsr: Utilities for MEPS data

#' Read a 2d field from NetCDF file.
#'
#' \code{read_netcdf} returns a 2d field from a NetCDF file for a given
#' member and lead time. This function has only been tested on MEPS data.
#' The levels argument will currently only get data on pressure levels.
#'
#'
#' @param filename NetCDF file name to read from.
#' @param parameter Parameter to read - may be a parameter in the netcdf file
#'   or a standard HARP parameter name.
#' @param member The member to read from the file.
#' @param lead_time The lead time to read from the file.
#' @param level The pressure level for the data
#'
#' @return A 2-d array.
#' @export
#'
#' @examples
#' file_name <- "/lustre/storeB/immutable/archive/projects/metproduction/MEPS/2018/01/01/meps_extracted_2_5km_20180101T00Z.nc"
#' model_field <- read_netcdf(file_name, "air_temperature_2m", 0, 0)
#' model_field <- read_netcdf(file_name, "T2m", 0, 0)
#'
read_netcdf <- function(filename, parameter, member, lead_time, level = NULL, oper_sfx = FALSE) {
  #
  if (!requireNamespace("ncdf4", quietly = TRUE)) {
    stop("Package ncdf4 required for read_netcdf() - Please install from CRAN",
      call. = FALSE
    )
  }
  #
  if (!requireNamespace("stringr", quietly = TRUE)) {
    stop("Package stringr required for read_netcdf() - Please install from CRAN",
      call. = FALSE
    )
  }
  #
  ncID      <- ncdf4::nc_open(filename)
  nc_fields <- names(ncID$var)
  if (parameter %in% nc_fields) {
    nc_param <- parameter
  } else {
    nc_param  <- get_netcdf_param_MET(parameter)
  }
  if (is.na(nc_param)) stop("Requested parameter ", parameter, " is unknown")
  if (stringr::str_detect(nc_param, "_pl") && is.null(level)) {
    stop("Pressure level must be supplied for ", parameter, "with e.g. level = 1000")
  }
  #
  nc_lead_times	  <- ncdf4::ncvar_get(ncID, "time")
  nc_lead_times   <- (nc_lead_times - nc_lead_times[1]) / 3600
  nc_members      <- ncdf4::ncvar_get(ncID, "ensemble_member")

  lead_time_index <- which(nc_lead_times %in% lead_time)
  if (length(lead_time_index) == 0) stop("Requested lead time ", lead_time, " not found")
  member_index    <- which(nc_members %in% member)
  if (length(member_index) == 0) stop("Requested member ", member, " not found")
  if (!is.null(level)) {
    nc_pressure     <- try(ncdf4::ncvar_get(ncID, "pressure"), silent = TRUE)
    if (inherits(nc_pressure, "try-error")) {
      print("Trying to read var name: surface")
      nc_pressure <- ncdf4::ncvar_get(ncID, "surface")
      if (!inherits(nc_pressure, "try-error")) print("success!")
    }

    level_index   <- which(nc_pressure %in% level)
    if(length(level_index) == 0) stop("Requested pressure level ", level, "not found")
  } else {
    level_index   <- 1
  }
  #
  if (oper_sfx) {
    nc_start <- c(1, 1, member_index, lead_time_index)
    nc_count <- c(-1, -1, 1, 1)
  } else {
    nc_start <- c(1, 1, member_index, level_index, lead_time_index)
    nc_count <- c(-1, -1, 1, 1, 1)
  }
  nc_data <- ncdf4::ncvar_get(ncID, nc_param, start = nc_start, count = nc_count)
  ncdf4::nc_close(ncID)
  nc_data
}