R/griddap_split.R

In rerddapUtils: Miscellaneous Utilities for 'rerddap'

#' Title Get ERDDAP gridded data restricted to a given season of the year
#'
#' \code{griddap_season} uses the R program 'rerddap::griddap()' to extract environmental data
#' from an 'ERDDAP' server in an (time, z, y ,x) bounding box where time is restricted to a
#' given season of the year  (see below).  Arguments are the same as in 'rerddap::griddap()'
#' except for the added 'season' parameter.  'read' and 'fmt' options are ignored.
#' @param datasetx Anything coercable to an object of class info. So the output of a
#' call to \code{\link[rerddap]{info}}, or a datasetid, which will internally be passed
#' through \code{\link[rerddap]{info}}
#' @param ... Dimension arguments. See examples. Can be any 1 or more of the
#' dimensions for the particular dataset - and the dimensions vary by dataset.
#' For each dimension, pass in a vector of length two, with min and max value
#' desired. at least 1 required.
#' @param fields (character) Fields to return, in a character vector.
#' @param stride (integer) How many values to get. 1 = get every value, 2 = get every other value, etc.
#'                Default: 1 (i.e., get every value)
#' @param request_split A numeric vector indicating the number of splits for each dimension, used
#'   to segment the request into manageable chunks. This is particularly useful for large datasets.
#' @param fmt (character) One of:
#' -  nc: save output to a netcdf file
#' -  memory:  save output in a dataframe in memory
#' -  duckdb:  save data in a duckdb database
#' -  parquet:  save data as a partitioned parquet dataset in a directory
#' @param url A URL for an ERDDAP server. Default:
#' https://upwell.pfeg.noaa.gov/erddap/ - See [rerddap::eurl()] for
#' more information
#' @param store ignored
#' @param read ignored
#' @param callopts Curl options passed on to \code{\link[crul]{verb-GET}}
#' @param aggregate_file A string specifying:
#'  - if is.null(aggregate_file) then a temporary file is created that does not
#'       overwrite user space
#'  - if "fmt = nc" the path to the NetCDF file to store the results
#'  - if "fmt = duckdb" the path to the duckdb file to store the results
#'  - if "fmt = memory" the value is ignored
#' @return Varies by format requested:
#'  - if "fmt = nc" the path to the NetCDF file where the results are stored
#'  - if "fmt = duckdb" the path to the duckdb file where the results are stored
#'  - if "fmt = memory" the usual rerddap::griddap dataframe
#' @export
#'
#' @examples
#' myURL <- "https://coastwatch.pfeg.noaa.gov/erddap/"
#' response <- try(httr::HEAD(myURL, httr::timeout(20)), silent = TRUE)
#' if (inherits(response, "try-error")) {
#'   message("The ERDDAP\u2122 server is not responding")
#' } else
#' {
#' data(wind_info)
#' request_split <- list(time = 2, altitude = 1, latitude = 1, longitude = 1)
#' res <- try(griddap_split(wind_info,
#'                      time = c('2015-12-31','2016-01-01'),
#'                      latitude = c(20, 21),
#'                      longitude = c(220, 221),
#'                      fields = 'mod_current',
#'                      request_split = request_split
#'                      ), silent = TRUE)
#' if (inherits(res, "try-error")) {
#'    message("Unable to retrieve data from the ERDDAP\u2122 server")
#'    }
#' }
griddap_split <- function(datasetx, ..., fields = 'all', stride = 1, request_split = NULL, fmt = "nc",
                    url = rerddap::eurl(), store = rerddap::disk(),
                    read = TRUE, callopts = list(), aggregate_file = NULL) {

  x <- datasetx
  if (is.null(request_split)) {
    #print('no split is given')
    #print('this must be a vector the same length of the number of dimensions')
    #print('each elementof the vector is the number of splits in that dimension')
    #stop('stopped on error')
    cli::cli_warn(c(
      "{.arg split} must be provided.",
      "i" = "It must be a vector the same length as the number of dimensions.",
      "i" = "Each element is the number of splits for that dimension."
    ))
    return(NULL)
  }
  if ( (fmt == 'nc') | (fmt == 'duckdb')) {
    if (!is.null(aggregate_file )) {
      file_check <- file.exists(aggregate_file)
      if (file_check) {
        #print('aggregate file already exists')
        #print('either rename the exising file or change aggregate_file')
        #stop('program is stopping')
        cli::cli_warn(c(
          "Aggregate file already exists: {.path {aggregate_file}}",
          "i" = "Rename the existing file or choose a different {.arg aggregate_file}."
        ))
        return(NULL)
      }
    }
  }
  if ( fmt == 'parquet') {
    if (!is.null(aggregate_file )) {
      dir_check <- dir.exists(aggregate_file)
      if (dir_check) {
        #print('aggregate file already exists')
        #print('either rename the exising file or change aggregate_file')
        #stop('program is stopping')
        cli::cli_warn(c(
          "Directpry already exists: {.path {aggregate_file}}",
          "i" = "Rename the existing Directory or choose a different {.arg aggregate_file}."
        ))
        return(NULL)
      }
    }
  }

  dimargs <- list(...)
  if (length(dimargs) == 0) stop("no dimension arguments passed, see ?griddap")
  if (inherits(x, "info")) {
    url <- x$base_url
    message("info() output passed to x; setting base url to: ", url)
  } else {
    x <- rerddap::as.info(x, url)
  }
  if (attr(x, "type") != "griddap")
    stop("datasetid '", attr(x, "datasetid"), "' not of type griddap")
  check_dims(dimargs, x)
  if (!is.null(dimargs$time)) {
    check_time_range(dimargs, x)
  }
  call_list <- extract_rerddap_call(x, dimargs, stride, fields, url)
  dimargs <- call_list$dimargs
  dim_args <- call_list$dim_args
  dimVals <- call_list$dimVals
  url_base <- call_list$url_base
  fields <- call_list$fields
  nc_file <- call_list$nc_file
  split_dim <- define_split(dimVals, request_split)
  extract <- split_griddap_request(x, url, stride,
                                    request_split, split_dim, fields, fmt,
                                    callopts, nc_file, aggregate_file)
  return(extract)
}



