diff --git a/R/gefs.R b/R/gefs.R index 73e265af..18c7bf6d 100644 --- a/R/gefs.R +++ b/R/gefs.R @@ -58,24 +58,24 @@ #' #' #Get forecast for a certain variable. #' forecast <- gefs("Total_precipitation_surface_6_Hour_Accumulation_ens", -#' lat, lon) +#' lat, lon, ens = 0, time = 12) #' #' #Fetch a different date (available up to 10 days prior to today) #' forecast_yesterday_prec <- gefs( #' "Total_precipitation_surface_6_Hour_Accumulation_ens", -#' lat, lon, date=format(as.Date(Sys.time()) - 1, "%Y%m%d")) +#' lat, lon, ens = 1, time = 6, date=format(as.Date(Sys.time()) - 1, "%Y%m%d")) #' #' #specific ensemble and times, for the 1800 forecast. #' # here ensembles 1-3 (ensembles are numbered starting with 0) -#' # and time for 2 days from today at 1800 +#' # and two time periods: c(1800, 2400) #' date <- format(as.Date(Sys.time()) - 1, "%Y%m%d") #' var <- "Temperature_height_above_ground_ens" -#' gefs(var, lat, lon, date = date, forecast_time = "1800", ens_idx=2:4, -#' time_idx=1:8) +#' gefs(var, lat, lon, date = date, forecast_time = "1800", ens=1:3, +#' time=6*(3:4)) #' #' #One ensemble, all latitudes and longitudes (this is a big file) for the #' # next 3 days. -#' # gefs(var, ens=1, time=1:12) +#' # gefs(var, ens=1, time=6*(1:12)) #' } #' gefs <- function(var, lat, lon, ...) { @@ -107,11 +107,11 @@ gefs_CONNECT <- function(date = format(Sys.time(), "%Y%m%d"), } #' @rdname gefs -gefs_GET <- function(var, lat, lon, +gefs_GET <- function(var, lat = NULL, lon = NULL, ens = NULL, time = NULL, date = format(Sys.time(), "%Y%m%d"), forecast_time = c("0000", "0600", "1200", "1800"), - ens_idx = 1:21, - time_idx = 1:65, + ens_idx = NULL, # will be removed in future version + time_idx = NULL, # will be removed in future version dims = NULL, raw = FALSE, ...) { @@ -147,24 +147,42 @@ gefs_GET <- function(var, lat, lon, n_time <- varsize[ndims] #time is always the last dimension # Set the indices for each dimension + dim_names <- sapply(v$dim, function(d) { d$name }) dim_idxs <- list() for (i in 1:length(v$dim)) { dn <- v$dim[[i]]$name if(dn == "lon") { - dim_idxs[[i]] <- if(!missing(lon)) which(v$dim[[i]]$vals %in% (round(lon,0) %% 360)) else 1:v$dim[[i]]$len + dim_idxs[[i]] <- if(!is.null(lon)) which(v$dim[[i]]$vals %in% (round(lon,0) %% 360)) else 1:v$dim[[i]]$len } else if (dn == "lat") { - dim_idxs[[i]] <- if(!missing(lat)) which(v$dim[[2]]$vals %in% round(lat, 0)) else 1:v$dim[[2]]$len + dim_idxs[[i]] <- if(!is.null(lat)) which(v$dim[[i]]$vals %in% round(lat, 0)) else 1:v$dim[[i]]$len } else if (dn == "ens") { - dim_idxs[[i]] <- ens_idx - } else if (dn %in% c("time1", "time2")) { - dim_idxs[[i]] <- time_idx + dim_idxs[[i]] <- if(!is.null(ens)) which(v$dim[[i]]$vals %in% ens) else 1:v$dim[[i]]$len + } else if (dn %in% c("time", "time1", "time2")) { + dim_idxs[[i]] <- if(!is.null(ens)) which(v$dim[[i]]$vals %in% time) else 1:v$dim[[i]]$len } else if (dn %in% names(additional_dims)) { dim_idxs[[i]] <- which(v$dim[[j]]$vals %in% additional_dims[[dn]]) } else { dim_idxs[[i]] <- 1:v$dim[[i]]$len } } - names(dim_idxs) <- lapply(1:length(v$dim), function(i) { v$dim[[i]]$name }) + names(dim_idxs) <- dim_names + + # Assign ens_idx and time_idx if used + if(!is.null(ens_idx)) { + d_idx <- which(dim_names == "ens") + if(!is.null(ens_idx)) message("'ens_idx' is depreciated and will be removed in future versions, please specify values (not indices) in 'ens' instead.") + if(!any(ens_idx %in% 1:v$dim[[d_idx]]$len)) stop("'ens_idx' is out of bounds, check the dimension values with 'gefs_dimension_values(dim = 'ens').") + dim_idxs[['ens']] <- ens_idx + } + if(!is.null(time_idx)) { + time_name <- grep("time", dim_names, value = TRUE) + d_idx <- which(dim_names %in% time_name) + if(!is.null(time_idx)) message("'time_idx' is depreciated and will be removed in future versions, please specify values (not indices) in 'time' instead.") + if(!any(time_idx %in% 1:v$dim[[d_idx]]$len)) stop("'time_idx' is out of bounds, check the dimension values with 'gefs_dimension_values(dim = 'time').") + dim_idxs[[time_name]] <- time_idx + } + + #start indices of dimensions to read from data start <- sapply(dim_idxs, function(d) { min(d) }) @@ -189,7 +207,7 @@ gefs_GET <- function(var, lat, lon, forecast_time <- strsplit(fname, ".grib2")[[8]] list(forecast_date = date, forecast_time = forecast_time, - dimensions = dims, + dimensions = names(dim_idxs), data = d) } @@ -228,9 +246,31 @@ gefs_dimensions <- function(con = NULL, ...) { #' #' @param dim (character) the dimension. #' @rdname gefs -gefs_dimension_values <- function(dim, con = NULL, ...) { - if (is.null(dim) || missing(dim)) stop("dim cannot be NULL or missing.") +gefs_dimension_values <- function(dim, var = NULL, con = NULL, ...) { + if (missing(dim)) stop("dim cannot be NULL or missing.") if (is.null(con)) con = gefs_CONNECT(...) - con$dim[[dim]]$vals + + if (!(dim %in% names(con$dim))) { + stop(paste0(dim, " is not a valid GEFS dimension. Check with 'gefs_dimensions()'.")) + } + + if (!is.null(var)) { + v <- con$var[[var]] + dim_names <- sapply(v$dim, function(d) { d$name }) + + # there are multiple "time" dimensions, so get the one for this variable + if(dim == "time") dim <- grep("time", dim_names, value = TRUE) + + if(!(dim %in% dim_names)) stop(paste0(dim, + " is not in variable dimensions: ", + paste0(dim_names, collapse = ", "), + ".")) + + dim_idx <- which(dim_names == dim) + res = con$var[[var]]$dim[[dim_idx]]$vals + } else { + res <- con$dim[[dim]]$vals + } + res } diff --git a/man/gefs.Rd b/man/gefs.Rd index e093e5c2..fef1ed9b 100644 --- a/man/gefs.Rd +++ b/man/gefs.Rd @@ -16,9 +16,10 @@ gefs(var, lat, lon, ...) gefs_CONNECT(date = format(Sys.time(), "\%Y\%m\%d"), forecast_time = c("0000", "0600", "1200", "1800")) -gefs_GET(var, lat, lon, date = format(Sys.time(), "\%Y\%m\%d"), - forecast_time = c("0000", "0600", "1200", "1800"), ens_idx = 1:21, - time_idx = 1:65, dims = NULL, raw = FALSE, ...) +gefs_GET(var, lat = NULL, lon = NULL, ens = NULL, time = NULL, + date = format(Sys.time(), "\%Y\%m\%d"), forecast_time = c("0000", + "0600", "1200", "1800"), ens_idx = NULL, time_idx = NULL, + dims = NULL, raw = FALSE, ...) gefs_latitudes(con = NULL, ...) @@ -28,7 +29,7 @@ gefs_variables(con = NULL, ...) gefs_dimensions(con = NULL, ...) -gefs_dimension_values(dim, con = NULL, ...) +gefs_dimension_values(dim, var = NULL, con = NULL, ...) } \arguments{ \item{var}{the variable to get. Must be one of the variables listed in @@ -97,24 +98,24 @@ lon <- -118.2188 #Get forecast for a certain variable. forecast <- gefs("Total_precipitation_surface_6_Hour_Accumulation_ens", - lat, lon) + lat, lon, ens = 0, time = 12) #Fetch a different date (available up to 10 days prior to today) forecast_yesterday_prec <- gefs( "Total_precipitation_surface_6_Hour_Accumulation_ens", - lat, lon, date=format(as.Date(Sys.time()) - 1, "\%Y\%m\%d")) + lat, lon, ens = 1, time = 6, date=format(as.Date(Sys.time()) - 1, "\%Y\%m\%d")) #specific ensemble and times, for the 1800 forecast. # here ensembles 1-3 (ensembles are numbered starting with 0) -# and time for 2 days from today at 1800 +# and two time periods: c(1800, 2400) date <- format(as.Date(Sys.time()) - 1, "\%Y\%m\%d") var <- "Temperature_height_above_ground_ens" -gefs(var, lat, lon, date = date, forecast_time = "1800", ens_idx=2:4, - time_idx=1:8) +gefs(var, lat, lon, date = date, forecast_time = "1800", ens=1:3, + time=6*(3:4)) #One ensemble, all latitudes and longitudes (this is a big file) for the # next 3 days. -# gefs(var, ens=1, time=1:12) +# gefs(var, ens=1, time=6*(1:12)) } } diff --git a/tests/testthat/test-gefs.R b/tests/testthat/test-gefs.R index 9f248064..8a3614dc 100644 --- a/tests/testthat/test-gefs.R +++ b/tests/testthat/test-gefs.R @@ -13,12 +13,14 @@ test_that("gefs errors", { skip_on_travis() expect_error(gefs(lat=lat, lon=lon), "Need to specify the variable to get. A list of variables is available from gefs_variables().") - expect_error(gefs(var = temp, lat = c(-43, -41), lon = lons, ens_idx = 1, time_idx = 1), "Latitudes must be sequential.", fixed = TRUE) - expect_error(gefs(var = temp, lat = lats, lon = c(213, 211), ens_idx = 1, time_idx = 1), "Longitudes must be sequential.", fixed = TRUE) - expect_error(gefs(var = temp, lat = -91, lon = lons, ens_idx = 1, time_idx = 1), "Latitudes must be in c(-90,90).", fixed = TRUE) - expect_error(gefs(var = temp, lat = 91, lon = lons, ens_idx = 1, time_idx = 1), "Latitudes must be in c(-90,90).", fixed = TRUE) - expect_error(gefs(var = temp, lat = lats, lon = 361, ens_idx = 1, time_idx = 1), "Longitudes must be in c(-180,180) or c(0,360).", fixed = TRUE) - expect_error(gefs(var = temp, lat = lats, lon = -181, ens_idx = 1, time_idx = 1), "Longitudes must be in c(-180,180) or c(0,360).", fixed = TRUE) + expect_error(gefs(var = temp, lat = c(-43, -41), lon = lons, ens = 1, time = 6), "Latitudes must be sequential.", fixed = TRUE) + expect_error(gefs(var = temp, lat = lats, lon = c(213, 211), ens = 1, time = 6), "Longitudes must be sequential.", fixed = TRUE) + expect_error(gefs(var = temp, lat = -91, lon = lons, ens = 1, time = 6), "Latitudes must be in c(-90,90).", fixed = TRUE) + expect_error(gefs(var = temp, lat = 91, lon = lons, ens = 1, time = 6), "Latitudes must be in c(-90,90).", fixed = TRUE) + expect_error(gefs(var = temp, lat = lats, lon = 361, ens = 1, time = 6), "Longitudes must be in c(-180,180) or c(0,360).", fixed = TRUE) + expect_error(gefs(var = temp, lat = lats, lon = -181, ens = 1, time = 6), "Longitudes must be in c(-180,180) or c(0,360).", fixed = TRUE) + expect_error(gefs(var = temp, lat = lats, lon = lons, ens_idx = 22, time_idx = 1), "'ens_idx' is out of bounds, check the dimension values with 'gefs_dimension_values(dim = 'ens').", fixed = TRUE) + expect_error(gefs(var = temp, lat = lats, lon = lons, ens_idx = 1, time_idx = 67), "'time_idx' is out of bounds, check the dimension values with 'gefs_dimension_values(dim = 'time').", fixed = TRUE) }) test_that("gefs HTTP requests", { @@ -28,13 +30,13 @@ test_that("gefs HTTP requests", { ### Get raw and processed data d_raw <- gefs(var = temp, - ens_idx = 1:2, - time_idx = 1:2, + ens = 1:2, + time = c(6,12), forecast_time = "0000", lon = lons, lat = lats, raw = TRUE) d <- gefs(var = temp, - ens_idx = 1:2, - time_idx = 1:2, + ens = 0:1, + time = c(6,12), forecast_time = "0000", lon = lons, lat = lats) @@ -50,7 +52,7 @@ test_that("gefs HTTP requests", { expect_true(all(sort(unique(d$data$lon)) == sort(round(lons %% 360, 0)))) expect_true(all(sort(unique(d$data$lat)) == sort(round(lats, 0)))) - ### Tests of data transformation from multidimensional array to data frame + ### Tests of data transformation from multidimensional array to data frame grid <- expand.grid(lon = round(lons %% 360, 0), lat = round(lats, 0), ens = 1:2, time = 1:2) @@ -70,7 +72,29 @@ test_that("gefs HTTP requests", { d$data$time2 == 1,][[temp]])) }) - + +test_that("ens_idx and time_idx replace ens and time values", { + skip_on_cran() + skip_on_travis() + skip_on_appveyor() + + d <- gefs(var = temp, + ens = 0, + time = 6, + ens_idx = 1:2, + forecast_time = "0000", + lon = lons, lat = lats) + expect_true(all(0:1 %in% unique(d$data$ens))) + + d <- gefs(var = temp, + ens = 0, + time = 6, + time_idx = 3:4, + forecast_time = "0000", + lon = lons, lat = lats) + expect_true(all(c(12,18) %in% unique(d$data$time1))) +}) + test_that("gefs_variables returns characters.", { skip_on_cran() skip_on_travis() @@ -112,6 +136,12 @@ test_that("gefs_dimensions returns character list.", { expect_is(dims[1], "character") }) +test_that("gefs_dimension_values errors", { + expect_error(gefs_dimension_values(dim = "time2", var = temp), "time2 is not in variable dimensions: lon, lat, height_above_ground, ens, time1.", fixed = TRUE) + expect_error(gefs_dimension_values(), "dim cannot be NULL or missing.", fixed = TRUE) + expect_error(gefs_dimension_values(dim = "ens1"), "ens1 is not a valid GEFS dimension. Check with 'gefs_dimensions()'.", fixed = TRUE) +}) + test_that("gefs_dimension_values returns numeric array.", { skip_on_cran() skip_on_travis() @@ -120,6 +150,4 @@ test_that("gefs_dimension_values returns numeric array.", { vals = gefs_dimension_values("lat") expect_is(vals, "array") expect_is(vals[1], "numeric") - - expect_error(gefs_dimension_values(dim = NULL), "dim cannot be NULL or missing.") })