# Split and Process 'rerddap::griddap()' Requests
#
# This function handles splitting and processing requests to the `rerddap::griddap()` service,
# facilitating the handling of large datasets by dividing the requests into manageable parts.
# It supports output in different formats, including NetCDF ('nc') and DuckDB ('duckdb'), and
# performs the necessary operations to aggregate the results into a single file or database.
#
# @param info An object containing the dataset information, typically obtained from a call to
#  `rerddap::info()` or similar function that provides metadata about the dataset.
# @param url The base URL of the ERDDAP server from which to fetch the data.
# @param stride An integer indicating the stride (step size) for the query, allowing for
#   data thinning to reduce the volume of data retrieved.
# @param request_split A numeric vector indicating how to split the dataset query along each
#   dimension to manage large datasets more effectively.
# @param split_dim A named list where each name corresponds to a dimension in the dataset and
#   the values indicate how to split each dimension for the request.
# @param fields A character vector specifying which fields (variables) to retrieve from the
#   dataset. A value of 'all' indicates that all available fields should be included.
# @param fmt The format of the output, with supported values including 'nc' for NetCDF files
#   and 'duckdb' for DuckDB database format.
# @param aggregate_file The path to the output file or database where the aggregated results
#   from the split requests should be stored.
#
# @return Depending on the `fmt`, either the path to the aggregated NetCDF file or the DuckDB
#   database file containing the data retrieved and processed from the split requests.
#
split_griddap_request <- function(info, url, stride,
                                  request_split, split_dim, fields,  fmt,
                                  callopts, nc_file, aggregate_file) {
  # if fmt is duckdb,  open the connection
  con_db <- NULL
  if (fmt == 'duckdb') {
    if (is.null(aggregate_file)) {
      aggregate_file <- tempfile("extract", tmpdir = tempdir(), fileext = '.duckdb')
    }
    drv <- duckdb::duckdb()
    con_db <- duckdb::dbConnect(drv, aggregate_file)
  }
  # if fmt is netcdf,  create new netcdf file and copy attributes
  if (fmt == 'nc') {
    if (is.null(aggregate_file)) {
      aggregate_file <- tempfile("extract", tmpdir = tempdir(), fileext = '.nc')
    }
    return <- create_nc_file(info, fields, nc_file, aggregate_file )
    return <- copy_attributes(info, fields, nc_file, aggregate_file)
  }
  # if fmt is parquet,  create new partitioned parquet directory
  if (fmt == 'parquet') {
    if (is.null(aggregate_file)) {
      aggregate_file <- file.path(tempdir(), "partitioned_parquet_files")
      dir.create(aggregate_file, recursive = TRUE, showWarnings = FALSE)
    }
  }
  griddapOptsNames <- c('datasetx',  names(split_dim), 'fields', "stride", "callopts")
  split_names <- names(split_dim)
  final_result <- NULL
  griddapOpts <- list(info)
  final_result <- recursive_extract(1, request_split, split_dim, griddapOpts,
                                    stride, callopts, final_result,
                                    info, fields, fmt, con_db, aggregate_file)
  if (fmt == 'duckdb') {
    duckdb::dbDisconnect(con_db, shutdown = TRUE)
    duckdb::duckdb_shutdown(drv)
    final_result <- aggregate_file
  } else if ((fmt == 'nc') | (fmt == 'parquet')) {
    final_result <- aggregate_file
  }
  return(final_result)
}


# Recursively Split and Extract Data Requests
#
# This function is designed to recursively split and process data requests based on specified
# dimensions and constraints. It is particularly useful for managing large datasets or requests
# that need to be broken down into smaller chunks to comply with server limitations or to
# optimize data retrieval. The function iterates through each dimension, applies constraints,
# and accumulates the results.
#
# @param level The current level of recursion, indicating the dimension being processed.
# @param request_split A vector specifying the number of splits for each dimension of the data
#   request. The length of this vector should match the number of dimensions in the dataset.
# @param split_dim A named list where each name corresponds to a dimension and each value
#   specifies the splitting criteria for that dimension.
# @param griddapOpts A list of options for the `rerddap::griddap()` call, modified recursively
#   to reflect the current split criteria.
# @param stride The stride (step size) for the data retrieval, which can help reduce the volume
#   of data retrieved by skipping over points.
# @param callopts A list of additional options to pass to the HTTP request function, useful for
#   configuring details such as timeouts or proxy settings.
# @param final_result An accumulator for the results from each recursive call, which can be
#   aggregated or processed further.
# @param info Dataset information, typically obtained from a call to `rerddap::info()` or a
#   similar function, providing metadata about the dataset.
# @param fields Specifies which fields (variables) to retrieve from the dataset. Can be 'all'
#   for retrieving all fields or a vector of specific field names.
# @param fmt The format of the output, which can influence how results are stored or aggregated.
# @param con_db A connection object for database formats (e.g., DuckDB) when `fmt` is set to
#   a database format. This parameter is used to manage database connections during data extraction.
# @param aggregate_file The path to the file where aggregated results should be stored, applicable
#   when data is being compiled into a single file.
#
# @return The final aggregated result of the data extraction process, the specifics of which
#   depend on the `fmt` parameter and the nature of the `final_result` accumulation.
#
recursive_extract <- function(level, request_split, split_dim, griddapOpts,
                              stride, callopts, final_result,
                              info, fields, fmt, con_db, aggregate_file) {
  split_names <- names(split_dim)

  # Base case: If the level exceeds the length of request_split, perform the extraction
  if (level > length(request_split)) {
    final_result <- partial_extract(fields, fmt, griddapOpts,  stride, callopts,
                                    final_result, con_db, aggregate_file)
    return(final_result)
  }

  # Recursive step: Iterate over the current dimension
  temp_name <- split_names[level]
  for (i in seq(1, request_split[[level]])) {
    temp_value <- get_dim_constraint(request_split, split_dim, level, i)
    griddapOpts[[temp_name]] <- temp_value
    # Recurse to the next level with updated griddapOpts
    final_result <- recursive_extract(level + 1, request_split, split_dim, griddapOpts,
                                      stride, callopts, final_result,
                                      info, fields, fmt, con_db, aggregate_file)
  }

  return(final_result)
}


# Perform a Partial Data Extraction and Aggregation
#
# This function facilitates a portion of the data extraction process from a `rerddap::griddap()` call,
# handling specific fields and formatting options. It supports aggregating the extracted data into
# various formats such as in-memory objects, DuckDB, or NetCDF files, based on the specified format.
#
# @param fields The fields (variables) to be extracted from the dataset.
# @param fmt The format for aggregating the extracted data: 'memory', 'duckdb', or 'nc' (NetCDF).
# @param griddapOpts A list of options for the `rerddap::griddap()` call, including dataset ID,
#   dimensions, and any other parameters necessary for the data request.
# @param stride An integer specifying the stride (step size) for the data retrieval, useful for
#   thinning the data.
# @param callopts A list of additional options to customize the HTTP request made by `rerddap::griddap()`.
# @param final_result An accumulator for the results, which varies based on the specified format. For
#   'memory', this would be a data structure holding the aggregated data; for 'duckdb' or 'nc', it might
#   not be directly used.
# @param con_db A connection object for the DuckDB database, used when `fmt` is 'duckdb'.
# @param aggregate_file The file path for storing the aggregated results when `fmt` is 'nc'.
#
# @return Depending on the specified `fmt`, returns the aggregated data in the chosen format. For
#   'memory', it returns the aggregated data structure; for 'duckdb', it may return the connection
#   object or an indication of success; for 'nc', it may return the path to the NetCDF file or an
#   indication of success.
#
partial_extract <- function(fields, fmt,  griddapOpts, stride, callopts,
                            final_result, con_db, aggregate_file) {
  names(griddapOpts)[1] <- 'datasetx'
  griddapOpts$fields <- fields
  griddapOpts$stride <- stride
  griddapOpts$callopts <- callopts
  if (fmt == "parquet") {
    griddapOpts$fmt = 'parquet'
  }
  extract <- suppressMessages(do.call(rerddap::griddap, griddapOpts))
  if (fmt == 'memory') {
    final_result <- aggregate_memory(extract, final_result)
  } else if (fmt == 'duckdb') {
    final_result <- aggregate_duckdb(extract, con_db)
  } else if (fmt == 'nc') {
    final_result <- aggregate_netcdf(extract, aggregate_file)
  } else if (fmt == 'parquet') {
    final_result <- aggregate_parquet(extract, aggregate_file)
  }
return(final_result)
}



# Calculate Dimension Constraints for Data Retrieval Requests
#
# This function determines the constraints for a specific dimension of a dataset, based on the
# provided splitting criteria. It's designed to support segmented data retrieval operations by
# calculating the subset bounds for each segment of the request. This is particularly useful in
# the context of large dataset queries where requests may need to be split to avoid server
# limitations or to manage memory more efficiently.
#
# @param request_split A numeric vector indicating the number of splits for each dimension of
#   the data request. The length of this vector should match the number of dimensions in the dataset.
# @param split_dim A list where each element corresponds to a dimension in the dataset and contains
#   the full range of values for that dimension. The list's names should match the dimension names.
# @param var_index The index (integer) of the variable within `split_dim` for which the constraint
#   is being calculated.
# @param loop_index The index (integer) within the specified split of the dimension to calculate
#   the constraint for. This helps identify the specific segment of the split dimension.
#
# @return A numeric vector of length two, representing the lower and upper bounds of the constraint
#   for the specified dimension segment. If the dimension is not split (`request_split[var_index] == 1`),
#   it returns the full range for that dimension.
#
get_dim_constraint <- function(request_split, split_dim, var_index,  loop_index){
  split_names <- names(split_dim)
  temp_name <- split_names[var_index]
  if(request_split[var_index] == 1) {
    temp_value <- split_dim[[var_index]]
    temp_value <- unlist(temp_value)
    temp_value <- c(temp_value[1], temp_value[length(temp_value)])
  } else {
    temp_value <- split_dim[[var_index]]
    temp_value <- temp_value[loop_index]
    temp_value <- unlist(temp_value, use.names=FALSE)
    temp_value <- c(temp_value[1], temp_value[length(temp_value)])
  }
  return(temp_value)
}

# Define Splits for Dimension Values Based on Requested Segments
#
# This function calculates and organizes the splits for each dimension of a dataset based on the
# specified number of segments. It is useful for preparing segmented data retrieval requests,
# particularly when working with large datasets or when server limitations necessitate breaking down
# requests into smaller chunks.
#
# @param dimVals A list where each element corresponds to a dimension in the dataset and contains
#   the full range of values for that dimension. The list's names should match the dimension names.
# @param request_split A numeric vector indicating the number of segments into which each dimension
#   should be split. The length of this vector should match the number of dimensions in `dimVals`.
#
# @return A list of the same structure as `dimVals`, where each dimension's values are further
#   organized into the specified number of segments. If a dimension is not to be split (`request_split`
#   value of 1), its original range is preserved.
#
define_split <- function(dimVals, request_split) {
  indices <- request_split > 1
  dim_splits <- list()
  for (i in seq_along(dimVals)) {
    temp_name <- names(dimVals)[i]
    temp_dim <- dimVals[[i]]
    no_cut <- request_split[[i]]
    if (no_cut == 1) {
      dim_splits[[temp_name]] <- temp_dim
    } else {
      dim_splits[[temp_name]] <- split(temp_dim, cut(seq_along(temp_dim), no_cut, labels = FALSE))
    }
  }
  return(dim_splits)
}

# Aggregate Extracted Data in Memory
#
# This function aggregates newly extracted data with an existing data structure in memory.
# It is designed to sequentially combine data extracts, typically used when data is being
# retrieved in segments and needs to be compiled into a single dataset.
#
# @param extract A data frame or list representing the newly extracted data to be aggregated.
#   This parameter is expected to have a structure compatible with `final_result`, particularly
#   in terms of column names and types if it's a data frame.
# @param final_result The current aggregate of previously extracted data. If `NULL`,
#   `extract` becomes the initial dataset for aggregation. This parameter should be a data
#   structure (e.g., a data frame or list) that can hold the combined result of successive
#   extractions.
#
# @return Returns an updated version of `final_result` that includes the data from `extract`.
#   If `final_result` was initially `NULL`, returns `extract`.
aggregate_memory <- function(extract, final_result){
  if (is.null(final_result)) {
    final_result <- extract
  } else {
    final_result$data <- rbind(final_result$data, extract$data)
  }
  return(final_result)
}

#Aggregate Extracted Data into a DuckDB Database
#
#This function facilitates the aggregation of newly extracted data into a DuckDB database.
#It checks if a table named 'extract' exists in the database and, depending on its presence,
#either creates the table and inserts the data or appends the data to an existing table.
#
#@param extract A data frame or list representing the newly extracted data to be aggregated
#  into the DuckDB database. This parameter should contain the actual data in a format
#  compatible with DuckDB table structures.
#@param con_db A DuckDB connection object created via `DBI::dbConnect()`. This connection
#  is used to interact with the DuckDB database for writing the extracted data.
#
#@return This function does not return a value but performs an operation that inserts or
#  appends data into a DuckDB database.
#
aggregate_duckdb <- function(extract, con_db){
  # Corrected function body with appropriate variable names
  query <- "SELECT table_name FROM information_schema.tables WHERE table_name = 'extract';"
  result <- DBI::dbGetQuery(con_db, query) # Corrected the connection object variable name
  table_exists <- nrow(result) > 0

  if (!table_exists) {
    DBI::dbWriteTable(con_db, "extract", extract$data)
  } else {
    DBI::dbWriteTable(con_db, "extract1", extract$data, append = TRUE)
    DBI::dbExecute(con_db, "INSERT INTO extract SELECT * FROM extract1")
    # Assuming cleanup is needed for the temporary table 'extract1'
    DBI::dbExecute(con_db, "DROP TABLE extract1")
  }
}



# Aggregate Extracted Data into a NetCDF File
#
# This function aggregates extracted data into an existing NetCDF file. It is particularly useful
# for sequentially adding data from multiple extracts into a single NetCDF file, handling
# alignment of dimensions and variables between the extract and the aggregate file. It supports
# time dimensions and other variables, converting time values appropriately if necessary.
#
# @param extract A list representing the newly extracted data, including metadata about the
#   extraction (e.g., file names, dimensions, variables). The `summary` element of this list
#   should contain the filename, dimensions (`dim`), and variables (`var`) of the extracted data.
# @param aggregate_file A string specifying the path to the NetCDF file into which the extracted
#   data should be aggregated. This file should already exist and be structured to accommodate
#   the extracted data.
#
# @return This function does not return a value but modifies the specified NetCDF file by
#   adding or updating data from the extract.
#
# @examples
# # Assuming `extract` contains the necessary structure and `aggregate_file` is the path to
# # your NetCDF file:
# aggregate_netcdf(extract, "path/to/your/aggregate_file.nc")
#
# # Note: This example assumes that `extract` and `aggregate_file.nc` are correctly prepared
# # and compatible. Actual usage will require specific setup of the NetCDF file and extracts.
#

aggregate_netcdf <- function(extract, aggregate_file) {
  extract_nc <- ncdf4::nc_open(extract$summary$filename)
  on.exit(ncdf4::nc_close(extract_nc), add = TRUE)
  root <- ncdf4::nc_open(aggregate_file, write = TRUE)
  on.exit(ncdf4::nc_close(root), add = TRUE)

  # Compute start/count indices for each dimension
  dim_indices <- lapply(names(extract$summary$dim), function(dim_name) {
    extract_vals <- extract$data[[dim_name]]
    if (dim_name == "time") {
      extract_vals <- as.numeric(lubridate::as_datetime(extract_vals))
    }
    file_vals <- ncdf4::ncvar_get(root, dim_name)
    idx <- which(file_vals >= min(extract_vals) & file_vals <= max(extract_vals))
    list(start = idx[1], count = length(idx))
  })

  #start <- rev(sapply(dim_indices, `[[`, "start"))
  #count <- rev(sapply(dim_indices, `[[`, "count"))
  start <- rev(vapply(dim_indices, `[[`, FUN.VALUE = numeric(1), "start"))
  count <- rev(vapply(dim_indices, `[[`, FUN.VALUE = numeric(1), "count"))
  # Write each field into the aggregate file
  lapply(names(extract$summary$var), function(field) {
    ncdf4::ncvar_put(root, field, ncdf4::ncvar_get(extract_nc, field),
                     start = start, count = count)
  })

  invisible(NULL)
}

aggregate_parquet <- function(extract, aggregate_file){
  download_file <- attr(extract, "path")
  no_files <- length(list.files(aggregate_file))
  filename <- sprintf("extract%02d.parquet", no_files + 1)
  file.copy(
    from = download_file,
    to = file.path(aggregate_file, basename(filename))
  )
  invisible(NULL)
}