exactextractr coverage - 91.30%

Files
Source

# Copyright (c) 2018-2022 ISciences, LLC.
# All rights reserved.
#
# This software is licensed under the Apache License, Version 2.0 (the "License").
# You may not use this file except in compliance with the License. You may
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

.resultColNames <- function(...) {
  .resultColumns(...)$colname
}

#' Return column names to be used for summary operations
#'
#' @param value_names names of value raster layers
#' @param weight_names names of weighting raster layers
#' @param fun functions or names of summary operations
#' @param full_colnames return a complete column name even when there is no
#'                      ambiguity?
#' @param quantiles quantiles to use when \code{stat_names} contains \code{quantile}
#' @return character vector of column names
#' @keywords internal
.resultColumns <- function(value_names, weight_names, fun, full_colnames, quantiles=numeric(), unique_values=numeric(), colname_fun = NULL) {
  if (inherits(fun, 'standardGeneric')) {
    stat_names <- fun@generic[1]
  } else if (is.function(fun)) {
    stat_names <- 'fun'
  } else if (is.null(fun)) {
    stat_names <- ''
  } else {
    stat_names <- fun
  }

  quantile_index = which(stat_names == 'quantile')
  if (length(quantile_index) != 0) {
    stat_names <- splice(stat_names, quantile_index, rep('quantile', length(quantiles)))
  }

  for (stat in c('frac', 'weighted_frac')) {
    frac_index = which(stat_names == stat)
    if (length(frac_index) != 0) {
      stat_names <- splice(stat_names, frac_index, rep(stat, length(unique_values)))
    }
  }

  ind <- .valueWeightIndexes(length(value_names), length(weight_names))
  vn <- value_names[ind$values]
  wn <- weight_names[ind$weights]

  # determine all combinations of index and stat
  z <- expand.grid(index=seq_along(vn),
                   stat_name=stat_names, stringsAsFactors=FALSE)

  z$values <- vn[z$index]
  z$base_value <- NA_real_

  for (stat in c('frac', 'weighted_frac')) {
    ifrac <- which(z$stat_name == stat)
    z$base_value[ifrac] <- rep(unique_values, each = length(ifrac) / length(unique_values))
  }

  iquantile <- which(z$stat_name == 'quantile')
  z$base_value[iquantile] <- rep(quantiles, each = length(iquantile) / length(quantiles))

  if (is.null(wn)) {
    z$weights <- NA
  } else {
    z$weights <- wn[z$index]
  }
  z$weights[!.includeWeightInColName(z$stat_name)] <- NA

  if (is.null(colname_fun)) {
    colname_fun <- function(...) {
      .makeColname(full_colnames = full_colnames, ...)
    }
  }

  z$colname <- mapply(colname_fun,
                      fun_name = z$stat_name,
                      values = z$values,
                      weights = z$weights,
                      fun_value = z$base_value,
                      MoreArgs = list(nvalues = length(value_names),
                                      nweights = length(weight_names)),
                      USE.NAMES = FALSE)

  return(z)
}

.makeColname <- function(fun_name, values, weights, fun_value, full_colnames, nvalues, nweights) {
  # construct column names for each index, stat
  # add weight layer name only if layer is ambiguously weighted
  if (fun_name == 'quantile') {
    fun_component <- sprintf('q%02d', as.integer(100 * fun_value))
  } else if (fun_name %in% c('frac', 'weighted_frac')) {
    fun_component <- sprintf('%s_%s', fun_name, fun_value)
  } else {
    fun_component <- fun_name
  }

  ret <- fun_component
  if (full_colnames || nvalues > 1) {
    ret <- paste(ret, values, sep='.')
  }
  if ((!is.na(weights)) && ((full_colnames & nweights > 0) || nweights > 1)) {
    ret <- paste(ret, weights, sep='.')
  }

  return(ret)
}

.includeWeightInColName <- function(fun) {
  .isWeighted(fun) | fun == 'fun'
}

.isWeighted <- function(stat_name) {
  stat_name %in% c('weighted_mean', 'weighted_sum', 'weighted_frac')
}

#' Compute indexes for the value and weight layers that should be
#' processed together
#'
#' @param num_values number of layers in value raster
#' @param num_weights number of layers in weighting raster
#' @return list with \code{values} and \code{weights} elements
#'         providing layer indexes
#' @keywords internal
.valueWeightIndexes <- function(num_values, num_weights) {
  if (num_weights == 0) {
    vi <- seq_len(num_values)
    wi <- NA
  } else if (num_values == num_weights) {
    # process in parallel
    vi <- seq_len(num_values)
    wi <- seq_len(num_weights)
  } else if (num_values == 1 && num_weights > 1) {
    # recycle values
    vi <- rep.int(1, num_weights)
    wi <- seq_len(num_weights)
  } else if (num_values > 1 && num_weights == 1) {
    # recycle weights
    vi <- seq_len(num_values)
    wi <- rep.int(1, num_values)
  }

  list(values = vi, weights = wi)
}

.areaMethod <- function(crs_obj) {
  if (!(is.na(crs_obj)) && crs_obj$units_gdal == 'degree') {
    return('spherical')
  } else {
    return('cartesian')
  }
}

.validateFlag <- function(value, name) {
  if(!(is.logical(value) && length(value) == 1 && !is.na(value))) {
    stop(name, ' must be TRUE or FALSE')
  }
}

.validateNumericScalar <- function(value, name) {
  if (!(is.numeric(value) && length(value) == 1 && !is.na(value))) {
    stop(name, ' must be a single numeric value')
  }
}

.validateNumericScalarOrNA <- function(value, name) {
  if (!(is.numeric(value) && length(value) == 1)) {
    stop(name, ' must be a single numeric value')
  }
}

.validateUniqueNames <- function(x) {
  nm <- names(x)
  if (!is.null(nm)) {
    if (length(nm) != length(unique(nm))) {
      stop('names of input rasters must be unique')
    }
  }
}

# faster replacement for as.data.frame when input is a named list
# with equal-length columns
# from Advanced R, sec. 24.4.2
.quickDf <- function(lst) {
  class(lst) <- 'data.frame'
  attr(lst, 'row.names') <- .set_row_names(length(lst[[1]]))
  lst
}

.singleColumnToVector <- function(df) {
  if (ncol(df) == 1) {
    df[, 1]
  } else {
    df
  }
}

# Return the number of standard (non-...) arguments in a supplied function that
# do not have a default value. This is used to fail if the summary function
# provided by the user cannot accept arguments of values and weights.
.num_expected_args <- function(fun) {
  a <- formals(args(fun))
  a <- a[names(a) != '...']
  sum(sapply(a, nchar) == 0)
}

.startReading <- function(r) {
  if(inherits(r, 'BasicRaster')) {
    return(raster::readStart(r))
  } else if (inherits(r, 'SpatRaster')) {
    terra::readStart(r)
  }

  return(r)
}

.stopReading <- function(r) {
  if(inherits(r, 'BasicRaster')) {
    return(raster::readStop(r))
  } else if (inherits(r, 'SpatRaster')) {
    terra::readStop(r)
  }

  return(r)
}

.crs <- function(r) {
  if(inherits(r, 'BasicRaster')) {
    if (utils::packageVersion('raster') < numeric_version('3.5')) {
      return(raster::crs(r))
    } else {
      return(terra::crs(r))
    }
  } else if (inherits(r, 'SpatRaster')) {
    return(terra::crs(r))
  } else {
    stop('Unknown type: ', class(r))
  }
}

.setValues <- function(r, x) {
  if(inherits(r, 'BasicRaster')) {
    if (utils::packageVersion('raster') < numeric_version('3.5')) {
      raster::values(r) <- x
    } else {
      terra::values(r) <- x
    }
  } else if (inherits(r, 'SpatRaster')) {
    raster::values(r) <- x
  } else {
    stop('Unknown type: ', class(r))
  }

  return(r)
}

.numLayers <- function(r) {
  if(inherits(r, 'BasicRaster')) {
    return(raster::nlayers(r))
  } else if (inherits(r, 'SpatRaster')) {
    return(terra::nlyr(r))
  } else {
    stop('Unknown type: ', class(r))
  }
}

.isRaster <- function(r) {
  inherits(r, 'BasicRaster') | inherits(r, 'SpatRaster')
}

.xFromCol <- function(r, col) {
  if (inherits(r, 'BasicRaster')) {
    raster::xFromCol(r, col)
  } else if (inherits(r, 'SpatRaster')) {
    terra::xFromCol(r, col)
  } else {
    stop('Unknown type: ', class(r))
  }
}

.colFromX <- function(r, x) {
  if (inherits(r, 'BasicRaster')) {
    raster::colFromX(r, x)
  } else if (inherits(r, 'SpatRaster')) {
    terra::colFromX(r, x)
  } else {
    stop('Unknown type: ', class(r))
  }
}

.yFromRow <- function(r, row) {
  if (inherits(r, 'BasicRaster')) {
    raster::yFromRow(r, row)
  } else if (inherits(r, 'SpatRaster')) {
    terra::yFromRow(r, row)
  } else {
    stop('Unknown type: ', class(r))
  }
}

.rowFromY <- function(r, y) {
  if (inherits(r, 'BasicRaster')) {
    raster::rowFromY(r, y)
  } else if (inherits(r, 'SpatRaster')) {
    terra::rowFromY(r, y)
  } else {
    stop('Unknown type: ', class(r))
  }
}

.cellFromRowCol <- function(r, row, col) {
  if (inherits(r, 'BasicRaster')) {
    raster::cellFromRowCol(r, row, col)
  } else if (inherits(r, 'SpatRaster')) {
    terra::cellFromRowCol(r, row, col)
  } else {
    stop('Unknown type: ', class(r))
  }
}

.extent <- function(r) {
  if (inherits(r, 'BasicRaster')) {
    ex <- r@extent
    c(ex@xmin, ex@ymin, ex@xmax, ex@ymax)
  } else if (inherits(r, 'SpatRaster')) {
    ex <- terra::ext(r)
    c(ex$xmin, ex$ymin, ex$xmax, ex$ymax)
  } else {
    stop('Unknown type: ', class(r))
  }
}

.res <- function(r) {
  if (inherits(r, 'BasicRaster')) {
    raster::res(r)
  } else if (inherits(r, 'SpatRaster')) {
    terra::res(r)
  } else {
    stop('Unknown type: ', class(r))
  }
}

.getValuesBlock <- function(r, row, nrows, col, ncols) {
  if (inherits(r, 'BasicRaster')) {
    raster::getValuesBlock(r, row, nrows, col, ncols, format = 'm')
  } else if (inherits(r, 'SpatRaster')) {
    terra::readValues(r, row, nrows, col, ncols, mat = TRUE)
  } else {
    stop('Unknown type: ', class(r))
  }
}

.createProgress <- function(progress, n) {
  if (progress && n > 1) {
    pb <- utils::txtProgressBar(min = 0, max = n, initial=0, style=3)
    update_progress <- function() {
      i <- 1 + utils::getTxtProgressBar(pb)
      utils::setTxtProgressBar(pb, i)
      if (i == n) {
        close(pb)
      }
    }
  } else {
    update_progress <- function() {}
  }

  return(update_progress)
}

.isInMemory <- function(r) {
  if (inherits(r, 'BasicRaster')) {
    return(raster::inMemory(r))
  } else if (inherits(r, 'SpatRaster')) {
    return(terra::inMemory(r[[1]]))
  } else {
    stop('Unknown type: ', class(r))
  }
}

.netCDFBlockSize <- function(fname, varname) {
  nc <- NULL
  sz <- NA

  tryCatch({
    nc <- ncdf4::nc_open(fname)
    sz <- nc$var[[varname]]$chunksizes

    dim_index <- nc$var[[varname]]$dimids + 1L
    dim_names <- sapply(dim_index, function(i) nc$dim[[i]]$name)

    if (all(is.na(sz))) {
      # file is not compressed
      sz <- rep.int(1, length(dim_names))
    }

    names(sz) <- dim_names
  }, finally = {
    if (!is.null(nc)) {
      ncdf4::nc_close(nc)
    }
  })

  # flip dimensions 1 and 2 so we return row/col
  return(sz[c(2, 1, seq_along(sz)[-(1:2)])])
}

.blockSize <- function(r) {
  # set default return value in case file is uncompressed and has
  # has no block size, or we simply can't figure it out
  ret <- c(1, 1)

  if (inherits(r, 'BasicRaster')) {
    if (r[[1]]@file@driver == 'netcdf') {
      ret <- .netCDFBlockSize(r[[1]]@file@name, attr(r[[1]]@data, 'zvar'))
    } else if (r@file@driver == 'gdal') {
      ret <- c(r@file@blockrows, r@file@blockcols)
    }
  } else if (inherits(r, 'SpatRaster')) {
    ret <- terra::fileBlocksize(r)[1, ]
  }

  unname(ret)
}

.eagerLoad <- function(r, geoms, max_cells_in_memory, message_on_fail) {
  if (is.null(r)) {
    return(NULL)
  }

  cells_required <- .numCells(r, geoms)

  if (cells_required <= max_cells_in_memory) {
    box <- sf::st_bbox(geoms)
    geom_ext <- terra::ext(box[c('xmin', 'xmax', 'ymin', 'ymax')])

    if (!inherits(r, 'SpatRaster')) {
      # current CRAN version of terra (1.4-22) does not preserve
      # names on conversion (https://github.com/rspatial/terra/issues/430)
      nm <- names(r)
      r <- terra::rast(r)
      names(r) <- nm
    }

    overlap_ext <- terra::intersect(terra::ext(r), geom_ext)
    if (is.null(overlap_ext)) {
      # Extents do not overlap, and terra::crop will throw an error
      # if we try to crop. Return the input raster as-is; nothing will be
      # read from it anyway.
      return(r)
    }

    r <- terra::crop(r, geom_ext, snap = 'out')
  } else if (message_on_fail) {
    message('Cannot preload entire working area of ', cells_required,
            ' cells with max_cells_in_memory = ', max_cells_in_memory, '.',
            ' Raster values will be read for each feature individually.')
  }

  return(r)
}

.numCells <- function(r, g) {
  if (is.null(r)) {
    return(0)
  }

  box <- sf::st_bbox(g)

  top <- .rowFromY(r, box['ymax'])
  bottom <- .rowFromY(r, box['ymin'])
  left <- .colFromX(r, box['xmin'])
  right <- .colFromX(r, box['xmax'])

  if (is.na(top) && is.na(bottom)) {
    return(0L)
  }
  if (is.na(left) && is.na(right)) {
    return(0L)
  }
  if (is.na(top)) {
    top <- 1
  }
  if (is.na(bottom)) {
    bottom <- nrow(r)
  }
  if (is.na(left)) {
    left <- 1
  }
  if (is.na(right)) {
    right <- ncol(r)
  }

  return( (bottom - top + 1) * (right - left + 1) * .numLayers(r) )
}

splice <- function(x, i, replacement) {
  c(x[seq_along(x) < i],
    replacement,
    x[seq_along(x) > i])
}

# Copyright (c) 2018-2022 ISciences, LLC.
# All rights reserved.
#
# This software is licensed under the Apache License, Version 2.0 (the "License").
# You may not use this file except in compliance with the License. You may
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

setGeneric("exact_extract", function(x, y, ...)
	standardGeneric("exact_extract"))

#' Extract or summarize values from rasters
#'
#' Extracts the values of cells in a raster (`RasterLayer`, `RasterStack`
#' `RasterBrick`, or `SpatRaster`) that are covered by polygons in a
#' simple feature collection (`sf` or `sfc`) or `SpatialPolygonsDataFrame`.
#' Returns either a summary of the extracted values or the extracted values
#' themselves.
#'
#' @details
#' `exact_extract` extracts the values of cells in a raster that are covered
#' by polygonal features in a simple feature collection (`sf` or `sfc`) or
#' `SpatialPolygonDataFrame`, as well as the fraction or area of each cell that
#' is covered by the feature. Pixels covered by all parts of the polygon are
#' considered. If an (invalid) multipart polygon covers the same pixels more
#' than once, the pixel may have a coverage fraction greater than one.
#'
#' The function can either return pixel values directly to the caller, or can
#' return the result of a predefined summary operation or user-defined R
#' function applied to the values. These three approaches are described in the
#' subsections below.
#'
#' ## Returning extracted values directly
#'
#' If `fun` is not specified, `exact_extract` will return a list with
#' one data frame for each feature in the input feature collection. The data
#' frame will contain a column with cell values from each layer in the input
#' raster (and optional weighting raster) and a column indicating
#' the fraction or area of the cell that is covered by the polygon.
#'
#' If the input rasters have only one layer, the value and weight columns in the
#' data frame will be named `values` or `weights`. When the input rasters have
#' more than one layer, the columns will be named according to `names(x)` and
#' `names(weights)`. The column containing pixel coverage will be called
#' `coverage_fraction` when `coverage_area = FALSE`, or `coverage_area` when
#' `coverage_area = TRUE`. Additional columns can be added to the returned data
#' frames with the `include_area`, `include_cell`, and `include_xy` arguments.
#'
#' If the output data frames for multiple features are to be combined (e.g.,
#' with `rbind`), it may be useful to include identifying column(s) from the
#' input features in the returned data frames using `include_cols`.
#'
#' ## Predefined summary operations
#'
#' Often the individual pixel values are not needed; only one or more summary
#' statistics (e.g., mean, sum) is required for each feature. Common summary
#' statistics can be calculated by `exact_extract` directly using a predefined
#' summary operation. Where possible, this approach is advantageous because it
#' allows the package to calculate the statistics incrementally, avoiding the
#' need to store all pixel values in memory at the same time. This allows the
#' package to process arbitrarily large data with a small amount of memory. (The
#' `max_pixels_in_memory` argument can be used to fine-tune the amount of memory
#' made available to `exact_extract`.)
#'
#' To summarize pixel values using a predefined summary option, `fun` should be
#' set to a character vector of one or more operation names. If the input raster
#' has a single layer and a single summary operation is specified,
#' `exact_extract` will return a vector with the result of the summary operation
#' for each feature in the input. If the input raster has multiple layers, or if
#' multiple summary operations are specified, `exact_extract` will return a data
#' frame with a row for each feature and a column for each summary operation /
#' layer combination. (The `force_df` option can be used to always return a data
#' frame instead of a vector.)
#'
#' The following summary operations are supported:
#'
#'  * `min` - the minimum non-`NA` value in any raster cell wholly or
#'            partially covered by the polygon
#'  * `max` - the maximum non-`NA` value in any raster cell wholly or
#'            partially covered by the polygon
#'  * `count` - the sum of fractions of raster cells with non-`NA`
#'              values covered by the polygon
#'  * `sum`   - the sum of non-`NA` raster cell values, multiplied by
#'              the fraction of the cell that is covered by the polygon
#'  * `mean` - the mean cell value, weighted by the fraction of each cell
#'             that is covered by the polygon
#'  * `median` - the median cell value, weighted by the fraction of each cell
#'               that is covered by the polygon
#'  * `quantile` - arbitrary quantile(s) of cell values, specified in
#'                 `quantiles`, weighted by the fraction of each cell that is
#'                  covered by the polygon
#'  * `mode` - the most common cell value, weighted by the fraction of
#'             each cell that is covered by the polygon. Where multiple
#'             values occupy the same maximum number of weighted cells,
#'             the largest value will be returned.
#'  * `majority` - synonym for `mode`
#'  * `minority` - the least common cell value, weighted by the fraction
#'                 of each cell that is covered by the polygon. Where
#'                 multiple values occupy the same minimum number of
#'                 weighted cells, the smallest value will be returned.
#'  * `variety` - the number of distinct values in cells that are wholly or
#'                partially covered by the polygon.
#'  * `variance` - the population variance of cell values, weighted by the
#'                 fraction of each cell that is covered by the polygon.
#'  * `stdev` - the population standard deviation of cell values, weighted by
#'              the fraction of each cell that is covered by the polygon.
#'  * `coefficient_of_variation` - the population coefficient of variation of
#'                                 cell values, weighted by the fraction of each
#'                                 cell that is covered by the polygon.
#'  * `weighted_mean` - the mean cell value, weighted by the product of
#'                      the fraction of each cell covered by the polygon
#'                      and the value of a second weighting raster provided
#'                      as `weights`
#'  * `weighted_sum` - the sum of defined raster cell values, multiplied by
#'                     the fraction of each cell that is covered by the polygon
#'                     and the value of a second weighting raster provided
#'                     as `weights`
#'  * `frac` - returns one column for each possible value of `x`, with the
#'             the fraction of defined raster cells that are equal to that
#'             value.
#'  * `weighted_frac` - returns one column for each possible value of `x`,
#'                      with the fraction of defined cells that are equal
#'                      to that value, weighted by `weights.`
#'
#' In all of the summary operations, `NA` values in the the primary raster (`x`)
#' raster are ignored (i.e., `na.rm = TRUE`.) If `NA` values occur in the
#' weighting raster, the result of the weighted operation will be `NA`. `NA`
#' values in both `x` and `weights` can be replaced on-the-fly using the
#' `default_value` and `default_weight` arguments.
#'
#' ## User-defined summary functions
#'
#' If no predefined summary operation is suitable, a user-defined R function may
#' be provided as `fun`. The function will be called once for each feature and
#' must return either a single value or a data frame. The results of the
#' function for each feature will be combined and returned by `exact_extract`.
#'
#' The simplest way to write a summary function is to set
#' argument `summarize_df = TRUE`. (For backwards compatibility, this is not the
#' default.) In this mode, the summary function takes the signature
#' `function(df, ...)` where `df` is the same data frame that would be returned
#' by `exact_extract` with `fun = NULL`.
#'
#' With `summarize_df = FALSE`, the function must have the signature
#' `function(values, coverage_fractions, ...)` when weights are not used, and
#' `function(values, coverage_fractions, weights, ...)` when weights are used.
#' If the value and weight rasters each have a single layer, the function arguments
#' will be vectors; if either has multiple layers, the function arguments will
#' be data frames, with column names taken from the names of the value/weight
#' rasters. Values brought in through the `include_xy`, `include_area`,
#' `include_cell`, and `include_cols` arguments will be added to the `values`
#' data frame. For most applications, it is simpler to set `summarize_df = TRUE`
#' and work with all inputs in a single data frame.
#'
#' @param     x a `RasterLayer`, `RasterStack`, `RasterBrick`, or `SpatRaster`
#' @param     y a `sf`, `sfc`, `SpatialPolygonsDataFrame`, or `SpatialPolygons`
#'            object with polygonal geometries
#' @param     fun an optional function or character vector, as described below
#' @param     weights  a weighting raster to be used with the `weighted_mean`
#'                     and `weighted_sum` summary operations or a user-defined
#'                     summary function. When `weights` is set to `'area'`, the
#'                     cell areas of `x` will be calculated and used as weights.
#' @param     coverage_area  if `TRUE`, output pixel `coverage_area`
#'                           instead of `coverage_fraction`
#' @param     default_value  an optional value to use instead of `NA` in `x`
#' @param     default_weight an optional value to use instead of `NA` in `weights`
#' @param     quantiles   quantiles to be computed when `fun = 'quantile'`
#' @param     append_cols when `fun` is not `NULL`, an optional character vector
#'                        of columns from `y` to be included in returned data frame.
#' @param     force_df always return a data frame instead of a vector, even if
#'                     `x` has only one layer and `fun` has length 1
#' @param     summarize_df  pass values, coverage fraction/area, and weights to
#'                          `fun` as a single data frame instead of
#'                          separate arguments.
#' @param     full_colnames include the names of `x` and `weights` in
#'                          the names of the data frame for each feature, even if
#'                          `x` or `weights` has only one layer.
#'                          This is useful when the results of multiple
#'                          calls to `exact_extract` are combined with
#'                          `cbind`.
#' @param     colname_fun an optional function used to construct column names.
#'                        Should accept arguments `values` (name of value layer),
#'                        `weights` (name of weight layer), `fun_name` (value of
#'                        `fun`), `fun_value` (value associated with `fun`, for
#'                        `fun %in% c('quantile', 'frac', 'weighted_frac)`
#'                        `nvalues` (number of value layers), `weights`
#'                        (number of weight layers)
#' @param     include_area if `TRUE`, and `fun` is `NULL`, augment
#'                       the data frame for each feature with a column
#'                       for the cell area. If the units of the raster CRS are
#'                       degrees, the area in square meters will be calculated
#'                       based on a spherical approximation of Earth. Otherwise,
#'                       a Cartesian area will be calculated (and will be the
#'                       same for all pixels.) If `TRUE` and `fun` is
#'                       not `NULL`, add `area` to the data frame passed
#'                       to `fun` for each feature.
#' @param     include_cell if `TRUE`, and `fun` is `NULL`, augment
#'                       the data frame for each feature with a column
#'                       for the cell index (`cell`). If `TRUE` and
#'                       `fun` is not `NULL`, add `cell` to the
#'                       data frame passed to `fun` for each feature.
#' @param     include_cols an optional character vector of column names in
#'                         `y` to be added to the data frame for each
#'                         feature that is either returned (when `fun` is
#'                         `NULL`) or passed to `fun`.
#' @param     include_xy if `TRUE`, and `fun` is `NULL`, augment
#'                       the returned data frame for each feature with columns
#'                       for cell center coordinates (`x` and `y`). If
#'                       `TRUE` and `fun` is not `NULL`, add
#'                       `x` and `y` to the data frame passed to `fun`
#'                       for each feature.
#' @param     stack_apply   if `TRUE`, apply `fun` independently to
#'                          each layer or `x` (and its corresponding layer
#'                          of `weights`, if provided.) The number of
#'                          layers in `x` and `weights` must equal
#'                          each other or `1`, in which case the
#'                          single layer raster will be recycled.
#'                          If `FALSE`, apply `fun` to all layers of
#'                          `x` (and `weights`) simultaneously.
#' @param     max_cells_in_memory the maximum number of raster cells to load at
#'                                a given time when using a named summary operation
#'                                for `fun` (as opposed to a function defined using
#'                                R code). If a polygon covers more than `max_cells_in_memory`
#'                                raster cells, it will be processed in multiple chunks.
#' @param     grid_compat_tol     require value and weight grids to align within
#'                                `grid_compat_tol` times the smaller of the two
#'                                grid resolutions.
#' @param     progress if `TRUE`, display a progress bar during processing
#' @param     ... additional arguments to pass to `fun`
#' @return a vector, data frame, or list of data frames, depending on the type
#'         of `x` and the value of `fun` (see Details)
#' @examples
#' rast <- raster::raster(matrix(1:100, ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
#' poly <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))')
#'
#' # named summary operation on RasterLayer, returns vector
#' exact_extract(rast, poly, 'mean')
#'
#' # two named summary operations on RasterLayer, returns data frame
#' exact_extract(rast, poly, c('min', 'max'))
#'
#' # named summary operation on RasterStack, returns data frame
#' stk <- raster::stack(list(a=rast, b=sqrt(rast)))
#' exact_extract(stk, poly, 'mean')
#'
#' # named weighted summary operation, returns vector
#' weights <- raster::raster(matrix(runif(100), ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
#' exact_extract(rast, poly, 'weighted_mean', weights=weights)
#'
#' # custom summary function, returns vector
#' exact_extract(rast, poly, function(value, cov_frac) length(value[cov_frac > 0.9]))
#'
#' @name exact_extract
#' @md
NULL

.exact_extract <- function(x, y, fun=NULL, ...,
                           weights=NULL,
                           append_cols=NULL,
                           coverage_area=FALSE,
                           default_value=NA_real_,
                           default_weight=NA_real_,
                           include_area=FALSE,
                           include_cell=FALSE,
                           include_cols=NULL,
                           include_xy=FALSE,
                           force_df=FALSE,
                           full_colnames=FALSE,
                           stack_apply=FALSE,
                           summarize_df=FALSE,
                           quantiles=NULL,
                           progress=TRUE,
                           max_cells_in_memory=30000000,
                           grid_compat_tol=1e-3,
                           colname_fun=NULL
                           ) {
  area_weights <- is.character(weights) && length(weights) == 1 && weights == 'area'
  if (area_weights) {
    weights <- NULL
  }

  .validateFlag(coverage_area, 'coverage_area')
  .validateFlag(force_df, 'force_df')
  .validateFlag(full_colnames, 'full_colnames')
  .validateFlag(include_area, 'include_area')
  .validateFlag(include_cell, 'include_cell')
  .validateFlag(include_xy, 'include_xy')
  .validateFlag(progress, 'progress')
  .validateFlag(stack_apply, 'stack_apply')
  .validateFlag(summarize_df, 'summarize_df')
  .validateNumericScalar(max_cells_in_memory, 'max_cells_in_memory')
  .validateNumericScalar(grid_compat_tol, 'grid_compat_tol')
  .validateNumericScalarOrNA(default_value, 'default_value')
  .validateNumericScalarOrNA(default_weight, 'default_weight')

  .validateUniqueNames(x)
  .validateUniqueNames(weights)

  if(!is.null(append_cols)) {
    if (!inherits(y, 'sf')) {
      stop(sprintf('append_cols only supported for sf arguments (received %s)',
                   paste(class(y), collapse = ' ')))
    }

    force_df <- TRUE
  }

  if(!is.null(include_cols)) {
    if (!inherits(y, 'sf')) {
      stop(sprintf('include_cols only supported for sf arguments (received %s)',
                   paste(class(y), collapse = ' ')))
    }
  }

  if(sf::st_geometry_type(y, by_geometry = FALSE) == 'GEOMETRY') {
    if (!all(sf::st_dimension(y) == 2)) {
      stop("Features in sfc_GEOMETRY must be polygonal")
    }
  }

  if(!is.null(weights)) {
    if (!(.isRaster(weights) || weights == 'area')) {
      stop("Weights must be a Raster object or \"area\".")
    }

    if (is.character(fun) && !any(startsWith(fun, "weighted"))) {
      warning("Weights provided but no requested operations use them.")
    }

    if (inherits(weights, 'BasicRaster') && !is.na(sf::st_crs(x))) {
      if (is.na(sf::st_crs(weights))) {
        warning("No CRS specified for weighting raster; assuming it has the same CRS as the value raster.")
      } else if (sf::st_crs(x) != sf::st_crs(weights)) {
        stop("Weighting raster does not have the same CRS as value raster.")
      }
    }
  }

  analysis_crs <- sf::st_crs(x)
  if(is.na(sf::st_crs(x)) && !is.na(sf::st_crs(y))) {
    warning("No CRS specified for raster; assuming it has the same CRS as the polygons.")
    analysis_crs <- sf::st_crs(y)
  } else if(is.na(sf::st_crs(y)) && !is.na(sf::st_crs(x))) {
    warning("No CRS specified for polygons; assuming they have the same CRS as the raster.")
  } else if(sf::st_crs(x) != sf::st_crs(y)) {
    y <- sf::st_transform(y, sf::st_crs(x))
    warning("Polygons transformed to raster CRS (EPSG:", sf::st_crs(x)$epsg, ")")
  }
  if(area_weights || include_area || coverage_area) {
    area_method <- .areaMethod(analysis_crs)
  } else {
    area_method <- NULL
  }

  if (is.character(fun)) {
    if (length(fun) == 0) {
      stop("No summary operations provided.")
    }

    if (length(list(...)) > 0) {
      stop("exact_extract was called with a named summary operation that ",
           "does not accept additional arguments ...")
    }

    if (include_xy) {
      stop("include_xy must be FALSE for named summary operations")
    }

    if (include_cell) {
      stop("include_cell must be FALSE for named summary operations")
    }

    if (include_area) {
      stop("include_area must be FALSE for named summary operations")
    }

    if (!is.null(include_cols)) {
      stop("include_cols not supported for named_summary operations (see argument append_cols)")
    }

    if (summarize_df) {
      stop("summarize_df can only be used when `fun` is an R function")
    }
  } else if (is.function(fun)) {
    if (summarize_df) {
      if (.num_expected_args(fun) < 1) {
        stop("exact_extract was called with a function that does not appear to ",
             "be of the form `function(df, ...`).")
      }
    } else if (is.null(weights)) {
      if (.num_expected_args(fun) < 2) {
        stop("exact_extract was called with a function that does not appear to ",
             "be of the form `function(values, coverage_fractions, ...)`. If ",
             "the summary function should accept a single data frame argument, ",
             "set `summarize_df = TRUE`.")
      }
    } else if (.num_expected_args(fun) < 3) {
        stop("exact_extract was called with a function that does not appear to ",
             "be of the form `function(values, coverage_fractions, weights, ...)`.",
             "If the summary function should accept a single data frame argument, ",
             "set `summarize_df = TRUE`.")
    }
  } else if (is.null(fun)) {
    if (length(list(...)) > 0) {
      stop("Unexpected arguments: ", paste(names(list(...)), collapse = ','))
    }

    if (summarize_df) {
      stop("summarize_df can only be used when `fun` is an R function")
    }
    if (stack_apply) {
      stop("stack_apply can only be used when `fun` is a summary operation or function")
    }
    if (!is.null(append_cols)) {
      stop("append_cols can only be used when `fun` is a summary operation or function. See `include_cols`.")
    }
  } else {
    stop("fun must be a character vector, function, or NULL")
  }

  geoms <- sf::st_geometry(y)

  if (requireNamespace('terra', quietly = TRUE)) {
    strategies <- c('check_block_size', 'eager_load')
  } else {
    strategies = c()
  }

  x_orig <- NULL
  if ('eager_load' %in% strategies && length(geoms) > 1 && (!.isInMemory(x))) {
    # Eagerly load the entire area to be processed into memory. If the raster
    # block sizes are large, this potentially allows us to cache data in memory
    # that would not fit within the GDAL block cache (because we can cache an
    # area smaller than a single block). We only do this when the terra package
    # is available, because raster::crop throws an error when we set snap =
    # 'out'.
    if (include_cell) {
      # retain a reference to the original so that we can use it for include_cell
      x_orig <- x
    }
    x <- .eagerLoad(x, geoms, max_cells_in_memory, message_on_fail = progress)
    weights <- .eagerLoad(weights, geoms, max_cells_in_memory, message_on_fail = progress)
  }

  # At this point, if the data have not been preloaded into memory, either:
  # - we don't have the terra package available
  # - we can't fit the whole processing area in memory
  # Whether this is a problem depends on how the raster is chunked. If the
  # raster is poorly chunked, we can't do anything to improve the perormance,
  # but we try to detect the situation so we can alert the user.
  if (!.isInMemory(x)) {
    # - Do we have a RasterStack? Suggest using terra instead.
    if (inherits(x, 'RasterStack')) {
        message('exact_extract may perform poorly when extracting from a RasterStack. ',
                'It is recommended to use a SpatRaster from the terra package instead. ',
                'Convert using terra::rast(x)')
    }

    if (('check_block_size' %in% strategies) && inherits(x, 'SpatRaster')) {
      # Can we fit at least one block from every layer in GDAL's block cache?
      # If not, warn the user of expected poor performance.
      block_dim <- .blockSize(x)
      bytes_per_pixel <- 8L # terra has no datatype function?
      block_sz_mb <- prod(block_dim[1:2])* bytes_per_pixel / 1024 / 1024
      cache_sz_mb <- terra::gdalCache()

      min_cache_needed_mb <- ceiling(block_sz_mb) * .numLayers(x)


      if (min_cache_needed_mb >= cache_sz_mb) {
        message('Loading one raster block of data per layer requires approximately ',
                min_cache_needed_mb, ' MB but the GDAL block size cache is only ',
                cache_sz_mb, ' MB. This is likely to result in very slow performance. ',
                'It is recommended to increase the cache size using ',
                'terra::gdalCache(new_size), reduce the number of layers in the ',
                'input raster, or rewrite the input raster using a smaller chunk size.')
      }
    }
  }

  update_progress <- .createProgress(progress, length(geoms))
  tryCatch({
    x <- .startReading(x)
    weights <- .startReading(weights)

    if (is.character(fun)) {
      # Compute all stats in C++.
      # CPP_stats returns a matrix, which gets turned into a column by sapply
      # Results has one column per feature and one row per stat/raster layer
      if((is.null(weights) && !area_weights) && any(.isWeighted(fun))) {
        stop("Weighted stat requested but no weights provided.")
      }

      results <- lapply(sf::st_as_binary(geoms, EWKB=TRUE), function(wkb) {
        ret <- CPP_stats(x, weights, wkb, default_value, default_weight, coverage_area, area_method, fun, max_cells_in_memory, grid_compat_tol, quantiles)
        update_progress()
        return(ret)
      })

      if ('frac' %in% fun || 'weighted_frac' %in% fun) {
        unique_values <- unique(Reduce(c, lapply(results, function(r) attr(r, 'unique_values'))))
      } else {
        unique_values <- numeric()
      }

      result_cols <- .resultColumns(names(x), names(weights), fun, full_colnames, quantiles, unique_values, colname_fun = colname_fun)
      consistent_cols <- !(result_cols$stat_name %in% c('frac', 'weighted_frac'))

      # Construct an empty matrix that we can populate with stats returned
      # for each geometry, with 0 as a default where a geometry does not
      # return a stat (e.g., 'frac' for a specific value). Then convert the
      # filled matrix to a data frame. (This is faster than populating the
      # data frame directly, because the data frame assignment operator is
      # slow.)
      result_mat <- matrix(0.0, nrow = length(geoms), ncol = nrow(result_cols))
      for (i in seq_along(results)) {
        cols <- consistent_cols | result_cols$base_value %in% attr(results[[i]], 'unique_values')
        result_mat[i, cols] <- results[[i]]
      }
      colnames(result_mat) <- result_cols$colname

      result_df <- as.data.frame(result_mat)

      if (ncol(result_df) == 1 && !force_df) {
        return(result_df[[1]])
      }

      # drop duplicated columns (occurs when an unweighted stat is
      # requested alongside a weighted stat, with a stack of weights)
      result_df <- result_df[, unique(names(result_df)), drop = FALSE]

      if (!is.null(append_cols)) {
        result_df <- cbind(sf::st_drop_geometry(y[, append_cols]), result_df)
      }

      return(result_df)
    } else {
      num_values <- .numLayers(x)
      value_names <- names(x)
      if (.isRaster(weights)) {
        num_weights <- .numLayers(weights)
        weight_names <- names(weights)
      } else if (area_weights) {
        num_weights <- 1
        weight_names <- 'area'
        weights <- NULL
      } else {
        num_weights <- 0
        weight_names <- NULL
      }

      if (stack_apply || (num_values == 1 && num_weights <= 1)) {
        apply_layerwise <- TRUE

        if (num_values > 1 && num_weights > 1 && num_values != num_weights) {
          stop(sprintf("Can't apply function layerwise with stacks of %d value layers and %d layers", num_values, num_weights))
        }

        result_names <- .resultColNames(value_names, weight_names, fun, full_colnames, colname_fun = colname_fun)
        num_results <- max(num_weights, num_values)
        ind <- .valueWeightIndexes(num_values, num_weights)
      } else {
        apply_layerwise <- FALSE
      }

      # For R summary functions and data frame output, we avoid using
      # input layer names if there is only one value/weight layer
      if (length(value_names) == 1) {
        value_names <- 'value'
      }
      if (length(weight_names) == 1) {
        weight_names <- 'weight'
      }

      ret <- lapply(seq_along(geoms), function(feature_num) {
        wkb <- sf::st_as_binary(geoms[[feature_num]], EWKB=TRUE)

        if (!is.null(include_cols)) {
          include_col_values <- sf::st_drop_geometry(y[feature_num, include_cols])
        } else {
          include_col_values <- NULL
        }

        # only raise a disaggregation warning for the first feature
        warn_on_disaggregate <- feature_num == 1

        col_list <- CPP_exact_extract(x,
                                      x_orig,
                                      weights,
                                      wkb,
                                      default_value,
                                      default_weight,
                                      include_xy,
                                      include_cell,
                                      include_area,
                                      area_weights,
                                      coverage_area,
                                      area_method,
                                      include_col_values,
                                      value_names,
                                      weight_names,
                                      warn_on_disaggregate,
                                      grid_compat_tol)
        if (coverage_area) {
          coverage_col <- 'coverage_area'
        } else {
          coverage_col <- 'coverage_fraction'
        }

        if (!is.null(include_cols)) {
          # Replicate the include_cols vectors to be as long as the other columns,
          # so we can use quickDf
          nrow <- length(col_list[[coverage_col]])
          col_list[include_cols] <- lapply(col_list[include_cols], rep, nrow)
        }
        df <- .quickDf(col_list)

        update_progress()

        # No summary function? Return the whole data frame.
        if (is.null(fun)) {
          return(df)
        }

        # Summary function accepts a data frame? Pass it along.
        if (summarize_df && !apply_layerwise) {
          return(fun(df, ...))
        }

        # Break the data frame into components that we can pass
        # separately to the summary functions.
        included_cols_df <- df[, !(names(df) %in% c(value_names, weight_names, coverage_col)), drop = FALSE]
        vals_df <- df[, value_names, drop = FALSE]
        weights_df <- df[, weight_names, drop = FALSE]
        cov_fracs <- df[[coverage_col]]

        if (apply_layerwise) {
          result <- lapply(seq_len(num_results), function(i) {
            if (summarize_df) {
              # Pack everything into a single data frame and pass it to `fun`
              arg_df <- cbind(value = vals_df[[ind$values[i]]],
                              included_cols_df)
              if (num_weights > 0) {
                arg_df$weight <- weights_df[[ind$weights[i]]]
              }
              arg_df[[coverage_col]] <- df[[coverage_col]]

              return(fun(arg_df, ...))
            } else {
              # Pull values and weights out into vectors, unless we have
              # included columns (x/y/cell, etc.) in which case values
              # remain a data frame. Retained for backwards compat.

              vx <- vals_df[, ind$values[i]]
              if (ncol(included_cols_df) > 0) {
                vx <- cbind(data.frame(value = vx), included_cols_df)
              }
              if (num_weights == 0) {
                return(fun(vx, cov_fracs, ...))
              } else {
                return(fun(vx, cov_fracs, weights_df[, ind$weights[i]], ...))
              }
            }
          })

          if (num_results == 1) {
            return(result[[1]])
          }

          names(result) <- result_names

          return(.quickDf(result))
        } else {
          # Pass all layers to callback, to be handled together
          # Included columns (x/y/cell) are passed with the values.
          # Pass single-column data frames as vectors.
          vals_df <- .singleColumnToVector(cbind(vals_df, included_cols_df))
          weights_df <- .singleColumnToVector(weights_df)

          if (num_weights == 0) {
            return(fun(vals_df, cov_fracs, ...))
          } else {
            return(fun(vals_df, cov_fracs, weights_df, ...))
          }
        }
      })

      # if we have columns to append, iterate over the results and
      # append the columns to each, coercing into a data frame if
      # necessary. We need to iterate over the results before binding
      # them together because we don't know the legnth of each result
      # or even if their lengths are constant.
      if (!is.null(append_cols)) {
        for (i in seq_along(ret)) {
          if (!is.data.frame(ret[[i]])) {
            ret[[i]] = data.frame(result = ret[[i]])
          }

          if (nrow(ret[[i]]) >= 1) {
            append_row <- i
          } else {
            # cbinding row 0 cleanly brings in zero-length columns
            # of the correct names and types, in the correct positions
            append_row <- 0
          }

          ret[[i]] <- cbind(sf::st_drop_geometry(y[append_row, append_cols]),
                            ret[[i]],
                            row.names = NULL)
        }
      }

      if (!is.null(fun)) {
        if (all(sapply(ret, is.data.frame))) {
          # function returned a data frame for each polygon? rbind them
          if (requireNamespace('dplyr', quietly = TRUE)) {
            ret <- dplyr::bind_rows(ret) # handle column name mismatches
          } else {
            ret <- do.call(rbind, ret)
          }
        } else {
          # function returned something else; combine the somethings into
          # an array
          ret <- simplify2array(ret)

          if (force_df) {
            ret <- data.frame(result = ret)
          }
        }
      }

      return(ret)
    }
  }, finally={
    .stopReading(x)
    .stopReading(weights)
  })
}

# faster replacement for as.data.frame when input is a named list
# with equal-length columns
# from Advanced R, sec. 24.4.2
.quickDf <- function(lst) {
  class(lst) <- 'data.frame'
  attr(lst, 'row.names') <- .set_row_names(length(lst[[1]]))
  lst
}

.singleColumnToVector <- function(df) {
  if (ncol(df) == 1) {
    df[, 1]
  } else {
    df
  }
}

#' @import sf
#' @import raster
#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='Raster', y='sf'),
          .exact_extract)

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='Raster', y='SpatialPolygonsDataFrame'),
          function(x, y, ...) .exact_extract(x, sf::st_as_sf(y), ...))

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='Raster', y='SpatialPolygons'),
          function(x, y, ...) .exact_extract(x, sf::st_as_sf(y), ...))

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='Raster', y='sfc_MULTIPOLYGON'), .exact_extract)

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='Raster', y='sfc_POLYGON'), .exact_extract)

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='Raster', y='sfc_GEOMETRY'), .exact_extract)

# CRAN version of sf does not explicitly declare this class. If we do not do it
# ourselves, documentation is not generated for the `sfc_GEOMETRYCOLLECTION`
# overload.
setOldClass('sfc_GEOMETRYCOLLECTION')

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='Raster', y='sfc_GEOMETRYCOLLECTION'), .exact_extract)

#' @import sf
#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='SpatRaster', y='sf'),
          .exact_extract)

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='SpatRaster', y='SpatialPolygonsDataFrame'),
          function(x, y, ...) .exact_extract(x, sf::st_as_sf(y), ...))

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='SpatRaster', y='SpatialPolygons'),
          function(x, y, ...) .exact_extract(x, sf::st_as_sf(y), ...))

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='SpatRaster', y='sfc_MULTIPOLYGON'), .exact_extract)

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='SpatRaster', y='sfc_POLYGON'), .exact_extract)

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='SpatRaster', y='sfc_GEOMETRY'), .exact_extract)

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='SpatRaster', y='sfc_GEOMETRYCOLLECTION'), .exact_extract)

# Copyright (c) 2020-2022 ISciences, LLC.
# All rights reserved.
#
# This software is licensed under the Apache License, Version 2.0 (the "License").
# You may not use this file except in compliance with the License. You may
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

setGeneric("exact_resample", function(x, y, ...)
	standardGeneric("exact_resample"))

#' Resample a raster to a new grid
#'
#' @param x a \code{RasterLayer} or \code{SpatRaster} to be resampled
#' @param y a raster of the same class as \code{x} with a grid definition to
#'          which \code{x} should be resampled
#' @param fun a named summary operation or R function to be used for the resampling
#' @param coverage_area use cell coverage areas instead of coverage fractions
#'                      in \code{fun}
#' @return a resampled version of \code{x}, returned as a \code{RasterLayer} or
#'         \code{SpatRaster}, depending on the values of \code{x} and \code{y}
#'
#' @name exact_resample
NULL

.exact_resample <- function(x, y, fun, coverage_area = FALSE) {
  analysis_crs <- sf::st_crs(x)
  if (sf::st_crs(x) != sf::st_crs(y)) {
    if (is.na(sf::st_crs(x))) {
      warning("No CRS specified for source raster; assuming it has the same CRS as destination raster.")
      analysis_crs <- sf::st_crs(y)
    } else if (is.na(sf::st_crs(y))) {
      warning("No CRS specified for destination raster; assuming it has the same CRS as source raster.")
    } else {
      stop('Destination raster must have same CRS as source.')
    }
  }

  if (is.character(fun)) {
    if (length(fun) != 1) {
      stop("Only a single operation may be used for resampling.")
    }

    if (startsWith(fun, 'weighted')) {
      stop("Weighted operations cannot be used for resampling.")
    }

    if (.numLayers(x) != 1) {
      stop("Raster to be resampled must have a single layer.")
    }

    summary_stat <- fun
    summary_fun <- NULL
  } else {
    if (!is.function(fun)) {
      stop("fun must be a named summary operation or an R function")
    }

    if (.num_expected_args(fun) != 2) {
      stop("exact_extract was called with a function that does not appear to ",
           "be of the form `function(values, coverage_fractions)`.")
    }

    summary_stat <- NULL
    summary_fun <- fun
  }

  .validateFlag(coverage_area, 'coverage_area')

  area_method <- .areaMethod(analysis_crs)

  x <- .startReading(x)
  tryCatch({
    ret <- CPP_resample(x,
                        y,
                        summary_stat,
                        summary_fun,
                        coverage_area,
                        area_method)
    if (inherits(x, 'SpatRaster')) {
      terra::rast(ret)
    } else {
      ret
    }
  }, finally={
    .stopReading(x)
  })
}

#' @import sf
#' @import raster
#' @useDynLib exactextractr
#' @rdname exact_resample
#' @export
setMethod('exact_resample', signature(x='RasterLayer', y='RasterLayer'), .exact_resample)

#' @useDynLib exactextractr
#' @rdname exact_resample
#' @export
setMethod('exact_resample', signature(x='SpatRaster', y='SpatRaster'), .exact_resample)

# Copyright (c) 2018-2021 ISciences, LLC.
# All rights reserved.
#
# This software is licensed under the Apache License, Version 2.0 (the "License").
# You may not use this file except in compliance with the License. You may
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

setGeneric('coverage_fraction', function(x, y, crop=FALSE, ...)
	standardGeneric('coverage_fraction'))

.coverage_fraction <- function(x, y, crop) {
  if(is.na(sf::st_crs(x)) && !is.na(sf::st_crs(y))) {
    warning("No CRS specified for raster; assuming it has the same CRS as the polygons.")
  } else if(is.na(sf::st_crs(y)) && !is.na(sf::st_crs(x))) {
    warning("No CRS specified for polygons; assuming they have the same CRS as the raster.")
  } else if(sf::st_crs(x) != sf::st_crs(y)) {
    y <- sf::st_transform(y, sf::st_crs(x))
    warning("Polygons transformed to raster CRS (EPSG:", sf::st_crs(x)$epsg, ")")
  }

  lapply(sf::st_as_binary(y, EWKB=TRUE), function(wkb) {
    CPP_coverage_fraction(x, wkb, crop)
  })
}

#' Compute the fraction of raster cells covered by a polygon
#'
#' @param     x a (possibly empty) \code{RasterLayer} whose resolution and
#'            extent will be used for the generated \code{RasterLayer}.
#' @param     y a \code{sf} object with polygonal geometries
#' @param     crop if \code{TRUE}, each generated \code{RasterLayer} will be
#'                 cropped to the extent of its associated feature.
#' @return    a list with a \code{RasterLayer} for each feature in \code{y}.
#'            Values of the raster represent the fraction of each
#'            cell in \code{x} that is covered by \code{y}.
#' @examples
#' rast <- raster::raster(matrix(1:100, ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
#' poly <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))')
#'
#' cov_frac <- coverage_fraction(rast, poly)[[1]]
#' @name coverage_fraction
NULL

#' @import sf
#' @import raster
#' @useDynLib exactextractr
#' @rdname coverage_fraction
#' @export
setMethod('coverage_fraction', signature(x='RasterLayer', y='sf'), function(x, y, crop=FALSE) {
  coverage_fraction(x, sf::st_geometry(y), crop)
})

#' @rdname coverage_fraction
#' @export
setMethod('coverage_fraction', signature(x='RasterLayer', y='sfc_MULTIPOLYGON'), .coverage_fraction)

#' @rdname coverage_fraction
#' @export
setMethod('coverage_fraction', signature(x='RasterLayer', y='sfc_POLYGON'), .coverage_fraction)

#' @rdname coverage_fraction
#' @export
setMethod('coverage_fraction', signature(x='SpatRaster', y='sf'), function(x, y, crop=FALSE) {
  coverage_fraction(x, sf::st_geometry(y), crop)
})

#' @rdname coverage_fraction
#' @export
setMethod('coverage_fraction', signature(x='SpatRaster', y='sfc_MULTIPOLYGON'), .coverage_fraction)

#' @rdname coverage_fraction
#' @export
setMethod('coverage_fraction', signature(x='SpatRaster', y='sfc_POLYGON'), .coverage_fraction)

// Copyright (c) 2018-2021 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#include <Rcpp.h>

#include "geos_r.h"

#include "raster_utils.h"

#include "exactextract/src/raster.h"
#include "exactextract/src/raster_cell_intersection.h"

using exactextract::RasterView;
using exactextract::raster_cell_intersection;

// [[Rcpp::export]]
Rcpp::S4 CPP_coverage_fraction(Rcpp::S4 & rast,
                               const Rcpp::RawVector & wkb,
                               bool crop)
{
  try {
    GEOSAutoHandle geos;
    Rcpp::Environment xx = Rcpp::Environment::namespace_env("exactextractr");

    auto grid = make_grid(rast);
    auto coverage_fraction = raster_cell_intersection(grid, geos.handle, read_wkb(geos.handle, wkb).get());

    if (crop) {
      grid = coverage_fraction.grid();
    }
    RasterView<float> coverage_view(coverage_fraction, grid);

    Rcpp::NumericMatrix weights{static_cast<int>(grid.rows()),
                                static_cast<int>(grid.cols())};

    for (size_t i = 0; i < grid.rows(); i++) {
      for (size_t j = 0; j < grid.cols(); j++) {
        weights(i, j) = coverage_view(i, j);
        if (!crop && std::isnan(weights(i, j))) {
          weights(i, j) = 0;
        }
      }
    }

    if (rast.inherits("SpatRaster")) {
      Rcpp::Environment terra = Rcpp::Environment::namespace_env("terra");
      Rcpp::Function rastFn = terra["rast"];
      Rcpp::Function extFn = terra["ext"];
      Rcpp::Function crsFn = terra["crs"];

      Rcpp::S4 ext = extFn(grid.xmin(), grid.xmax(), grid.ymin(), grid.ymax());

      return rastFn(weights,
                    Rcpp::Named("ext") = ext,
                    Rcpp::Named("crs") = crsFn(rast));
    } else {
      Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster");
      Rcpp::Function rasterFn = raster["raster"];
      Rcpp::Function crsFn = xx[".crs"];

      return rasterFn(weights,
                      Rcpp::Named("xmn")=grid.xmin(),
                      Rcpp::Named("xmx")=grid.xmax(),
                      Rcpp::Named("ymn")=grid.ymin(),
                      Rcpp::Named("ymx")=grid.ymax(),
                      Rcpp::Named("crs")=crsFn(rast));
    }
  } catch (std::exception & e) {
    // throw predictable exception class
#ifdef __SUNPRO_CC
    // Rcpp::stop crashes CRAN Solaris build
    // https://github.com/RcppCore/Rcpp/issues/1159
    Rf_error(e.what());
    return R_NilValue;
#else
    Rcpp::stop(e.what());
#endif
  }
}

// Copyright (c) 2018-2020 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef EXACTEXTRACT_MATRIX_H
#define EXACTEXTRACT_MATRIX_H

#include <iomanip>
#include <iterator>
#include <memory>
#include <cstring>
#include <vector>

namespace exactextract {

template<typename T>
class Matrix {

    public:
        using value_type = T;

        Matrix(size_t rows, size_t cols) :
            m_rows{rows},
            m_cols{cols}
        {
            if (m_rows > 0 && m_cols > 0) {
                // new T[]() initializes to zero
                m_data = std::unique_ptr<T[]>(new T[m_rows * m_cols]());
            }
        }

        Matrix(size_t rows, size_t cols, T value) :
                m_rows{rows},
                m_cols{cols}
        {
            if (m_rows > 0 && m_cols > 0) {
                // new T[] does not initialize
                m_data = std::unique_ptr<T[]>(new T[m_rows * m_cols]);
            }

            std::fill(m_data.get(), m_data.get() + m_rows*m_cols, value);
        }

        explicit Matrix(const std::vector<std::vector<T>> & data) :
            m_rows{data.size()},
            m_cols{data[0].size()}
        {
            m_data = std::unique_ptr<T[]>(new T[m_rows*m_cols]());

            auto lastpos = m_data.get();
            for (auto& row : data) {
                lastpos = std::copy(row.begin(), row.end(), lastpos);
            }
        }

        Matrix(Matrix<T>&& m) noexcept :
                m_rows{m.rows()},
                m_cols{m.cols()}
        {
            m_data = std::move(m.m_data);
        }

        T& operator()(size_t row, size_t col) {
            check(row, col);
            return m_data[row*m_cols + col];
        }

        T operator()(size_t row, size_t col) const {
            check(row, col);
            return m_data[row*m_cols + col];
        }

        bool operator==(const Matrix<T> & other) const {
            if (m_rows != other.m_rows) {
                return false;
            }
            if (m_cols != other.m_cols) {
                return false;
            }

            return 0 == memcmp(m_data.get(), other.m_data.get(), m_rows*m_cols*sizeof(T));
        }

        void increment(size_t row, size_t col, const T & val) {
            check(row, col);
            m_data[row*m_cols + col] += val;
        }

        size_t rows() const { return m_rows; }
        size_t cols() const { return m_cols; }

        T* row(size_t row) {
            return &(m_data[row*m_cols]);
        }

        T* data() {
            return m_data.get();
        }

#ifdef MATRIX_CHECK_BOUNDS
        void check(size_t row, size_t col) const {
                if (row + 1 > m_rows) {
                    throw std::out_of_range("Row " + std::to_string(row) + " is out of range.");
                }
                if (col + 1 > m_cols) {
                    throw std::out_of_range("Col " + std::to_string(col) + " is out of range.");
                }
        }
#else
        void check(size_t, size_t) const {}
#endif

    private:
        std::unique_ptr<T[]> m_data;

        size_t m_rows;
        size_t m_cols;

};

template<typename T>
std::ostream& operator<<(std::ostream & os, const Matrix<T> & m) {
    for (size_t i = 0; i < m.rows(); i++) {
        for (size_t j = 0; j < m.cols(); j++) {
            if (m(i, j) != 0) {
                os << std::right << std::fixed << std::setw(10) << std::setprecision(6) <<
                    m(i, j) << " ";
            } else {
                os << "           ";
            }
        }
        os << std::endl;
    }

    return os;
}

}

#endif

// Copyright (c) 2018-2021 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef EXACTEXTRACT_RASTER_H
#define EXACTEXTRACT_RASTER_H

#include <array>
#include <cmath>
#include <limits>

#include "grid.h"
#include "matrix.h"

namespace exactextract {
    template<typename T>
    class AbstractRaster {
    public:
        using value_type = T;

        explicit AbstractRaster(const Grid<bounded_extent> & ex) :
            m_grid{ex},
            m_nodata{std::is_floating_point<T>::value ? std::numeric_limits<T>::quiet_NaN() : std::numeric_limits<T>::min()},
            m_has_nodata{false}
            {}

        AbstractRaster(const Grid<bounded_extent> & ex, const T& nodata_val) : m_grid{ex}, m_nodata{nodata_val}, m_has_nodata{true}  {}

        virtual ~AbstractRaster() = default;

        size_t rows() const {
            return m_grid.rows();
        }

        size_t cols() const {
            return m_grid.cols();
        }

        size_t size() const {
            return rows() * cols();
        }

        double xres() const {
            return m_grid.dx();
        }

        double yres() const {
            return m_grid.dy();
        }

        double xmin() const {
            return m_grid.xmin();
        }

        double ymin() const {
            return m_grid.ymin();
        }

        double xmax() const {
            return m_grid.xmax();
        }

        double ymax() const {
            return m_grid.ymax();
        }

        const Grid<bounded_extent>& grid() const {
            return m_grid;
        }

        virtual T operator()(size_t row, size_t col) const = 0;

        bool has_nodata() const { return m_has_nodata; }

        T nodata() const { return m_nodata; }

        bool get(size_t row, size_t col, T & val) const {
            val = operator()(row, col);

            if (m_has_nodata && val == m_nodata) {
                return false;
            }
            if (std::is_floating_point<T>::value && std::isnan(val)) {
                return false;
            }

            return true;
        }

        void set_nodata(const T & val) {
            m_has_nodata = true;
            m_nodata = val;
        }

        bool operator==(const AbstractRaster<T> & other) const {
            if (rows() != other.rows())
                return false;

            if (cols() != other.cols())
                return false;

            if (xres() != other.xres())
                return false;

            if (yres() != other.yres())
                return false;

            if (xmin() != other.xmin())
                return false;

            if (ymin() != other.ymin())
                return false;

            // Do the rasters differ in their definition of NODATA? If so, we need to do some more detailed
            // checking below.
            bool nodata_differs = has_nodata() != other.has_nodata() ||
                    (has_nodata() && other.has_nodata() && nodata() != other.nodata());

            for (size_t i = 0; i < rows(); i++) {
                for (size_t j = 0; j < cols(); j++) {
                    if(operator()(i, j) != other(i, j)) {
                        // Override default behavior of NAN != NAN
                        if (!std::isnan(operator()(i, j)) || !std::isnan(other(i, j))) {
                            return false;
                        }
                    } else if (nodata_differs && (operator()(i, j) == nodata() || other(i, j) == other.nodata())) {
                        // For data types that do not have NAN, or even for floating point types where the user has
                        // selected some other value to represent NODATA, we need to reverse a positive equality test
                        // where the value is considered to be NODATA in one raster but not the other.
                        return false;
                    }
                }
            }

            return true;
        }

        bool operator!=(const AbstractRaster<T> & other) const {
            return !(operator==(other));
        }

        class Iterator : public std::iterator<std::forward_iterator_tag, T> {
        public:
            Iterator(const AbstractRaster<T>* r, size_t i, size_t j) :
                    m_rast(r), m_row(i), m_col(j) {}

            const T& operator*() const {
                m_val = m_rast->operator()(m_row, m_col);
                return m_val;
            }

            Iterator& operator++() {
                m_col++;
                if (m_col == m_rast->cols()) {
                    m_col = 0;
                    m_row++;
                }
                return *this;
            }

            friend bool operator==(const Iterator& a,
                                   const Iterator& b) {
                return a.m_rast == b.m_rast && a.m_row == b.m_row && a.m_col == b.m_col;
            }

            friend bool operator!=(const Iterator& a,
                                   const Iterator& b) {
                return a.m_row != b.m_row || a.m_col != b.m_col || a.m_rast != b.m_rast;
            }

        private:
            const AbstractRaster<T>* m_rast;
            mutable T m_val;
            size_t m_row;
            size_t m_col;
        };

        Iterator begin() const {
            return Iterator(this, 0, 0);
        }

        Iterator end() const {
            return Iterator(this, rows(), 0);
        }
    private:
        Grid<bounded_extent> m_grid;
        T m_nodata;
        bool m_has_nodata;
    };

    template<typename T>
    class Raster : public AbstractRaster<T> {
    public:
        Raster(Matrix<T>&& values, const Box & box) : AbstractRaster<T>(
                Grid<bounded_extent>(
                        box,
                        (box.xmax - box.xmin) / values.cols(),
                        (box.ymax - box.ymin) / values.rows())),
                m_values{std::move(values)} {}

        Raster(Matrix<T>&& values, const Grid<bounded_extent> & g) :
            AbstractRaster<T>(g),
            m_values{std::move(values)} {}

        Raster(const Box & box, size_t nrow, size_t ncol) :
                AbstractRaster<T>(Grid<bounded_extent>(box, (box.xmax-box.xmin) / ncol, (box.ymax-box.ymin) / nrow)),
                m_values{nrow, ncol}
                {}

        explicit Raster(const Grid<bounded_extent> & ex) :
            AbstractRaster<T>(ex),
            m_values{ex.rows(), ex.cols()}
            {}

        Matrix<T>& data() {
            return m_values;
        }

        T& operator()(size_t row, size_t col) {
            return m_values(row, col);
        }

        T operator()(size_t row, size_t col) const override {
            return m_values(row, col);
        }

    private:
        Matrix<T> m_values;
    };

    template<typename T>
    class RasterView : public AbstractRaster<T>{
    public:
        // Construct a view of a raster r at an extent ex that is larger
        // and/or of finer resolution than r
        RasterView(const AbstractRaster<T> & r, Grid<bounded_extent> ex) :
            AbstractRaster<T>(ex),
                    m_raster{r},
                    m_x_off{0},
                    m_y_off{0},
                    m_rx{1},
                    m_ry{1} {
            if (!this->grid().empty()) {
                double disaggregation_factor_x = r.xres() / ex.dx();
                double disaggregation_factor_y = r.yres() / ex.dy();

                if (std::abs(disaggregation_factor_x - std::round(disaggregation_factor_x)) > 1e-6 ||
                    std::abs(disaggregation_factor_y - std::round(disaggregation_factor_y)) > 1e-6) {
                    throw std::runtime_error("Must construct view at resolution that is an integer multiple of original.");
                }

                if (disaggregation_factor_x < 0 || disaggregation_factor_y < 0) {
                    throw std::runtime_error("Must construct view at equal or higher resolution than original.");
                }

                m_x_off = static_cast<long>(std::round((ex.xmin() - r.xmin()) / ex.dx()));
                m_y_off = static_cast<long>(std::round((r.ymax() - ex.ymax()) / ex.dy()));
                m_rx = static_cast<size_t>(std::round(disaggregation_factor_x));
                m_ry = static_cast<size_t>(std::round(disaggregation_factor_y));
            }

            if (r.has_nodata()) {
                this->set_nodata(r.nodata());
            }
        }

        T operator()(size_t row, size_t col) const override {
            if (m_raster.grid().empty()) {
                return this->nodata();
            }

            if (m_x_off < 0 && static_cast<size_t>(-m_x_off) > col) {
                return this->nodata();
            }
            if (m_y_off < 0 && static_cast<size_t>(-m_y_off) > row) {
                return this->nodata();
            }

            size_t i0 = (row + m_y_off) / m_ry;
            size_t j0 = (col + m_x_off) / m_rx;

            if (i0 > m_raster.rows() - 1 || j0 > m_raster.cols() - 1) {
                return this->nodata();
            }

            return m_raster(i0, j0);
        }

    private:
        const AbstractRaster<T>& m_raster;

        long m_x_off;
        long m_y_off;
        size_t m_rx;
        size_t m_ry;
    };


    template<typename T>
    std::ostream& operator<<(std::ostream & os, const AbstractRaster<T> & m) {
        for (size_t i = 0; i < m.rows(); i++) {
            for (size_t j = 0; j < m.cols(); j++) {
                if (m(i, j) != 0) {
                    os << std::right << std::fixed << std::setw(10) << std::setprecision(6) <<
                       m(i, j) << " ";
                } else {
                    os << "           ";
                }
            }
            os << std::endl;
        }

        return os;
    }
}

#endif //EXACTEXTRACT_RASTER_H

// Copyright (c) 2018-2020 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

// [[Rcpp::plugins("cpp14")]]

#pragma once

#include <memory>

#include <Rcpp.h>

#include <geos_c.h>

using geom_ptr = std::unique_ptr<GEOSGeometry, std::function<void(GEOSGeometry*)>>;
using wkb_reader_ptr = std::unique_ptr<GEOSWKBReader, std::function<void(GEOSWKBReader*)>>;

// GEOS warning handler
static void geos_warn(const char* fmt, ...) {
  char buf[BUFSIZ] = { '\0' };

  va_list msg;
  va_start(msg, fmt);
  vsnprintf(buf, BUFSIZ*sizeof(char), fmt, msg);
  va_end(msg);

  Rcpp::Function warning("warning");
  warning(buf);
}

// GEOS error handler
static void geos_error(const char* fmt, ...) {
  char buf[BUFSIZ] = { '\0' };

  va_list msg;
  va_start(msg, fmt);
  vsnprintf(buf, BUFSIZ*sizeof(char), fmt, msg);
  va_end(msg);

  Rcpp::stop(buf);
}

// GEOSContextHandle wrapper to ensure finishGEOS is called.
struct GEOSAutoHandle {
  GEOSAutoHandle() {
    handle = initGEOS_r(geos_warn, geos_error);
  }

  ~GEOSAutoHandle() {
    finishGEOS_r(handle);
  }

  GEOSContextHandle_t handle;
};

// Return a smart pointer to a Geometry, given WKB input
static inline geom_ptr read_wkb(const GEOSContextHandle_t & context, const Rcpp::RawVector & wkb) {
  wkb_reader_ptr wkb_reader{ GEOSWKBReader_create_r(context), [context](GEOSWKBReader* r) { GEOSWKBReader_destroy_r(context, r); } };

  geom_ptr geom{ GEOSWKBReader_read_r(context,
                                      wkb_reader.get(),
                                      std::addressof(wkb[0]),
                                      wkb.size()),
                 [context](GEOSGeometry* g) { GEOSGeom_destroy_r(context, g); } };

  if (geom.get() == nullptr) {
    Rcpp::stop("Failed to parse WKB geometry");
  }

  return geom;
}

// Copyright (c) 2018-2022 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

// [[Rcpp::plugins("cpp14")]]
#include <Rcpp.h>

#include <memory>

#include "geos_r.h"
#include "raster_utils.h"

#include "s4_raster_source.h"

#include "exactextract/src/geos_utils.h"
#include "exactextract/src/grid.h"
#include "exactextract/src/matrix.h"
#include "exactextract/src/raster_cell_intersection.h"
#include "exactextract/src/raster_source.h"
#include "exactextract/src/raster_stats.h"

using exactextract::Box;
using exactextract::Matrix;
using exactextract::Grid;
using exactextract::subdivide;
using exactextract::bounded_extent;
using exactextract::Raster;
using exactextract::RasterView;
using exactextract::raster_cell_intersection;
using exactextract::RasterStats;
using exactextract::RasterSource;

#include <vector>
#include <string>

// [[Rcpp::export]]
Rcpp::List CPP_exact_extract(Rcpp::S4 & rast,
                             Rcpp::Nullable<Rcpp::S4> & rast_uncropped,
                             Rcpp::Nullable<Rcpp::S4> & weights,
                             const Rcpp::RawVector & wkb,
                             double default_value,
                             double default_weight,
                             bool include_xy,
                             bool include_cell_number,
                             bool include_area,
                             bool area_weights,
                             bool coverage_areas,
                             Rcpp::Nullable<Rcpp::CharacterVector> & p_area_method,
                             Rcpp::Nullable<Rcpp::List> & include_cols,
                             Rcpp::CharacterVector & src_names,
                             Rcpp::Nullable<Rcpp::CharacterVector> & p_weights_names,
                             bool warn_on_disaggregate,
                             double grid_compat_tol) {
  try {
    GEOSAutoHandle geos;
    Rcpp::Function names("names");

    auto grid = make_grid(rast);
    auto weights_grid = exactextract::Grid<bounded_extent>::make_empty();
    auto common_grid = grid;

    S4RasterSource rsrc(rast, default_value);
    int src_nlayers = get_nlayers(rast);

    std::unique_ptr<S4RasterSource> rweights;
    int weights_nlayers = 0;
    Rcpp::CharacterVector weights_names;

    std::string area_method;
    if (p_area_method.isNotNull()) {
      Rcpp::CharacterVector amethod = p_area_method.get();
      area_method = amethod[0];
    }

    if (weights.isNotNull()) {
      Rcpp::S4 weights_s4 = weights.get();
      weights_nlayers = get_nlayers(weights_s4);
      weights_grid = make_grid(weights_s4);

      common_grid = grid.common_grid(weights_grid, grid_compat_tol);

      rweights = std::make_unique<S4RasterSource>(weights_s4, default_weight);
      weights_names = p_weights_names.get();

      if (warn_on_disaggregate && (common_grid.dx() < grid.dx() || common_grid.dy() < grid.dy())) {
        Rcpp::warning("value raster implicitly disaggregated to match higher resolution of weights");
      }
    } else if (area_weights) {
      weights_nlayers = 1;
      weights_names = p_weights_names.get();
    }

    auto geom = read_wkb(geos.handle, wkb);

    auto bbox = exactextract::geos_get_box(geos.handle, geom.get());

    common_grid = common_grid.crop(bbox);

    auto coverage_fractions = raster_cell_intersection(common_grid, geos.handle, geom.get());
    auto& cov_grid = coverage_fractions.grid();

    // We can't construct an Rcpp::DataFrame with a variable number of columns
    // because of a bug in Rcpp (see https://github.com/RcppCore/Rcpp/pull/1099)
    // Instead, we build up a list, and then create a data frame from the list.
    // Profiling shows that this ends up in as.data.frame, which is slow.
    // Once Rcpp 1.0.6 is released, this can be reworked to be simpler and faster.
    Rcpp::List cols;

    Rcpp::NumericVector coverage_vec = as_vector(coverage_fractions);
    if (coverage_areas) {
        auto areas = get_area_raster(area_method, cov_grid);
        Rcpp::NumericVector area_vec = as_vector(*areas);
        coverage_vec = coverage_vec * area_vec;
    }
    Rcpp::LogicalVector covered = coverage_vec > 0;

    if (include_cols.isNotNull()) {
      Rcpp::List include_cols_list = include_cols.get();
      Rcpp::CharacterVector include_names = include_cols_list.attr("names");

      for (int i = 0; i < include_names.size(); i++) {
        std::string name(include_names[i]);
        cols[name] = include_cols_list[name];
      }
    }

    for (int i = 0; i < src_nlayers; i++) {
      auto values = rsrc.read_box(cov_grid.extent(), i);
      const NumericVectorRaster* r = static_cast<NumericVectorRaster*>(values.get());

      // TODO Perhaps extend this to preserve types (integer, logical.)
      // A bit challenging, since we don't know the type of the raster
      // until we call getValuesBlock, and even then we can't be sure
      // that the returned type is correct (as in a RasterStack with mixed types.)
      // Since R integers are only 32-bit, we are not going to lose data by
      // converting everything to numeric, although we pay a storage penalty.
      Rcpp::NumericVector value_vec = r->vec();
      if (grid.dx() != cov_grid.dx() || grid.dy() != cov_grid.dy() ||
          value_vec.size() != covered.size()) {
        // Transform values to same grid as coverage fractions
        RasterView<double> rt(*r, cov_grid);
        value_vec = as_vector(rt);
      }

      value_vec = value_vec[covered];
      cols[std::string(src_names[i])] = value_vec;
    }

    for (int i = 0; i < weights_nlayers; i++) {
      Rcpp::NumericVector weight_vec;
      if (area_weights) {
        auto weights = get_area_raster(area_method, cov_grid);
        weight_vec = as_vector(*weights);
      } else {
        auto values = rweights->read_box(cov_grid.extent(), i);
        const NumericVectorRaster* r = static_cast<NumericVectorRaster*>(values.get());
        weight_vec = r->vec();

        if (weights_grid.dx() != cov_grid.dx() || weights_grid.dy() != cov_grid.dy() ||
            weight_vec.size() != covered.size()) {
          // Transform weights to same grid as coverage fractions
          RasterView<double> rt (*r, cov_grid);

          weight_vec = as_vector(rt);
        }
      }

      weight_vec = weight_vec[covered];

      std::string colname(weights_names[i]);
      if (cols.containsElementNamed(colname.c_str())) {
        // append ".1" to the column name, to match the behavior of
        // data.frame(). We're safe to just add ".1" (instead of incrementing
        // .1 to .2, for example) because duplicated names within the values
        // or weight stack will already have been made unique by raster::stack()
        colname = colname + ".1"; // append .1
      }
      cols[colname] = weight_vec;
    }

    if (include_xy) {
      // Include xy values from whichever input raster has the higher resolution
      if (weights_nlayers > 0 && (weights_grid.dx() < grid.dx() || weights_grid.dy() < grid.dy())) {
        Rcpp::S4 weights_s4 = weights.get();
        cols["x"] = get_x_values(weights_s4, cov_grid)[covered];
        cols["y"] = get_y_values(weights_s4, cov_grid)[covered];
      } else {
        cols["x"] = get_x_values(rast, cov_grid)[covered];
        cols["y"] = get_y_values(rast, cov_grid)[covered];
      }
    }

    if (include_cell_number) {
      if (rast_uncropped.isNull()) {
        cols["cell"] = get_cell_numbers(rast, cov_grid)[covered];
      } else {
        Rcpp::S4 tmp = rast_uncropped.get();
        cols["cell"] = get_cell_numbers(tmp, cov_grid)[covered];
      }
    }

    if (include_area) {
      auto area_rast = get_area_raster(area_method, cov_grid);
      Rcpp::NumericVector area_vec = as_vector(*area_rast);
      cols["area"] = area_vec[covered];
    }

    if (coverage_areas) {
      cols["coverage_area"] = coverage_vec[covered];
    } else {
      cols["coverage_fraction"] = coverage_vec[covered];
    }

    return cols;
  } catch (std::exception & e) {
    // throw predictable exception class
#ifdef __SUNPRO_CC
    // Rcpp::stop crashes CRAN Solaris build
    // https://github.com/RcppCore/Rcpp/issues/1159
    Rf_error(e.what());
    return R_NilValue;
#else
    Rcpp::stop(e.what());
#endif
  }
}

enum class WeightingMethod {
  NONE,
  RASTER,
  AREA
};


static int get_num_stats(const Rcpp::StringVector & stats,
                  const std::size_t num_quantiles,
                  const std::size_t num_unique_values) {
  int num_stats = 0;
  for (const auto & stat : stats) {
      if (stat == std::string("frac") ||
          stat == std::string("weighted_frac")) {
        num_stats += static_cast<int>(num_unique_values);
      } else if (stat == std::string("quantile")) {
        num_stats += static_cast<int>(num_quantiles);
      } else {
        num_stats += 1;
      }
  }

  return num_stats;
}

// Return a matrix with one row per stat and one row per raster layer
// [[Rcpp::export]]
Rcpp::NumericMatrix CPP_stats(Rcpp::S4 & rast,
                              Rcpp::Nullable<Rcpp::S4> weights,
                              const Rcpp::RawVector & wkb,
                              double default_value,
                              double default_weight,
                              bool coverage_areas,
                              Rcpp::Nullable<Rcpp::CharacterVector> & p_area_method,
                              const Rcpp::StringVector & stats,
                              int max_cells_in_memory,
                              double grid_compat_tol,
                              const Rcpp::Nullable<Rcpp::NumericVector> & quantiles) {
  try {
    GEOSAutoHandle geos;

    if (max_cells_in_memory < 1) {
      Rcpp::stop("Invalid value for max_cells_in_memory: %d", max_cells_in_memory);
    }

    int nlayers = get_nlayers(rast);

    S4RasterSource rsrc(rast, default_value);

    std::unique_ptr<S4RasterSource> rweights;
    std::string area_method;
    if (p_area_method.isNotNull()) {
      area_method = ((Rcpp::CharacterVector) p_area_method.get())[0];
    }

    WeightingMethod weighting = WeightingMethod::NONE;
    int nweights = 0;
    if (weights.isNotNull()) {
      weighting = WeightingMethod::RASTER;
      Rcpp::S4 weights_s4 = weights.get();
      nweights = get_nlayers(weights_s4);

      if (nlayers > 1 && nweights > 1 && nlayers != nweights) {
        Rcpp::stop("Incompatible number of layers in value and weighting rasters");
      }

      rweights = std::make_unique<S4RasterSource>(weights_s4, default_weight);
    } else if (p_area_method.isNotNull()) {
      weighting = WeightingMethod::AREA;
      nweights = 1;
    }

    auto geom = read_wkb(geos.handle, wkb);
    auto bbox = exactextract::geos_get_box(geos.handle, geom.get());

    auto grid = weighting == WeightingMethod::RASTER ?
      rsrc.grid().common_grid(rweights->grid(), grid_compat_tol) : rsrc.grid();

    bool disaggregated = (grid.dx() < rsrc.grid().dx() || grid.dy() < rsrc.grid().dy());

    bool store_values = false;
    bool calc_value_set = false;
    std::size_t num_unique_values = 0;
    std::size_t num_quantiles = 0;

    for (const auto & stat : stats) {
      store_values = store_values || requires_stored_values(stat);

      if (disaggregated && (stat == std::string("count") ||
                            stat == std::string("sum"))) {
        Rcpp::stop("Cannot compute 'count' or 'sum' when value raster is disaggregated to resolution of weights.");
      }

      if (stat == std::string("frac") ||
          stat == std::string("weighted_frac")) {
        calc_value_set = true;
      }

      if (stat == std::string("quantile")) {
        if (quantiles.isNotNull()) {
          Rcpp::NumericVector qvec = quantiles.get();
          num_quantiles = qvec.size();
        }

        if (num_quantiles == 0) {
          Rcpp::stop("Quantiles not specified.");
        }
      }
    }

    int nresults = std::max(nlayers, nweights);

    std::vector<RasterStats<double>> raster_stats;
    raster_stats.reserve(nresults);
    for (int i = 0; i < nresults; i++) {
      raster_stats.emplace_back(store_values);
    }

    if (bbox.intersects(grid.extent())) {
      auto cropped_grid = grid.crop(bbox);

      for (const auto &subgrid : subdivide(cropped_grid, max_cells_in_memory)) {
        auto coverage_fraction = raster_cell_intersection(subgrid, geos.handle, geom.get());
        if (coverage_areas) {
          auto areas = get_area_raster(area_method, subgrid);
          for (size_t i = 0; i < coverage_fraction.rows(); i++) {
            for (size_t j = 0; j < coverage_fraction.cols(); j++) {
              coverage_fraction(i, j) = coverage_fraction(i, j) * (*areas)(i, j);
            }
          }
        }

        auto& cov_grid = coverage_fraction.grid();

        if (!cov_grid.empty()) {
          if (weighting == WeightingMethod::NONE) {
            for (int i = 0; i < nlayers; i++) {
              auto values = rsrc.read_box(cov_grid.extent(), i);
              raster_stats[i].process(coverage_fraction, *values);
            }
          } else if (weighting == WeightingMethod::AREA) {
            auto weights = get_area_raster(area_method, cov_grid);

            for (int i = 0; i < nlayers; i++) {
              auto values = rsrc.read_box(cov_grid.extent(), i);
              raster_stats[i].process(coverage_fraction, *values, *weights);
            }
          } else {
            if (nlayers > nweights) {
              // recycle weights
              auto weights = rweights->read_box(cov_grid.extent(), 0);

              for (int i = 0; i < nlayers; i++) {
                auto values = rsrc.read_box(cov_grid.extent(), i);
                raster_stats[i].process(coverage_fraction, *values, *weights);
              }
            } else if (nweights > nlayers) {
              // recycle values
              auto values = rsrc.read_box(cov_grid.extent(), 0);

              for (int i = 0; i < nweights; i++) {
                auto weights = rweights->read_box(cov_grid.extent(), i);
                raster_stats[i].process(coverage_fraction, *values, *weights);
              }
            } else {
              // process values and weights in parallel
              for (int i = 0; i < nlayers; i++) {
                auto values = rsrc.read_box(cov_grid.extent(), i);
                auto weights = rweights->read_box(cov_grid.extent(), i);

                raster_stats[i].process(coverage_fraction, *values, *weights);
              }
            }
          }
        }
      }
    }


    std::set<double> value_set;
    if (calc_value_set) {
      for (const auto& rs : raster_stats) {
        for (const auto& value : rs) {
          value_set.insert(value);
        }
      }
      num_unique_values = value_set.size();
    }

    auto stat_result_size = get_num_stats(stats, num_quantiles, num_unique_values);
    Rcpp::NumericMatrix stat_results = Rcpp::no_init(nresults, stat_result_size);

    if (calc_value_set) {
      Rcpp::NumericVector unique_values(value_set.begin(), value_set.end());
      stat_results.attr("unique_values") = unique_values;
    }

    // rows (j) represent layers
    // cols (i) represent stats
    for (int j = 0; j < nresults; j++) {
      const auto& rs = raster_stats[j];

      int i = 0;
      for(const auto& stat : stats) {
        if (stat == std::string("mean")) stat_results(j, i++) = rs.mean();

        else if (stat == std::string("sum")) stat_results(j, i++) = rs.sum();
        else if (stat == std::string("count")) stat_results(j, i++) = rs.count();

        else if (stat == std::string("min")) stat_results(j, i++) = rs.min().value_or(NA_REAL);
        else if (stat == std::string("max")) stat_results(j, i++) = rs.max().value_or(NA_REAL);

        else if (stat == std::string("median")) stat_results(j, i++) = rs.quantile(0.5).value_or(NA_REAL);

        else if (stat == std::string("mode")) stat_results(j, i++) = rs.mode().value_or(NA_REAL);
        else if (stat == std::string("majority")) stat_results(j, i++) = rs.mode().value_or(NA_REAL);
        else if (stat == std::string("minority")) stat_results(j, i++) = rs.minority().value_or(NA_REAL);

        else if (stat == std::string("variety")) stat_results(j, i++) = rs.variety();
        else if (stat == std::string("weighted_mean")) stat_results(j, i++) = rs.weighted_mean();
        else if (stat == std::string("weighted_sum")) stat_results(j, i++) = rs.weighted_sum();

        else if (stat == std::string("variance")) stat_results(j, i++) = rs.variance();
        else if (stat == std::string("stdev")) stat_results(j, i++) = rs.stdev();
        else if (stat == std::string("coefficient_of_variation")) stat_results(j, i++) = rs.coefficient_of_variation();

        else if (stat == std::string("quantile")) {
          Rcpp::NumericVector qvec = quantiles.get();
          for (double q : qvec) {
            stat_results(j, i++) = rs.quantile(q).value_or(NA_REAL);
          }
        }

        else if (stat == std::string("frac")) {
          for (double v : value_set) {
            stat_results(j, i++) = rs.frac(v).value_or(0);
          }
        }

        else if (stat == std::string("weighted_frac")) {
          for (double v : value_set) {
            stat_results(j, i++) = rs.weighted_frac(v).value_or(0);
          }
        }

        else Rcpp::stop("Unknown stat: " + stat);
      }
    }

    // FIXME set dimnames here

    return stat_results;
  } catch (std::exception & e) {
    // throw predictable exception class
#ifdef __SUNPRO_CC
    // Rcpp::stop crashes CRAN Solaris build
    // https://github.com/RcppCore/Rcpp/issues/1159
    Rf_error(e.what());
    return R_NilValue;
#else
    Rcpp::stop(e.what());
#endif
  }
}

// Copyright (c) 2018-2019 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef EXACTEXTRACT_BOX_H
#define EXACTEXTRACT_BOX_H

#include "coordinate.h"
#include "crossing.h"
#include "side.h"

#include <limits>

namespace exactextract {

    struct Box {
        double xmin;
        double ymin;
        double xmax;
        double ymax;

        Box(double p_xmin, double p_ymin, double p_xmax, double p_ymax) :
                xmin{p_xmin},
                ymin{p_ymin},
                xmax{p_xmax},
                ymax{p_ymax} {}

        static Box maximum_finite() {
            return {
                std::numeric_limits<double>::lowest(),
                std::numeric_limits<double>::lowest(),
                std::numeric_limits<double>::max(),
                std::numeric_limits<double>::max()
            };
        }

        static Box make_empty() {
            return {0, 0, 0, 0};
        }

        double width() const {
            return xmax - xmin;
        }

        double height() const {
            return ymax - ymin;
        }

        double area() const {
            return width() * height();
        }

        double perimeter() const {
            return 2 * width() + 2 * height();
        }

        bool intersects(const Box & other) const {
            if (other.ymin > ymax)
                return false;
            if (other.ymax < ymin)
                return false;
            if (other.xmin > xmax)
                return false;
            if (other.xmax < xmin)
                return false;

            return true;
        }

        Box intersection(const Box & other) const {
            return {
                std::max(xmin, other.xmin),
                std::max(ymin, other.ymin),
                std::min(xmax, other.xmax),
                std::min(ymax, other.ymax)
            };
        }

        Box translate(double dx, double dy) const {
            return {
                xmin + dx,
                ymin + dy,
                xmax + dx,
                ymax + dy
            };
        }

        Coordinate upper_left() const {
            return Coordinate{xmin, ymax};
        }

        Coordinate upper_right() const {
            return Coordinate{xmax, ymax};
        }

        Coordinate lower_left() const {
            return Coordinate{xmin, ymin};
        }

        Coordinate lower_right() const {
            return Coordinate{xmax, ymin};
        }

        Side side(const Coordinate &c) const;

        Crossing crossing(const Coordinate &c1, const Coordinate &c2) const;

        bool empty() const {
            return xmin >= xmax || ymin >= ymax;
        }

        Box expand_to_include(const Box & other) const {
            if (empty()) {
                return other;
            }

            if (other.empty()) {
                return *this;
            }

            return { std::min(xmin, other.xmin),
                     std::min(ymin, other.ymin),
                     std::max(xmax, other.xmax),
                     std::max(ymax, other.ymax) };
        }

        bool contains(const Box &b) const;

        bool contains(const Coordinate &c) const;

        bool strictly_contains(const Coordinate &c) const;

        bool operator==(const Box& other) const {
            return xmin == other.xmin && xmax == other.xmax && ymin == other.ymin && ymax == other.ymax;
        }

        friend std::ostream &operator<<(std::ostream &os, const Box &c);
    };

    std::ostream &operator<<(std::ostream &os, const Box &c);

}

#endif

// Copyright (c) 2021 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef EXACTEXTRACT_RASTER_AREA_H
#define EXACTEXTRACT_RASTER_AREA_H

#include "raster.h"

namespace exactextract {
    template<typename T>
    class CartesianAreaRaster : public AbstractRaster<T> {
    public:
        explicit CartesianAreaRaster(const Grid<bounded_extent> & ex) :
                AbstractRaster<T>(ex),
                m_area(grid_cell(AbstractRaster<T>::grid(), 0, 0).area())
        {}

        T operator()(size_t row, size_t column) const override {
            (void) row;
            (void) column;
            return m_area;
        }
    private:
        double m_area;
    };

    template<typename T>
    class SphericalAreaRaster : public AbstractRaster<T> {
    public:
        explicit SphericalAreaRaster(const Grid<bounded_extent> &ex) :
                AbstractRaster<T>(ex),
                m_areas(ex.rows())
        {
            const auto& g = AbstractRaster<T>::grid();

            double dlon = g.dx();
            double dlat = g.dy();

            for (size_t i = 0; i < g.rows(); i++) {
                double y = g.y_for_row(i);
                double ymin = y - 0.5 * dlat;
                double ymax = y + 0.5 * dlat;

                m_areas[i] = EARTH_RADIUS_SQ * PI_180 * std::abs(std::sin(ymin * PI_180) - std::sin(ymax * PI_180)) * dlon;
            }
        }

        T operator()(size_t row, size_t column) const override {
            (void) column;
            return m_areas[row];
        }

    private:
        static constexpr double EARTH_RADIUS = 6378137;
        static constexpr double EARTH_RADIUS_SQ = EARTH_RADIUS * EARTH_RADIUS;
        static constexpr double PI_180 = 3.141592653589793238462643383279502884197169399375105 / 180;

        std::vector<double> m_areas;
    };
}

#endif //EXACTEXTRACT_RASTER_AREA_H

// Copyright (c) 2018-2022 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef EXACTEXTRACT_RASTER_STATS_H
#define EXACTEXTRACT_RASTER_STATS_H

#include <algorithm>
#include <iostream>
#include <limits>
#include <unordered_map>

#include "raster_cell_intersection.h"
#include "weighted_quantiles.h"
#include "variance.h"

#include "../vend/optional.hpp"

namespace exactextract {

    template<typename T>
    class RasterStats {

    public:

1		# Copyright (c) 2018-2022 ISciences, LLC.
2		# All rights reserved.
3		#
4		# This software is licensed under the Apache License, Version 2.0 (the "License").
5		# You may not use this file except in compliance with the License. You may
6		# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		#
8		# Unless required by applicable law or agreed to in writing, software
9		# distributed under the License is distributed on an "AS IS" BASIS,
10		# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		# See the License for the specific language governing permissions and
12		# limitations under the License.
13
14		.resultColNames <- function(...) {
15	108x	.resultColumns(...)$colname
16		}
17
18		#' Return column names to be used for summary operations
19		#'
20		#' @param value_names names of value raster layers
21		#' @param weight_names names of weighting raster layers
22		#' @param fun functions or names of summary operations
23		#' @param full_colnames return a complete column name even when there is no
24		#' ambiguity?
25		#' @param quantiles quantiles to use when \code{stat_names} contains \code{quantile}
26		#' @return character vector of column names
27		#' @keywords internal
28		.resultColumns <- function(value_names, weight_names, fun, full_colnames, quantiles=numeric(), unique_values=numeric(), colname_fun = NULL) {
29	246x	if (inherits(fun, 'standardGeneric')) {
30	22x	stat_names <- fun@generic[1]
31	224x	} else if (is.function(fun)) {
32	32x	stat_names <- 'fun'
33	192x	} else if (is.null(fun)) {
34	42x	stat_names <- ''
35		} else {
36	150x	stat_names <- fun
37		}
38
39	246x	quantile_index = which(stat_names == 'quantile')
40	246x	if (length(quantile_index) != 0) {
41	6x	stat_names <- splice(stat_names, quantile_index, rep('quantile', length(quantiles)))
42		}
43
44	246x	for (stat in c('frac', 'weighted_frac')) {
45	492x	frac_index = which(stat_names == stat)
46	492x	if (length(frac_index) != 0) {
47	6x	stat_names <- splice(stat_names, frac_index, rep(stat, length(unique_values)))
48		}
49		}
50
51	246x	ind <- .valueWeightIndexes(length(value_names), length(weight_names))
52	246x	vn <- value_names[ind$values]
53	246x	wn <- weight_names[ind$weights]
54
55		# determine all combinations of index and stat
56	246x	z <- expand.grid(index=seq_along(vn),
57	246x	stat_name=stat_names, stringsAsFactors=FALSE)
58
59	246x	z$values <- vn[z$index]
60	246x	z$base_value <- NA_real_
61
62	246x	for (stat in c('frac', 'weighted_frac')) {
63	492x	ifrac <- which(z$stat_name == stat)
64	492x	z$base_value[ifrac] <- rep(unique_values, each = length(ifrac) / length(unique_values))
65		}
66
67	246x	iquantile <- which(z$stat_name == 'quantile')
68	246x	z$base_value[iquantile] <- rep(quantiles, each = length(iquantile) / length(quantiles))
69
70	246x	if (is.null(wn)) {
71	160x	z$weights <- NA
72		} else {
73	86x	z$weights <- wn[z$index]
74		}
75	246x	z$weights[!.includeWeightInColName(z$stat_name)] <- NA
76
77	246x	if (is.null(colname_fun)) {
78	245x	colname_fun <- function(...) {
79	387x	.makeColname(full_colnames = full_colnames, ...)
80		}
81		}
82
83	246x	z$colname <- mapply(colname_fun,
84	246x	fun_name = z$stat_name,
85	246x	values = z$values,
86	246x	weights = z$weights,
87	246x	fun_value = z$base_value,
88	246x	MoreArgs = list(nvalues = length(value_names),
89	246x	nweights = length(weight_names)),
90	246x	USE.NAMES = FALSE)
91
92	246x	return(z)
93		}
94
95		.makeColname <- function(fun_name, values, weights, fun_value, full_colnames, nvalues, nweights) {
96		# construct column names for each index, stat
97		# add weight layer name only if layer is ambiguously weighted
98	387x	if (fun_name == 'quantile') {
99	13x	fun_component <- sprintf('q%02d', as.integer(100 * fun_value))
100	374x	} else if (fun_name %in% c('frac', 'weighted_frac')) {
101	26x	fun_component <- sprintf('%s_%s', fun_name, fun_value)
102		} else {
103	348x	fun_component <- fun_name
104		}
105
106	387x	ret <- fun_component
107	387x	if (full_colnames \|\| nvalues > 1) {
108	137x	ret <- paste(ret, values, sep='.')
109		}
110	387x	if ((!is.na(weights)) && ((full_colnames & nweights > 0) \|\| nweights > 1)) {
111	54x	ret <- paste(ret, weights, sep='.')
112		}
113
114	387x	return(ret)
115		}
116
117		.includeWeightInColName <- function(fun) {
118	246x	.isWeighted(fun) \| fun == 'fun'
119		}
120
121		.isWeighted <- function(stat_name) {
122	334x	stat_name %in% c('weighted_mean', 'weighted_sum', 'weighted_frac')
123		}
124
125		#' Compute indexes for the value and weight layers that should be
126		#' processed together
127		#'
128		#' @param num_values number of layers in value raster
129		#' @param num_weights number of layers in weighting raster
130		#' @return list with \code{values} and \code{weights} elements
131		#' providing layer indexes
132		#' @keywords internal
133		.valueWeightIndexes <- function(num_values, num_weights) {
134	330x	if (num_weights == 0) {
135	223x	vi <- seq_len(num_values)
136	223x	wi <- NA
137	107x	} else if (num_values == num_weights) {
138		# process in parallel
139	90x	vi <- seq_len(num_values)
140	90x	wi <- seq_len(num_weights)
141	17x	} else if (num_values == 1 && num_weights > 1) {
142		# recycle values
143	8x	vi <- rep.int(1, num_weights)
144	8x	wi <- seq_len(num_weights)
145	9x	} else if (num_values > 1 && num_weights == 1) {
146		# recycle weights
147	9x	vi <- seq_len(num_values)
148	9x	wi <- rep.int(1, num_values)
149		}
150
151	330x	list(values = vi, weights = wi)
152		}
153
154		.areaMethod <- function(crs_obj) {
155	39x	if (!(is.na(crs_obj)) && crs_obj$units_gdal == 'degree') {
156	14x	return('spherical')
157		} else {
158	25x	return('cartesian')
159		}
160		}
161
162		.validateFlag <- function(value, name) {
163	2762x	if(!(is.logical(value) && length(value) == 1 && !is.na(value))) {
164	27x	stop(name, ' must be TRUE or FALSE')
165		}
166		}
167
168		.validateNumericScalar <- function(value, name) {
169	576x	if (!(is.numeric(value) && length(value) == 1 && !is.na(value))) {
170	4x	stop(name, ' must be a single numeric value')
171		}
172		}
173
174		.validateNumericScalarOrNA <- function(value, name) {
175	568x	if (!(is.numeric(value) && length(value) == 1)) {
176	4x	stop(name, ' must be a single numeric value')
177		}
178		}
179
180		.validateUniqueNames <- function(x) {
181	563x	nm <- names(x)
182	563x	if (!is.null(nm)) {
183	370x	if (length(nm) != length(unique(nm))) {
184	1x	stop('names of input rasters must be unique')
185		}
186		}
187		}
188
189		# faster replacement for as.data.frame when input is a named list
190		# with equal-length columns
191		# from Advanced R, sec. 24.4.2
192		.quickDf <- function(lst) {
193	162x	class(lst) <- 'data.frame'
194	162x	attr(lst, 'row.names') <- .set_row_names(length(lst[[1]]))
195	162x	lst
196		}
197
198		.singleColumnToVector <- function(df) {
199	14x	if (ncol(df) == 1) {
200	3x	df[, 1]
201		} else {
202	11x	df
203		}
204		}
205
206		# Return the number of standard (non-...) arguments in a supplied function that
207		# do not have a default value. This is used to fail if the summary function
208		# provided by the user cannot accept arguments of values and weights.
209		.num_expected_args <- function(fun) {
210	70x	a <- formals(args(fun))
211	70x	a <- a[names(a) != '...']
212	70x	sum(sapply(a, nchar) == 0)
213		}
214
215		.startReading <- function(r) {
216	521x	if(inherits(r, 'BasicRaster')) {
217	317x	return(raster::readStart(r))
218	204x	} else if (inherits(r, 'SpatRaster')) {
219	35x	terra::readStart(r)
220		}
221
222	204x	return(r)
223		}
224
225		.stopReading <- function(r) {
226	521x	if(inherits(r, 'BasicRaster')) {
227	317x	return(raster::readStop(r))
228	204x	} else if (inherits(r, 'SpatRaster')) {
229	35x	terra::readStop(r)
230		}
231
232	204x	return(r)
233		}
234
235		.crs <- function(r) {
236	12x	if(inherits(r, 'BasicRaster')) {
237	12x	if (utils::packageVersion('raster') < numeric_version('3.5')) {
238	!	return(raster::crs(r))
239		} else {
240	12x	return(terra::crs(r))
241		}
242	!	} else if (inherits(r, 'SpatRaster')) {
243	!	return(terra::crs(r))
244		} else {
245	!	stop('Unknown type: ', class(r))
246		}
247		}
248
249		.setValues <- function(r, x) {
250	13x	if(inherits(r, 'BasicRaster')) {
251	13x	if (utils::packageVersion('raster') < numeric_version('3.5')) {
252	!	raster::values(r) <- x
253		} else {
254	13x	terra::values(r) <- x
255		}
256	!	} else if (inherits(r, 'SpatRaster')) {
257	!	raster::values(r) <- x
258		} else {
259	!	stop('Unknown type: ', class(r))
260		}
261
262	13x	return(r)
263		}
264
265		.numLayers <- function(r) {
266	640x	if(inherits(r, 'BasicRaster')) {
267	520x	return(raster::nlayers(r))
268	120x	} else if (inherits(r, 'SpatRaster')) {
269	120x	return(terra::nlyr(r))
270		} else {
271	!	stop('Unknown type: ', class(r))
272		}
273		}
274
275		.isRaster <- function(r) {
276	191x	inherits(r, 'BasicRaster') \| inherits(r, 'SpatRaster')
277		}
278
279		.xFromCol <- function(r, col) {
280	30x	if (inherits(r, 'BasicRaster')) {
281	17x	raster::xFromCol(r, col)
282	13x	} else if (inherits(r, 'SpatRaster')) {
283	13x	terra::xFromCol(r, col)
284		} else {
285	!	stop('Unknown type: ', class(r))
286		}
287		}
288
289		.colFromX <- function(r, x) {
290	78x	if (inherits(r, 'BasicRaster')) {
291	30x	raster::colFromX(r, x)
292	48x	} else if (inherits(r, 'SpatRaster')) {
293	48x	terra::colFromX(r, x)
294		} else {
295	!	stop('Unknown type: ', class(r))
296		}
297		}
298
299		.yFromRow <- function(r, row) {
300	30x	if (inherits(r, 'BasicRaster')) {
301	17x	raster::yFromRow(r, row)
302	13x	} else if (inherits(r, 'SpatRaster')) {
303	13x	terra::yFromRow(r, row)
304		} else {
305	!	stop('Unknown type: ', class(r))
306		}
307		}
308
309		.rowFromY <- function(r, y) {
310	78x	if (inherits(r, 'BasicRaster')) {
311	30x	raster::rowFromY(r, y)
312	48x	} else if (inherits(r, 'SpatRaster')) {
313	48x	terra::rowFromY(r, y)
314		} else {
315	!	stop('Unknown type: ', class(r))
316		}
317		}
318
319		.cellFromRowCol <- function(r, row, col) {
320	24x	if (inherits(r, 'BasicRaster')) {
321	11x	raster::cellFromRowCol(r, row, col)
322	13x	} else if (inherits(r, 'SpatRaster')) {
323	13x	terra::cellFromRowCol(r, row, col)
324		} else {
325	!	stop('Unknown type: ', class(r))
326		}
327		}
328
329		.extent <- function(r) {
330	724x	if (inherits(r, 'BasicRaster')) {
331	572x	ex <- r@extent
332	572x	c(ex@xmin, ex@ymin, ex@xmax, ex@ymax)
333	152x	} else if (inherits(r, 'SpatRaster')) {
334	152x	ex <- terra::ext(r)
335	152x	c(ex$xmin, ex$ymin, ex$xmax, ex$ymax)
336		} else {
337	!	stop('Unknown type: ', class(r))
338		}
339		}
340
341		.res <- function(r) {
342	724x	if (inherits(r, 'BasicRaster')) {
343	572x	raster::res(r)
344	152x	} else if (inherits(r, 'SpatRaster')) {
345	152x	terra::res(r)
346		} else {
347	!	stop('Unknown type: ', class(r))
348		}
349		}
350
351		.getValuesBlock <- function(r, row, nrows, col, ncols) {
352	1054x	if (inherits(r, 'BasicRaster')) {
353	951x	raster::getValuesBlock(r, row, nrows, col, ncols, format = 'm')
354	103x	} else if (inherits(r, 'SpatRaster')) {
355	103x	terra::readValues(r, row, nrows, col, ncols, mat = TRUE)
356		} else {
357	!	stop('Unknown type: ', class(r))
358		}
359		}
360
361		.createProgress <- function(progress, n) {
362	252x	if (progress && n > 1) {
363	8x	pb <- utils::txtProgressBar(min = 0, max = n, initial=0, style=3)
364	8x	update_progress <- function() {
365	76x	i <- 1 + utils::getTxtProgressBar(pb)
366	76x	utils::setTxtProgressBar(pb, i)
367	76x	if (i == n) {
368	8x	close(pb)
369		}
370		}
371		} else {
372	244x	update_progress <- function() {}
373		}
374
375	252x	return(update_progress)
376		}
377
378		.isInMemory <- function(r) {
379	297x	if (inherits(r, 'BasicRaster')) {
380	272x	return(raster::inMemory(r))
381	25x	} else if (inherits(r, 'SpatRaster')) {
382	25x	return(terra::inMemory(r[[1]]))
383		} else {
384	!	stop('Unknown type: ', class(r))
385		}
386		}
387
388		.netCDFBlockSize <- function(fname, varname) {
389	1x	nc <- NULL
390	1x	sz <- NA
391
392	1x	tryCatch({
393	1x	nc <- ncdf4::nc_open(fname)
394	1x	sz <- nc$var[[varname]]$chunksizes
395
396	1x	dim_index <- nc$var[[varname]]$dimids + 1L
397	1x	dim_names <- sapply(dim_index, function(i) nc$dim[[i]]$name)
398
399	1x	if (all(is.na(sz))) {
400		# file is not compressed
401	!	sz <- rep.int(1, length(dim_names))
402		}
403
404	1x	names(sz) <- dim_names
405	1x	}, finally = {
406	1x	if (!is.null(nc)) {
407	1x	ncdf4::nc_close(nc)
408		}
409		})
410
411		# flip dimensions 1 and 2 so we return row/col
412	1x	return(sz[c(2, 1, seq_along(sz)[-(1:2)])])
413		}
414
415		.blockSize <- function(r) {
416		# set default return value in case file is uncompressed and has
417		# has no block size, or we simply can't figure it out
418	8x	ret <- c(1, 1)
419
420	8x	if (inherits(r, 'BasicRaster')) {
421	2x	if (r[[1]]@file@driver == 'netcdf') {
422	1x	ret <- .netCDFBlockSize(r[[1]]@file@name, attr(r[[1]]@data, 'zvar'))
423	1x	} else if (r@file@driver == 'gdal') {
424	1x	ret <- c(r@file@blockrows, r@file@blockcols)
425		}
426	6x	} else if (inherits(r, 'SpatRaster')) {
427	6x	ret <- terra::fileBlocksize(r)[1, ]
428		}
429
430	8x	unname(ret)
431		}
432
433		.eagerLoad <- function(r, geoms, max_cells_in_memory, message_on_fail) {
434	17x	if (is.null(r)) {
435	5x	return(NULL)
436		}
437
438	12x	cells_required <- .numCells(r, geoms)
439
440	12x	if (cells_required <= max_cells_in_memory) {
441	6x	box <- sf::st_bbox(geoms)
442	6x	geom_ext <- terra::ext(box[c('xmin', 'xmax', 'ymin', 'ymax')])
443
444	6x	if (!inherits(r, 'SpatRaster')) {
445		# current CRAN version of terra (1.4-22) does not preserve
446		# names on conversion (https://github.com/rspatial/terra/issues/430)
447	!	nm <- names(r)
448	!	r <- terra::rast(r)
449	!	names(r) <- nm
450		}
451
452	6x	overlap_ext <- terra::intersect(terra::ext(r), geom_ext)
453	6x	if (is.null(overlap_ext)) {
454		# Extents do not overlap, and terra::crop will throw an error
455		# if we try to crop. Return the input raster as-is; nothing will be
456		# read from it anyway.
457	!	return(r)
458		}
459
460	6x	r <- terra::crop(r, geom_ext, snap = 'out')
461	6x	} else if (message_on_fail) {
462	3x	message('Cannot preload entire working area of ', cells_required,
463	3x	' cells with max_cells_in_memory = ', max_cells_in_memory, '.',
464	3x	' Raster values will be read for each feature individually.')
465		}
466
467	12x	return(r)
468		}
469
470		.numCells <- function(r, g) {
471	12x	if (is.null(r)) {
472	!	return(0)
473		}
474
475	12x	box <- sf::st_bbox(g)
476
477	12x	top <- .rowFromY(r, box['ymax'])
478	12x	bottom <- .rowFromY(r, box['ymin'])
479	12x	left <- .colFromX(r, box['xmin'])
480	12x	right <- .colFromX(r, box['xmax'])
481
482	12x	if (is.na(top) && is.na(bottom)) {
483	!	return(0L)
484		}
485	12x	if (is.na(left) && is.na(right)) {
486	!	return(0L)
487		}
488	12x	if (is.na(top)) {
489	1x	top <- 1
490		}
491	12x	if (is.na(bottom)) {
492	!	bottom <- nrow(r)
493		}
494	12x	if (is.na(left)) {
495	!	left <- 1
496		}
497	12x	if (is.na(right)) {
498	1x	right <- ncol(r)
499		}
500
501	12x	return( (bottom - top + 1) * (right - left + 1) * .numLayers(r) )
502		}
503
504		splice <- function(x, i, replacement) {
505	12x	c(x[seq_along(x) < i],
506	12x	replacement,
507	12x	x[seq_along(x) > i])
508		}

1		# Copyright (c) 2020-2022 ISciences, LLC.
2		# All rights reserved.
3		#
4		# This software is licensed under the Apache License, Version 2.0 (the "License").
5		# You may not use this file except in compliance with the License. You may
6		# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		#
8		# Unless required by applicable law or agreed to in writing, software
9		# distributed under the License is distributed on an "AS IS" BASIS,
10		# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		# See the License for the specific language governing permissions and
12		# limitations under the License.
13
14	23x	setGeneric("exact_resample", function(x, y, ...)
15	23x	standardGeneric("exact_resample"))
16
17		#' Resample a raster to a new grid
18		#'
19		#' @param x a \code{RasterLayer} or \code{SpatRaster} to be resampled
20		#' @param y a raster of the same class as \code{x} with a grid definition to
21		#' which \code{x} should be resampled
22		#' @param fun a named summary operation or R function to be used for the resampling
23		#' @param coverage_area use cell coverage areas instead of coverage fractions
24		#' in \code{fun}
25		#' @return a resampled version of \code{x}, returned as a \code{RasterLayer} or
26		#' \code{SpatRaster}, depending on the values of \code{x} and \code{y}
27		#'
28		#' @name exact_resample
29		NULL
30
31		.exact_resample <- function(x, y, fun, coverage_area = FALSE) {
32	23x	analysis_crs <- sf::st_crs(x)
33	23x	if (sf::st_crs(x) != sf::st_crs(y)) {
34	3x	if (is.na(sf::st_crs(x))) {
35	1x	warning("No CRS specified for source raster; assuming it has the same CRS as destination raster.")
36	1x	analysis_crs <- sf::st_crs(y)
37	2x	} else if (is.na(sf::st_crs(y))) {
38	1x	warning("No CRS specified for destination raster; assuming it has the same CRS as source raster.")
39		} else {
40	1x	stop('Destination raster must have same CRS as source.')
41		}
42		}
43
44	22x	if (is.character(fun)) {
45	15x	if (length(fun) != 1) {
46	2x	stop("Only a single operation may be used for resampling.")
47		}
48
49	13x	if (startsWith(fun, 'weighted')) {
50	1x	stop("Weighted operations cannot be used for resampling.")
51		}
52
53	12x	if (.numLayers(x) != 1) {
54	1x	stop("Raster to be resampled must have a single layer.")
55		}
56
57	11x	summary_stat <- fun
58	11x	summary_fun <- NULL
59		} else {
60	7x	if (!is.function(fun)) {
61	!	stop("fun must be a named summary operation or an R function")
62		}
63
64	7x	if (.num_expected_args(fun) != 2) {
65	1x	stop("exact_extract was called with a function that does not appear to ",
66	1x	"be of the form `function(values, coverage_fractions)`.")
67		}
68
69	6x	summary_stat <- NULL
70	6x	summary_fun <- fun
71		}
72
73	17x	.validateFlag(coverage_area, 'coverage_area')
74
75	17x	area_method <- .areaMethod(analysis_crs)
76
77	17x	x <- .startReading(x)
78	17x	tryCatch({
79	17x	ret <- CPP_resample(x,
80	17x	y,
81	17x	summary_stat,
82	17x	summary_fun,
83	17x	coverage_area,
84	17x	area_method)
85	13x	if (inherits(x, 'SpatRaster')) {
86	7x	terra::rast(ret)
87		} else {
88	6x	ret
89		}
90	17x	}, finally={
91	17x	.stopReading(x)
92		})
93		}
94
95		#' @import sf
96		#' @import raster
97		#' @useDynLib exactextractr
98		#' @rdname exact_resample
99		#' @export
100		setMethod('exact_resample', signature(x='RasterLayer', y='RasterLayer'), .exact_resample)
101
102		#' @useDynLib exactextractr
103		#' @rdname exact_resample
104		#' @export
105		setMethod('exact_resample', signature(x='SpatRaster', y='SpatRaster'), .exact_resample)

1		# Copyright (c) 2018-2021 ISciences, LLC.
2		# All rights reserved.
3		#
4		# This software is licensed under the Apache License, Version 2.0 (the "License").
5		# You may not use this file except in compliance with the License. You may
6		# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		#
8		# Unless required by applicable law or agreed to in writing, software
9		# distributed under the License is distributed on an "AS IS" BASIS,
10		# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		# See the License for the specific language governing permissions and
12		# limitations under the License.
13
14	13x	setGeneric('coverage_fraction', function(x, y, crop=FALSE, ...)
15	13x	standardGeneric('coverage_fraction'))
16
17		.coverage_fraction <- function(x, y, crop) {
18	13x	if(is.na(sf::st_crs(x)) && !is.na(sf::st_crs(y))) {
19	1x	warning("No CRS specified for raster; assuming it has the same CRS as the polygons.")
20	12x	} else if(is.na(sf::st_crs(y)) && !is.na(sf::st_crs(x))) {
21	1x	warning("No CRS specified for polygons; assuming they have the same CRS as the raster.")
22	11x	} else if(sf::st_crs(x) != sf::st_crs(y)) {
23	1x	y <- sf::st_transform(y, sf::st_crs(x))
24	1x	warning("Polygons transformed to raster CRS (EPSG:", sf::st_crs(x)$epsg, ")")
25		}
26
27	13x	lapply(sf::st_as_binary(y, EWKB=TRUE), function(wkb) {
28	13x	CPP_coverage_fraction(x, wkb, crop)
29		})
30		}
31
32		#' Compute the fraction of raster cells covered by a polygon
33		#'
34		#' @param x a (possibly empty) \code{RasterLayer} whose resolution and
35		#' extent will be used for the generated \code{RasterLayer}.
36		#' @param y a \code{sf} object with polygonal geometries
37		#' @param crop if \code{TRUE}, each generated \code{RasterLayer} will be
38		#' cropped to the extent of its associated feature.
39		#' @return a list with a \code{RasterLayer} for each feature in \code{y}.
40		#' Values of the raster represent the fraction of each
41		#' cell in \code{x} that is covered by \code{y}.
42		#' @examples
43		#' rast <- raster::raster(matrix(1:100, ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
44		#' poly <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))')
45		#'
46		#' cov_frac <- coverage_fraction(rast, poly)[[1]]
47		#' @name coverage_fraction
48		NULL
49
50		#' @import sf
51		#' @import raster
52		#' @useDynLib exactextractr
53		#' @rdname coverage_fraction
54		#' @export
55		setMethod('coverage_fraction', signature(x='RasterLayer', y='sf'), function(x, y, crop=FALSE) {
56	!	coverage_fraction(x, sf::st_geometry(y), crop)
57		})
58
59		#' @rdname coverage_fraction
60		#' @export
61		setMethod('coverage_fraction', signature(x='RasterLayer', y='sfc_MULTIPOLYGON'), .coverage_fraction)
62
63		#' @rdname coverage_fraction
64		#' @export
65		setMethod('coverage_fraction', signature(x='RasterLayer', y='sfc_POLYGON'), .coverage_fraction)
66
67		#' @rdname coverage_fraction
68		#' @export
69		setMethod('coverage_fraction', signature(x='SpatRaster', y='sf'), function(x, y, crop=FALSE) {
70	!	coverage_fraction(x, sf::st_geometry(y), crop)
71		})
72
73		#' @rdname coverage_fraction
74		#' @export
75		setMethod('coverage_fraction', signature(x='SpatRaster', y='sfc_MULTIPOLYGON'), .coverage_fraction)
76
77		#' @rdname coverage_fraction
78		#' @export
79		setMethod('coverage_fraction', signature(x='SpatRaster', y='sfc_POLYGON'), .coverage_fraction)

1		// Copyright (c) 2018-2021 ISciences, LLC.
2		// All rights reserved.
3		//
4		// This software is licensed under the Apache License, Version 2.0 (the "License").
5		// You may not use this file except in compliance with the License. You may
6		// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		//
8		// Unless required by applicable law or agreed to in writing, software
9		// distributed under the License is distributed on an "AS IS" BASIS,
10		// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		// See the License for the specific language governing permissions and
12		// limitations under the License.
13
14		#include <Rcpp.h>
15
16		#include "geos_r.h"
17
18		#include "raster_utils.h"
19
20		#include "exactextract/src/raster.h"
21		#include "exactextract/src/raster_cell_intersection.h"
22
23		using exactextract::RasterView;
24		using exactextract::raster_cell_intersection;
25
26		// [[Rcpp::export]]
27	13x	Rcpp::S4 CPP_coverage_fraction(Rcpp::S4 & rast,
28		const Rcpp::RawVector & wkb,
29		bool crop)
30		{
31		try {
32	26x	GEOSAutoHandle geos;
33	39x	Rcpp::Environment xx = Rcpp::Environment::namespace_env("exactextractr");
34
35	13x	auto grid = make_grid(rast);
36	26x	auto coverage_fraction = raster_cell_intersection(grid, geos.handle, read_wkb(geos.handle, wkb).get());
37
38	13x	if (crop) {
39	2x	grid = coverage_fraction.grid();
40		}
41	26x	RasterView<float> coverage_view(coverage_fraction, grid);
42
43	13x	Rcpp::NumericMatrix weights{static_cast<int>(grid.rows()),
44	26x	static_cast<int>(grid.cols())};
45
46	904x	for (size_t i = 0; i < grid.rows(); i++) {
47	238204x	for (size_t j = 0; j < grid.cols(); j++) {
48	237313x	weights(i, j) = coverage_view(i, j);
49	237313x	if (!crop && std::isnan(weights(i, j))) {
50	209304x	weights(i, j) = 0;
51		}
52		}
53		}
54
55	13x	if (rast.inherits("SpatRaster")) {
56	3x	Rcpp::Environment terra = Rcpp::Environment::namespace_env("terra");
57	3x	Rcpp::Function rastFn = terra["rast"];
58	3x	Rcpp::Function extFn = terra["ext"];
59	3x	Rcpp::Function crsFn = terra["crs"];
60
61	1x	Rcpp::S4 ext = extFn(grid.xmin(), grid.xmax(), grid.ymin(), grid.ymax());
62
63		return rastFn(weights,
64	2x	Rcpp::Named("ext") = ext,
65	3x	Rcpp::Named("crs") = crsFn(rast));
66		} else {
67	36x	Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster");
68	36x	Rcpp::Function rasterFn = raster["raster"];
69	24x	Rcpp::Function crsFn = xx[".crs"];
70
71		return rasterFn(weights,
72	24x	Rcpp::Named("xmn")=grid.xmin(),
73	24x	Rcpp::Named("xmx")=grid.xmax(),
74	24x	Rcpp::Named("ymn")=grid.ymin(),
75	24x	Rcpp::Named("ymx")=grid.ymax(),
76	36x	Rcpp::Named("crs")=crsFn(rast));
77		}
78		} catch (std::exception & e) {
79		// throw predictable exception class
80		#ifdef __SUNPRO_CC
81		// Rcpp::stop crashes CRAN Solaris build
82		// https://github.com/RcppCore/Rcpp/issues/1159
83		Rf_error(e.what());
84		return R_NilValue;
85		#else
86		Rcpp::stop(e.what());
87		#endif
88		}
89		}

1		// Copyright (c) 2018 ISciences, LLC.
2		// All rights reserved.
3		//
4		// This software is licensed under the Apache License, Version 2.0 (the "License").
5		// You may not use this file except in compliance with the License. You may
6		// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		//
8		// Unless required by applicable law or agreed to in writing, software
9		// distributed under the License is distributed on an "AS IS" BASIS,
10		// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		// See the License for the specific language governing permissions and
12		// limitations under the License.
13
14		#ifndef EXACTEXTRACT_GRID_H
15		#define EXACTEXTRACT_GRID_H
16
17		#include <numeric>
18		#include <stdexcept>
19		#include <vector>
20
21		#include "box.h"
22
23		namespace exactextract {
24		struct infinite_extent {
25		static const size_t padding = 1;
26		};
27
28		struct bounded_extent {
29		static const size_t padding = 0;
30		};
31
32	391x	static inline bool is_integral(double d, double tol) {
33	391x	return std::abs(d - std::round(d)) <= tol;
34		}
35
36		template<typename extent_tag>
37		class Grid {
38
39		public:
40	1836196x	Grid(const Box & extent, double dx, double dy) :
41		m_extent{extent},
42		m_dx{dx},
43		m_dy{dy},
44	1836196x	m_num_rows{2*extent_tag::padding + (extent.ymax > extent.ymin ? static_cast<size_t>(std::round((extent.ymax - extent.ymin) / dy)) : 0)},
45	3672392x	m_num_cols{2*extent_tag::padding + (extent.xmax > extent.xmin ? static_cast<size_t>(std::round((extent.xmax - extent.xmin) / dx)) : 0)}
46		{}
47
48	118372x	static Grid make_empty() {
49	118372x	return Grid({0, 0, 0, 0}, 0, 0);
50		}
51
52	4794273x	size_t get_column(double x) const {
53		if (extent_tag::padding) {
54	6384784x	if (x < m_extent.xmin)
55	878x	return 0;
56	6383906x	if (x > m_extent.xmax)
57	1180x	return m_num_cols - 1;
58	6382726x	if (x == m_extent.xmax) // special case, returning the cell for which xmax is the right
59	3192320x	return m_num_cols - 2;
60		} else {
61	1601881x	if (x < m_extent.xmin \|\| x > m_extent.xmax)
62	!	throw std::out_of_range("x");
63
64	1601881x	if (x == m_extent.xmax)
65	801682x	return m_num_cols - 1;
66		}
67
68		// Since we've already range-checked our x value, make sure that
69		// the computed column index is no greater than the column index
70		// associated with xmax.
71	7186206x	return std::min(
72	2395402x	extent_tag::padding + static_cast<size_t>(std::floor((x - m_extent.xmin) / m_dx)),
73	4790804x	get_column(m_extent.xmax));
74		}
75
76	4793002x	size_t get_row(double y) const {
77		if (extent_tag::padding) {
78	6382114x	if (y > m_extent.ymax)
79	2246x	return 0;
80	6379868x	if (y < m_extent.ymin)
81	1896x	return m_num_rows - 1;
82	6377972x	if (y == m_extent.ymin) // special case, returning the cell for which ymin is the bottom
83	3190236x	return m_num_rows - 2;
84		} else {
85	1601945x	if (y < m_extent.ymin \|\| y > m_extent.ymax)
86	!	throw std::out_of_range("y");
87
88	1601945x	if (y == m_extent.ymin)
89	801682x	return m_num_rows - 1;
90		}
91
92		// Since we've already range-checked our y value, make sure that
93		// the computed row index is no greater than the column index
94		// associated with ymin.
95	7182393x	return std::min(
96	2394131x	extent_tag::padding + static_cast<size_t>(std::floor((m_extent.ymax - y) / m_dy)),
97	4788262x	get_row(m_extent.ymin));
98		}
99
100	29924212x	bool empty() const { return m_num_rows <= 2extent_tag::padding && m_num_cols <= 2extent_tag::padding; }
101
102	38787878x	size_t rows() const { return m_num_rows; }
103
104	115537109x	size_t cols() const { return m_num_cols; }
105
106	199x	size_t size() const { return rows()*cols(); }
107
108	5731163x	double xmin() const { return m_extent.xmin; }
109
110	3811117x	double xmax() const { return m_extent.xmax; }
111
112	3671541x	double ymin() const { return m_extent.ymin; }
113
114	5957169x	double ymax() const { return m_extent.ymax; }
115
116	6247212x	double dx() const { return m_dx; }
117
118	6475059x	double dy() const { return m_dy; }
119
120	2236374x	const Box& extent() const { return m_extent; }
121
122	400437x	size_t row_offset(const Grid & other) const { return static_cast<size_t>(std::round(std::abs(other.m_extent.ymax - m_extent.ymax) / m_dy)); }
123
124	400437x	size_t col_offset(const Grid & other) const { return static_cast<size_t>(std::round(std::abs(m_extent.xmin - other.m_extent.xmin) / m_dx)); }
125
126	26314x	double x_for_col(size_t col) const { return m_extent.xmin + ((col - extent_tag::padding) + 0.5) * m_dx; }
127
128	13356x	double y_for_row(size_t row) const { return m_extent.ymax - ((row - extent_tag::padding) + 0.5) * m_dy; }
129
130	1555x	Grid<extent_tag> crop(const Box & b) const {
131	1555x	if (extent().intersects(b)) {
132	1471x	return shrink_to_fit(b.intersection(extent()));
133		} else {
134	84x	return make_empty();
135		}
136		}
137
138	400841x	Grid<extent_tag> shrink_to_fit(const Box & b) const {
139	400841x	if (b.xmin < m_extent.xmin \|\| b.ymin < m_extent.ymin \|\| b.xmax > m_extent.xmax \|\| b.ymax > m_extent.ymax) {
140	!	throw std::range_error("Cannot shrink extent to bounds larger than original.");
141		}
142
143	400841x	size_t col0 = get_column(b.xmin);
144	400841x	size_t row1 = get_row(b.ymax);
145
146		// Shrink xmin and ymax to fit the upper-left corner of the supplied extent
147	400841x	double snapped_xmin = m_extent.xmin + (col0 - extent_tag::padding) * m_dx;
148	400841x	double snapped_ymax = m_extent.ymax - (row1 - extent_tag::padding) * m_dy;
149
150		// Make sure x0 and y1 are within the reduced extent. Because of
151		// floating point round-off errors, this is not always the case.
152	400841x	if (b.xmin < snapped_xmin) {
153	!	snapped_xmin -= m_dx;
154	!	col0--;
155		}
156
157	400841x	if (b.ymax > snapped_ymax) {
158	!	snapped_ymax += m_dy;
159	!	row1--;
160		}
161
162	400841x	size_t col1 = get_column(b.xmax);
163	400841x	size_t row0 = get_row(b.ymin);
164
165	400841x	size_t num_rows = 1 + (row0 - row1);
166	400841x	size_t num_cols = 1 + (col1 - col0);
167
168		// If xmax or ymin falls cleanly on a cell boundary, we don't
169		// need as many rows or columns as we otherwise would, because
170		// we assume that the rightmost cell of the grid is a closed
171		// interval in x, and the lowermost cell of the grid is a
172		// closed interval in y.
173	400841x	if (num_rows > 2 && (snapped_ymax - (num_rows-1)*m_dy <= b.ymin)) {
174	465x	num_rows--;
175		}
176	400841x	if (num_cols > 2 && (snapped_xmin + (num_cols-1)*m_dx >= b.xmax)) {
177	453x	num_cols--;
178		}
179
180		// Perform offsets relative to the new xmin, ymax origin
181		// points. If this is not done, then floating point roundoff
182		// error can cause progressive shrink() calls with the same
183		// inputs to produce different results.
184	400841x	Box reduced_box = {
185		snapped_xmin,
186	400841x	std::min(snapped_ymax - num_rows * m_dy, b.ymin),
187	400841x	std::max(snapped_xmin + num_cols * m_dx, b.xmax),
188		snapped_ymax
189		};
190
191		// Fudge computed xmax and ymin, if needed, to prevent extent
192		// from growing during a shrink operation.
193	400841x	if (reduced_box.xmax > m_extent.xmax) {
194	!	if (std::round((reduced_box.xmax - reduced_box.xmin)/m_dx) ==
195	!	std::round((m_extent.xmax - reduced_box.xmin)/m_dx)) {
196	!	reduced_box.xmax = m_extent.xmax;
197		} else {
198	!	throw std::runtime_error("Shrink operation failed.");
199		}
200		}
201	400841x	if (reduced_box.ymin < m_extent.ymin) {
202	3x	if (std::round((reduced_box.ymax - reduced_box.ymin)/m_dy) ==
203	3x	std::round((reduced_box.ymax - m_extent.ymin)/m_dy)) {
204	3x	reduced_box.ymin = m_extent.ymin;
205		} else {
206	!	throw std::runtime_error("Shrink operation failed.");
207		}
208		}
209
210	400841x	Grid<extent_tag> reduced{reduced_box, m_dx, m_dy};
211
212	!	if (b.xmin < reduced.xmin() \|\| b.ymin < reduced.ymin() \|\| b.xmax > reduced.xmax() \|\| b.ymax > reduced.ymax()) {
213	!	throw std::runtime_error("Shrink operation failed.");
214		}
215
216	801682x	return reduced;
217		}
218
219		template<typename extent_tag2>
220	98x	bool compatible_with(const Grid<extent_tag2> &b, double compatability_tol) const {
221		// Define a tolerance for grid compatibility, to be used in the following ways:
222		// Grid resolutions must be integer multiples of each other within a factor of 'compatibility_tol'
223		// Grid origin points must differ by an integer multiple of the smaller of the two grid resolutions,
224		// within a factor of 'compatibility_tol'.
225
226		if (empty() \|\| b.empty()) {
227	!	return true;
228		}
229
230	98x	double xtol = std::min(m_dx, b.m_dx) * compatability_tol;
231	98x	double ytol = std::min(m_dy, b.m_dy) * compatability_tol;
232
233		// Check x-resolution compatibility
234	98x	if (!is_integral(std::max(m_dx, b.m_dx) / std::min(m_dx, b.m_dx), xtol)) {
235	!	return false;
236		}
237
238		// Check y-resolution compatibility
239	98x	if (!is_integral(std::max(m_dy, b.m_dy) / std::min(m_dy, b.m_dy), ytol)) {
240	!	return false;
241		}
242
243		// Check left-hand boundary compatibility
244	98x	if (!is_integral(std::abs(b.m_extent.xmin - m_extent.xmin) / std::min(m_dx, b.m_dx), xtol)) {
245	1x	return false;
246		}
247
248		// Check upper boundary compatibility
249	97x	if (!is_integral(std::abs(b.m_extent.ymax - m_extent.ymax) / std::min(m_dy, b.m_dy), ytol)) {
250	!	return false;
251		}
252
253	97x	return true;
254		}
255
256		template<typename extent_tag2>
257	98x	Grid<extent_tag> common_grid(const Grid<extent_tag2> &b, double tol = 1e-6) const {
258	98x	if (!compatible_with(b, tol)) {
259	1x	throw std::runtime_error("Incompatible extents.");
260		}
261
262	97x	if (b.empty()) {
263	!	return *this;
264		}
265
266	97x	const double common_dx = std::min(m_dx, b.m_dx);
267	97x	const double common_dy = std::min(m_dy, b.m_dy);
268
269	97x	const double common_xmin = std::min(m_extent.xmin, b.m_extent.xmin);
270	97x	const double common_ymax = std::max(m_extent.ymax, b.m_extent.ymax);
271
272	97x	double common_xmax = std::max(m_extent.xmax, b.m_extent.xmax);
273	97x	double common_ymin = std::min(m_extent.ymin, b.m_extent.ymin);
274
275	97x	const long nx = static_cast<long>(std::round((common_xmax - common_xmin) / common_dx));
276	97x	const long ny = static_cast<long>(std::round((common_ymax - common_ymin) / common_dy));
277
278	97x	common_xmax = std::max(common_xmax, common_xmin + nx*common_dx);
279	97x	common_ymin = std::min(common_ymin, common_ymax - ny*common_dy);
280
281	97x	return {{ common_xmin, common_ymin, common_xmax, common_ymax}, common_dx, common_dy };
282		}
283
284		template<typename extent_tag2>
285		Grid<extent_tag> overlapping_grid(const Grid<extent_tag2> &b, double tol = 1e-6) const {
286		if (!compatible_with(b, tol)) {
287		throw std::runtime_error("Incompatible extents.");
288		}
289
290		if (empty() \|\| b.empty()) {
291		return make_empty();
292		}
293
294		const double common_dx = std::min(m_dx, b.m_dx);
295		const double common_dy = std::min(m_dy, b.m_dy);
296
297		const double common_xmin = std::max(m_extent.xmin, b.m_extent.xmin);
298		const double common_ymax = std::min(m_extent.ymax, b.m_extent.ymax);
299
300		double common_xmax = std::min(m_extent.xmax, b.m_extent.xmax);
301		double common_ymin = std::max(m_extent.ymin, b.m_extent.ymin);
302
303		const long nx = static_cast<long>(std::round((common_xmax - common_xmin) / common_dx));
304		const long ny = static_cast<long>(std::round((common_ymax - common_ymin) / common_dy));
305
306		common_xmax = std::max(common_xmax, common_xmin + nx*common_dx);
307		common_ymin = std::min(common_ymin, common_ymax - ny*common_dy);
308
309		return {{ common_xmin, common_ymin, common_xmax, common_ymax}, common_dx, common_dy };
310		}
311
312	399262x	bool operator==(const Grid<extent_tag> &b) const {
313		return
314	399262x	m_extent == b.m_extent &&
315	399509x	m_dx == b.m_dx &&
316	399509x	m_dy == b.m_dy;
317		}
318
319	399262x	bool operator!=(const Grid<extent_tag> &b) const {
320	399262x	return !(*this == b);
321		}
322
323		private:
324		Box m_extent;
325
326		double m_dx;
327		double m_dy;
328
329		size_t m_num_rows;
330		size_t m_num_cols;
331		};
332
333		Box grid_cell(const Grid<bounded_extent> & grid, size_t row, size_t col);
334		Box grid_cell(const Grid<infinite_extent> & grid, size_t row, size_t col);
335
336		Grid<infinite_extent> make_infinite(const Grid<bounded_extent> & grid);
337		Grid<bounded_extent> make_finite(const Grid<infinite_extent> & grid);
338
339		std::vector<Grid<bounded_extent>> subdivide(const Grid<bounded_extent> & grid, size_t max_size);
340
341		template<typename T>
342		Grid<bounded_extent> common_grid(T begin, T end) {
343		if (begin == end) {
344		return Grid<bounded_extent>::make_empty();
345		} else if (std::next(begin) == end) {
346		return begin->grid();
347		}
348		return std::accumulate(
349		std::next(begin),
350		end,
351		begin->grid(),
352		[](auto& acc, auto& op) {
353		return acc.common_grid(op.grid());
354		});
355		}
356
357		}
358
359		#endif //EXACTEXTRACT_INFINITEGRID_H

1		// Copyright (c) 2018-2020 ISciences, LLC.
2		// All rights reserved.
3		//
4		// This software is licensed under the Apache License, Version 2.0 (the "License").
5		// You may not use this file except in compliance with the License. You may
6		// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		//
8		// Unless required by applicable law or agreed to in writing, software
9		// distributed under the License is distributed on an "AS IS" BASIS,
10		// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		// See the License for the specific language governing permissions and
12		// limitations under the License.
13
14		#ifndef EXACTEXTRACT_MATRIX_H
15		#define EXACTEXTRACT_MATRIX_H
16
17		#include <iomanip>
18		#include <iterator>
19		#include <memory>
20		#include <cstring>
21		#include <vector>
22
23		namespace exactextract {
24
25		template<typename T>
26		class Matrix {
27
28		public:
29		using value_type = T;
30
31	1832812x	Matrix(size_t rows, size_t cols) :
32		m_rows{rows},
33	1832812x	m_cols{cols}
34		{
35	1832812x	if (m_rows > 0 && m_cols > 0) {
36		// new T[]() initializes to zero
37		m_data = std::unique_ptr<T[]>(new T[m_rows * m_cols]());
38		}
39		}
40
41	341x	Matrix(size_t rows, size_t cols, T value) :
42		m_rows{rows},
43	341x	m_cols{cols}
44		{
45	341x	if (m_rows > 0 && m_cols > 0) {
46		// new T[] does not initialize
47		m_data = std::unique_ptr<T[]>(new T[m_rows * m_cols]);
48		}
49
50	341x	std::fill(m_data.get(), m_data.get() + m_rows*m_cols, value);
51		}
52
53		explicit Matrix(const std::vector<std::vector<T>> & data) :
54		m_rows{data.size()},
55		m_cols{data[0].size()}
56		{
57		m_data = std::unique_ptr<T[]>(new T[m_rows*m_cols]());
58
59		auto lastpos = m_data.get();
60		for (auto& row : data) {
61		lastpos = std::copy(row.begin(), row.end(), lastpos);
62		}
63		}
64
65	517024x	Matrix(Matrix<T>&& m) noexcept :
66		m_rows{m.rows()},
67	517024x	m_cols{m.cols()}
68		{
69	517024x	m_data = std::move(m.m_data);
70		}
71
72	237684786x	T& operator()(size_t row, size_t col) {
73	237684786x	check(row, col);
74	237684786x	return m_data[row*m_cols + col];
75		}
76
77	54194793x	T operator()(size_t row, size_t col) const {
78	54194793x	check(row, col);
79	54194793x	return m_data[row*m_cols + col];
80		}
81
82		bool operator==(const Matrix<T> & other) const {
83		if (m_rows != other.m_rows) {
84		return false;
85		}
86		if (m_cols != other.m_cols) {
87		return false;
88		}
89
90		return 0 == memcmp(m_data.get(), other.m_data.get(), m_rowsm_colssizeof(T));
91		}
92
93	26984725x	void increment(size_t row, size_t col, const T & val) {
94	26984725x	check(row, col);
95	26984725x	m_data[row*m_cols + col] += val;
96		}
97
98	10022068x	size_t rows() const { return m_rows; }
99	115488381x	size_t cols() const { return m_cols; }
100
101		T* row(size_t row) {
102		return &(m_data[row*m_cols]);
103		}
104
105		T* data() {
106		return m_data.get();
107		}
108
109		#ifdef MATRIX_CHECK_BOUNDS
110		void check(size_t row, size_t col) const {
111		if (row + 1 > m_rows) {
112		throw std::out_of_range("Row " + std::to_string(row) + " is out of range.");
113		}
114		if (col + 1 > m_cols) {
115		throw std::out_of_range("Col " + std::to_string(col) + " is out of range.");
116		}
117		}
118		#else
119	318864304x	void check(size_t, size_t) const {}
120		#endif
121
122		private:
123		std::unique_ptr<T[]> m_data;
124
125		size_t m_rows;
126		size_t m_cols;
127
128		};
129
130		template<typename T>
131		std::ostream& operator<<(std::ostream & os, const Matrix<T> & m) {
132		for (size_t i = 0; i < m.rows(); i++) {
133		for (size_t j = 0; j < m.cols(); j++) {
134		if (m(i, j) != 0) {
135		os << std::right << std::fixed << std::setw(10) << std::setprecision(6) <<
136		m(i, j) << " ";
137		} else {
138		os << " ";
139		}
140		}
141		os << std::endl;
142		}
143
144		return os;
145		}
146
147		}
148
149		#endif

1		// Copyright (c) 2018-2022 ISciences, LLC.
2		// All rights reserved.
3		//
4		// This software is licensed under the Apache License, Version 2.0 (the "License").
5		// You may not use this file except in compliance with the License. You may
6		// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		//
8		// Unless required by applicable law or agreed to in writing, software
9		// distributed under the License is distributed on an "AS IS" BASIS,
10		// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		// See the License for the specific language governing permissions and
12		// limitations under the License.
13
14		// [[Rcpp::plugins("cpp14")]]
15		#include <Rcpp.h>
16
17		#include <memory>
18
19		#include "geos_r.h"
20		#include "raster_utils.h"
21
22		#include "s4_raster_source.h"
23
24		#include "exactextract/src/geos_utils.h"
25		#include "exactextract/src/grid.h"
26		#include "exactextract/src/matrix.h"
27		#include "exactextract/src/raster_cell_intersection.h"
28		#include "exactextract/src/raster_source.h"
29		#include "exactextract/src/raster_stats.h"
30
31		using exactextract::Box;
32		using exactextract::Matrix;
33		using exactextract::Grid;
34		using exactextract::subdivide;
35		using exactextract::bounded_extent;
36		using exactextract::Raster;
37		using exactextract::RasterView;
38		using exactextract::raster_cell_intersection;
39		using exactextract::RasterStats;
40		using exactextract::RasterSource;
41
42		#include <vector>
43		#include <string>
44
45		// [[Rcpp::export]]
46	153x	Rcpp::List CPP_exact_extract(Rcpp::S4 & rast,
47		Rcpp::Nullable<Rcpp::S4> & rast_uncropped,
48		Rcpp::Nullable<Rcpp::S4> & weights,
49		const Rcpp::RawVector & wkb,
50		double default_value,
51		double default_weight,
52		bool include_xy,
53		bool include_cell_number,
54		bool include_area,
55		bool area_weights,
56		bool coverage_areas,
57		Rcpp::Nullable<Rcpp::CharacterVector> & p_area_method,
58		Rcpp::Nullable<Rcpp::List> & include_cols,
59		Rcpp::CharacterVector & src_names,
60		Rcpp::Nullable<Rcpp::CharacterVector> & p_weights_names,
61		bool warn_on_disaggregate,
62		double grid_compat_tol) {
63		try {
64	306x	GEOSAutoHandle geos;
65	459x	Rcpp::Function names("names");
66
67	153x	auto grid = make_grid(rast);
68	153x	auto weights_grid = exactextract::Grid<bounded_extent>::make_empty();
69	153x	auto common_grid = grid;
70
71	306x	S4RasterSource rsrc(rast, default_value);
72	153x	int src_nlayers = get_nlayers(rast);
73
74	153x	std::unique_ptr<S4RasterSource> rweights;
75	153x	int weights_nlayers = 0;
76	306x	Rcpp::CharacterVector weights_names;
77
78	306x	std::string area_method;
79	153x	if (p_area_method.isNotNull()) {
80	11x	Rcpp::CharacterVector amethod = p_area_method.get();
81	11x	area_method = amethod[0];
82		}
83
84	153x	if (weights.isNotNull()) {
85	86x	Rcpp::S4 weights_s4 = weights.get();
86	43x	weights_nlayers = get_nlayers(weights_s4);
87	43x	weights_grid = make_grid(weights_s4);
88
89	43x	common_grid = grid.common_grid(weights_grid, grid_compat_tol);
90
91	43x	rweights = std::make_unique<S4RasterSource>(weights_s4, default_weight);
92	43x	weights_names = p_weights_names.get();
93
94	43x	if (warn_on_disaggregate && (common_grid.dx() < grid.dx() \|\| common_grid.dy() < grid.dy())) {
95	3x	Rcpp::warning("value raster implicitly disaggregated to match higher resolution of weights");
96		}
97	110x	} else if (area_weights) {
98	1x	weights_nlayers = 1;
99	1x	weights_names = p_weights_names.get();
100		}
101
102	306x	auto geom = read_wkb(geos.handle, wkb);
103
104	153x	auto bbox = exactextract::geos_get_box(geos.handle, geom.get());
105
106	153x	common_grid = common_grid.crop(bbox);
107
108	306x	auto coverage_fractions = raster_cell_intersection(common_grid, geos.handle, geom.get());
109	153x	auto& cov_grid = coverage_fractions.grid();
110
111		// We can't construct an Rcpp::DataFrame with a variable number of columns
112		// because of a bug in Rcpp (see https://github.com/RcppCore/Rcpp/pull/1099)
113		// Instead, we build up a list, and then create a data frame from the list.
114		// Profiling shows that this ends up in as.data.frame, which is slow.
115		// Once Rcpp 1.0.6 is released, this can be reworked to be simpler and faster.
116	306x	Rcpp::List cols;
117
118	306x	Rcpp::NumericVector coverage_vec = as_vector(coverage_fractions);
119	153x	if (coverage_areas) {
120	8x	auto areas = get_area_raster(area_method, cov_grid);
121	4x	Rcpp::NumericVector area_vec = as_vector(*areas);
122	4x	coverage_vec = coverage_vec * area_vec;
123		}
124	306x	Rcpp::LogicalVector covered = coverage_vec > 0;
125
126	153x	if (include_cols.isNotNull()) {
127	10x	Rcpp::List include_cols_list = include_cols.get();
128	10x	Rcpp::CharacterVector include_names = include_cols_list.attr("names");
129
130	10x	for (int i = 0; i < include_names.size(); i++) {
131	5x	std::string name(include_names[i]);
132	5x	cols[name] = include_cols_list[name];
133		}
134		}
135
136	341x	for (int i = 0; i < src_nlayers; i++) {
137	376x	auto values = rsrc.read_box(cov_grid.extent(), i);
138	188x	const NumericVectorRaster* r = static_cast<NumericVectorRaster*>(values.get());
139
140		// TODO Perhaps extend this to preserve types (integer, logical.)
141		// A bit challenging, since we don't know the type of the raster
142		// until we call getValuesBlock, and even then we can't be sure
143		// that the returned type is correct (as in a RasterStack with mixed types.)
144		// Since R integers are only 32-bit, we are not going to lose data by
145		// converting everything to numeric, although we pay a storage penalty.
146	188x	Rcpp::NumericVector value_vec = r->vec();
147	356x	if (grid.dx() != cov_grid.dx() \|\| grid.dy() != cov_grid.dy() \|\|
148	168x	value_vec.size() != covered.size()) {
149		// Transform values to same grid as coverage fractions
150	40x	RasterView<double> rt(*r, cov_grid);
151	40x	value_vec = as_vector(rt);
152		}
153
154	188x	value_vec = value_vec[covered];
155	188x	cols[std::string(src_names[i])] = value_vec;
156		}
157
158	212x	for (int i = 0; i < weights_nlayers; i++) {
159	118x	Rcpp::NumericVector weight_vec;
160	59x	if (area_weights) {
161	1x	auto weights = get_area_raster(area_method, cov_grid);
162	1x	weight_vec = as_vector(*weights);
163		} else {
164	116x	auto values = rweights->read_box(cov_grid.extent(), i);
165	58x	const NumericVectorRaster* r = static_cast<NumericVectorRaster*>(values.get());
166	58x	weight_vec = r->vec();
167
168	96x	if (weights_grid.dx() != cov_grid.dx() \|\| weights_grid.dy() != cov_grid.dy() \|\|
169	38x	weight_vec.size() != covered.size()) {
170		// Transform weights to same grid as coverage fractions
171	22x	RasterView<double> rt (*r, cov_grid);
172
173	22x	weight_vec = as_vector(rt);
174		}
175		}
176
177	59x	weight_vec = weight_vec[covered];
178
179	59x	std::string colname(weights_names[i]);
180	59x	if (cols.containsElementNamed(colname.c_str())) {
181		// append ".1" to the column name, to match the behavior of
182		// data.frame(). We're safe to just add ".1" (instead of incrementing
183		// .1 to .2, for example) because duplicated names within the values
184		// or weight stack will already have been made unique by raster::stack()
185	2x	colname = colname + ".1"; // append .1
186		}
187	59x	cols[colname] = weight_vec;
188		}
189
190	153x	if (include_xy) {
191		// Include xy values from whichever input raster has the higher resolution
192	30x	if (weights_nlayers > 0 && (weights_grid.dx() < grid.dx() \|\| weights_grid.dy() < grid.dy())) {
193	1x	Rcpp::S4 weights_s4 = weights.get();
194	1x	cols["x"] = get_x_values(weights_s4, cov_grid)[covered];
195	1x	cols["y"] = get_y_values(weights_s4, cov_grid)[covered];
196		} else {
197	29x	cols["x"] = get_x_values(rast, cov_grid)[covered];
198	29x	cols["y"] = get_y_values(rast, cov_grid)[covered];
199		}
200		}
201
202	153x	if (include_cell_number) {
203	24x	if (rast_uncropped.isNull()) {
204	12x	cols["cell"] = get_cell_numbers(rast, cov_grid)[covered];
205		} else {
206	12x	Rcpp::S4 tmp = rast_uncropped.get();
207	12x	cols["cell"] = get_cell_numbers(tmp, cov_grid)[covered];
208		}
209		}
210
211	153x	if (include_area) {
212	12x	auto area_rast = get_area_raster(area_method, cov_grid);
213	6x	Rcpp::NumericVector area_vec = as_vector(*area_rast);
214	6x	cols["area"] = area_vec[covered];
215		}
216
217	153x	if (coverage_areas) {
218	4x	cols["coverage_area"] = coverage_vec[covered];
219		} else {
220	149x	cols["coverage_fraction"] = coverage_vec[covered];
221		}
222
223	306x	return cols;
224		} catch (std::exception & e) {
225		// throw predictable exception class
226		#ifdef __SUNPRO_CC
227		// Rcpp::stop crashes CRAN Solaris build
228		// https://github.com/RcppCore/Rcpp/issues/1159
229		Rf_error(e.what());
230		return R_NilValue;
231		#else
232		Rcpp::stop(e.what());
233		#endif
234		}
235		}
236
237		enum class WeightingMethod {
238		NONE,
239		RASTER,
240		AREA
241		};
242
243
244	208x	static int get_num_stats(const Rcpp::StringVector & stats,
245		const std::size_t num_quantiles,
246		const std::size_t num_unique_values) {
247	208x	int num_stats = 0;
248	452x	for (const auto & stat : stats) {
249	724x	if (stat == std::string("frac") \|\|
250	480x	stat == std::string("weighted_frac")) {
251	10x	num_stats += static_cast<int>(num_unique_values);
252	234x	} else if (stat == std::string("quantile")) {
253	12x	num_stats += static_cast<int>(num_quantiles);
254		} else {
255	222x	num_stats += 1;
256		}
257		}
258
259	208x	return num_stats;
260		}
261
262		// Return a matrix with one row per stat and one row per raster layer
263		// [[Rcpp::export]]
264	214x	Rcpp::NumericMatrix CPP_stats(Rcpp::S4 & rast,
265		Rcpp::Nullable<Rcpp::S4> weights,
266		const Rcpp::RawVector & wkb,
267		double default_value,
268		double default_weight,
269		bool coverage_areas,
270		Rcpp::Nullable<Rcpp::CharacterVector> & p_area_method,
271		const Rcpp::StringVector & stats,
272		int max_cells_in_memory,
273		double grid_compat_tol,
274		const Rcpp::Nullable<Rcpp::NumericVector> & quantiles) {
275		try {
276	428x	GEOSAutoHandle geos;
277
278	214x	if (max_cells_in_memory < 1) {
279	1x	Rcpp::stop("Invalid value for max_cells_in_memory: %d", max_cells_in_memory);
280		}
281
282	213x	int nlayers = get_nlayers(rast);
283
284	426x	S4RasterSource rsrc(rast, default_value);
285
286	213x	std::unique_ptr<S4RasterSource> rweights;
287	426x	std::string area_method;
288	213x	if (p_area_method.isNotNull()) {
289	10x	area_method = ((Rcpp::CharacterVector) p_area_method.get())[0];
290		}
291
292	213x	WeightingMethod weighting = WeightingMethod::NONE;
293	213x	int nweights = 0;
294	213x	if (weights.isNotNull()) {
295	55x	weighting = WeightingMethod::RASTER;
296	55x	Rcpp::S4 weights_s4 = weights.get();
297	55x	nweights = get_nlayers(weights_s4);
298
299	55x	if (nlayers > 1 && nweights > 1 && nlayers != nweights) {
300	!	Rcpp::stop("Incompatible number of layers in value and weighting rasters");
301		}
302
303	55x	rweights = std::make_unique<S4RasterSource>(weights_s4, default_weight);
304	158x	} else if (p_area_method.isNotNull()) {
305	9x	weighting = WeightingMethod::AREA;
306	9x	nweights = 1;
307		}
308
309	426x	auto geom = read_wkb(geos.handle, wkb);
310	213x	auto bbox = exactextract::geos_get_box(geos.handle, geom.get());
311
312		auto grid = weighting == WeightingMethod::RASTER ?
313	213x	rsrc.grid().common_grid(rweights->grid(), grid_compat_tol) : rsrc.grid();
314
315	212x	bool disaggregated = (grid.dx() < rsrc.grid().dx() \|\| grid.dy() < rsrc.grid().dy());
316
317	212x	bool store_values = false;
318	212x	bool calc_value_set = false;
319	212x	std::size_t num_unique_values = 0;
320	212x	std::size_t num_quantiles = 0;
321
322	458x	for (const auto & stat : stats) {
323	250x	store_values = store_values \|\| requires_stored_values(stat);
324
325	508x	if (disaggregated && (stat == std::string("count") \|\|
326	258x	stat == std::string("sum"))) {
327	2x	Rcpp::stop("Cannot compute 'count' or 'sum' when value raster is disaggregated to resolution of weights.");
328		}
329
330	736x	if (stat == std::string("frac") \|\|
331	488x	stat == std::string("weighted_frac")) {
332	10x	calc_value_set = true;
333		}
334
335	248x	if (stat == std::string("quantile")) {
336	14x	if (quantiles.isNotNull()) {
337	13x	Rcpp::NumericVector qvec = quantiles.get();
338	13x	num_quantiles = qvec.size();
339		}
340
341	14x	if (num_quantiles == 0) {
342	2x	Rcpp::stop("Quantiles not specified.");
343		}
344		}
345		}
346
347	208x	int nresults = std::max(nlayers, nweights);
348
349	416x	std::vector<RasterStats<double>> raster_stats;
350	208x	raster_stats.reserve(nresults);
351	431x	for (int i = 0; i < nresults; i++) {
352	223x	raster_stats.emplace_back(store_values);
353		}
354
355	208x	if (bbox.intersects(grid.extent())) {
356	199x	auto cropped_grid = grid.crop(bbox);
357
358	437x	for (const auto &subgrid : subdivide(cropped_grid, max_cells_in_memory)) {
359	476x	auto coverage_fraction = raster_cell_intersection(subgrid, geos.handle, geom.get());
360	238x	if (coverage_areas) {
361	4x	auto areas = get_area_raster(area_method, subgrid);
362	6x	for (size_t i = 0; i < coverage_fraction.rows(); i++) {
363	12x	for (size_t j = 0; j < coverage_fraction.cols(); j++) {
364	8x	coverage_fraction(i, j) = coverage_fraction(i, j) * (*areas)(i, j);
365		}
366		}
367		}
368
369	238x	auto& cov_grid = coverage_fraction.grid();
370
371	238x	if (!cov_grid.empty()) {
372	237x	if (weighting == WeightingMethod::NONE) {
373	360x	for (int i = 0; i < nlayers; i++) {
374	368x	auto values = rsrc.read_box(cov_grid.extent(), i);
375	184x	raster_stats[i].process(coverage_fraction, *values);
376		}
377	61x	} else if (weighting == WeightingMethod::AREA) {
378	18x	auto weights = get_area_raster(area_method, cov_grid);
379
380	18x	for (int i = 0; i < nlayers; i++) {
381	18x	auto values = rsrc.read_box(cov_grid.extent(), i);
382	9x	raster_stats[i].process(coverage_fraction, values, weights);
383		}
384		} else {
385	52x	if (nlayers > nweights) {
386		// recycle weights
387	4x	auto weights = rweights->read_box(cov_grid.extent(), 0);
388
389	7x	for (int i = 0; i < nlayers; i++) {
390	10x	auto values = rsrc.read_box(cov_grid.extent(), i);
391	5x	raster_stats[i].process(coverage_fraction, values, weights);
392		}
393	50x	} else if (nweights > nlayers) {
394		// recycle values
395	2x	auto values = rsrc.read_box(cov_grid.extent(), 0);
396
397	4x	for (int i = 0; i < nweights; i++) {
398	6x	auto weights = rweights->read_box(cov_grid.extent(), i);
399	3x	raster_stats[i].process(coverage_fraction, values, weights);
400		}
401		} else {
402		// process values and weights in parallel
403	100x	for (int i = 0; i < nlayers; i++) {
404	102x	auto values = rsrc.read_box(cov_grid.extent(), i);
405	102x	auto weights = rweights->read_box(cov_grid.extent(), i);
406
407	51x	raster_stats[i].process(coverage_fraction, values, weights);
408		}
409		}
410		}
411		}
412		}
413		}
414
415
416	416x	std::set<double> value_set;
417	208x	if (calc_value_set) {
418	22x	for (const auto& rs : raster_stats) {
419	36x	for (const auto& value : rs) {
420	24x	value_set.insert(value);
421		}
422		}
423	10x	num_unique_values = value_set.size();
424		}
425
426	208x	auto stat_result_size = get_num_stats(stats, num_quantiles, num_unique_values);
427	416x	Rcpp::NumericMatrix stat_results = Rcpp::no_init(nresults, stat_result_size);
428
429	208x	if (calc_value_set) {
430	10x	Rcpp::NumericVector unique_values(value_set.begin(), value_set.end());
431	10x	stat_results.attr("unique_values") = unique_values;
432		}
433
434		// rows (j) represent layers
435		// cols (i) represent stats
436	427x	for (int j = 0; j < nresults; j++) {
437	223x	const auto& rs = raster_stats[j];
438
439	223x	int i = 0;
440	490x	for(const auto& stat : stats) {
441	271x	if (stat == std::string("mean")) stat_results(j, i++) = rs.mean();
442
443	233x	else if (stat == std::string("sum")) stat_results(j, i++) = rs.sum();
444	172x	else if (stat == std::string("count")) stat_results(j, i++) = rs.count();
445
446	155x	else if (stat == std::string("min")) stat_results(j, i++) = rs.min().value_or(NA_REAL);
447	149x	else if (stat == std::string("max")) stat_results(j, i++) = rs.max().value_or(NA_REAL);
448
449	143x	else if (stat == std::string("median")) stat_results(j, i++) = rs.quantile(0.5).value_or(NA_REAL);
450
451	142x	else if (stat == std::string("mode")) stat_results(j, i++) = rs.mode().value_or(NA_REAL);
452	107x	else if (stat == std::string("majority")) stat_results(j, i++) = rs.mode().value_or(NA_REAL);
453	104x	else if (stat == std::string("minority")) stat_results(j, i++) = rs.minority().value_or(NA_REAL);
454
455	100x	else if (stat == std::string("variety")) stat_results(j, i++) = rs.variety();
456	88x	else if (stat == std::string("weighted_mean")) stat_results(j, i++) = rs.weighted_mean();
457	41x	else if (stat == std::string("weighted_sum")) stat_results(j, i++) = rs.weighted_sum();
458
459	31x	else if (stat == std::string("variance")) stat_results(j, i++) = rs.variance();
460	30x	else if (stat == std::string("stdev")) stat_results(j, i++) = rs.stdev();
461	29x	else if (stat == std::string("coefficient_of_variation")) stat_results(j, i++) = rs.coefficient_of_variation();
462
463	28x	else if (stat == std::string("quantile")) {
464	28x	Rcpp::NumericVector qvec = quantiles.get();
465	39x	for (double q : qvec) {
466	27x	stat_results(j, i++) = rs.quantile(q).value_or(NA_REAL);
467		}
468		}
469
470	14x	else if (stat == std::string("frac")) {
471	34x	for (double v : value_set) {
472	24x	stat_results(j, i++) = rs.frac(v).value_or(0);
473		}
474		}
475
476	4x	else if (stat == std::string("weighted_frac")) {
477	6x	for (double v : value_set) {
478	4x	stat_results(j, i++) = rs.weighted_frac(v).value_or(0);
479		}
480		}
481
482	8x	else Rcpp::stop("Unknown stat: " + stat);
483		}
484		}
485
486		// FIXME set dimnames here
487
488	408x	return stat_results;
489	20x	} catch (std::exception & e) {
490		// throw predictable exception class
491		#ifdef __SUNPRO_CC
492		// Rcpp::stop crashes CRAN Solaris build
493		// https://github.com/RcppCore/Rcpp/issues/1159
494		Rf_error(e.what());
495		return R_NilValue;
496		#else
497	10x	Rcpp::stop(e.what());
498		#endif
499		}
500		}

1		// Copyright (c) 2018-2019 ISciences, LLC.
2		// All rights reserved.
3		//
4		// This software is licensed under the Apache License, Version 2.0 (the "License").
5		// You may not use this file except in compliance with the License. You may
6		// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		//
8		// Unless required by applicable law or agreed to in writing, software
9		// distributed under the License is distributed on an "AS IS" BASIS,
10		// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		// See the License for the specific language governing permissions and
12		// limitations under the License.
13
14		#ifndef EXACTEXTRACT_BOX_H
15		#define EXACTEXTRACT_BOX_H
16
17		#include "coordinate.h"
18		#include "crossing.h"
19		#include "side.h"
20
21		#include <limits>
22
23		namespace exactextract {
24
25		struct Box {
26		double xmin;
27		double ymin;
28		double xmax;
29		double ymax;
30
31	8934600x	Box(double p_xmin, double p_ymin, double p_xmax, double p_ymax) :
32		xmin{p_xmin},
33		ymin{p_ymin},
34		xmax{p_xmax},
35	8934600x	ymax{p_ymax} {}
36
37		static Box maximum_finite() {
38		return {
39		std::numeric_limits<double>::lowest(),
40		std::numeric_limits<double>::lowest(),
41		std::numeric_limits<double>::max(),
42		std::numeric_limits<double>::max()
43		};
44		}
45
46	404x	static Box make_empty() {
47	404x	return {0, 0, 0, 0};
48		}
49
50	6973418x	double width() const {
51	6973418x	return xmax - xmin;
52		}
53
54	6973418x	double height() const {
55	6973418x	return ymax - ymin;
56		}
57
58	6568936x	double area() const {
59	6568936x	return width() * height();
60		}
61
62	202241x	double perimeter() const {
63	202241x	return 2 * width() + 2 * height();
64		}
65
66	401621x	bool intersects(const Box & other) const {
67	401621x	if (other.ymin > ymax)
68	42x	return false;
69	401579x	if (other.ymax < ymin)
70	51x	return false;
71	401528x	if (other.xmin > xmax)
72	!	return false;
73	401528x	if (other.xmax < xmin)
74	!	return false;
75
76	401528x	return true;
77		}
78
79	4101192x	Box intersection(const Box & other) const {
80		return {
81	4101192x	std::max(xmin, other.xmin),
82	4101192x	std::max(ymin, other.ymin),
83	4101192x	std::min(xmax, other.xmax),
84	4101192x	std::min(ymax, other.ymax)
85		};
86		}
87
88		Box translate(double dx, double dy) const {
89		return {
90		xmin + dx,
91		ymin + dy,
92		xmax + dx,
93		ymax + dy
94		};
95		}
96
97		Coordinate upper_left() const {
98		return Coordinate{xmin, ymax};
99		}
100
101		Coordinate upper_right() const {
102		return Coordinate{xmax, ymax};
103		}
104
105		Coordinate lower_left() const {
106		return Coordinate{xmin, ymin};
107		}
108
109		Coordinate lower_right() const {
110		return Coordinate{xmax, ymin};
111		}
112
113		Side side(const Coordinate &c) const;
114
115		Crossing crossing(const Coordinate &c1, const Coordinate &c2) const;
116
117	517462x	bool empty() const {
118	517462x	return xmin >= xmax \|\| ymin >= ymax;
119		}
120
121	13x	Box expand_to_include(const Box & other) const {
122	13x	if (empty()) {
123	!	return other;
124		}
125
126	13x	if (other.empty()) {
127	!	return *this;
128		}
129
130	13x	return { std::min(xmin, other.xmin),
131	13x	std::min(ymin, other.ymin),
132	13x	std::max(xmax, other.xmax),
133	13x	std::max(ymax, other.ymax) };
134		}
135
136		bool contains(const Box &b) const;
137
138		bool contains(const Coordinate &c) const;
139
140		bool strictly_contains(const Coordinate &c) const;
141
142	400883x	bool operator==(const Box& other) const {
143	400883x	return xmin == other.xmin && xmax == other.xmax && ymin == other.ymin && ymax == other.ymax;
144		}
145
146		friend std::ostream &operator<<(std::ostream &os, const Box &c);
147		};
148
149		std::ostream &operator<<(std::ostream &os, const Box &c);
150
151		}
152
153		#endif

1		// Copyright (c) 2021 ISciences, LLC.
2		// All rights reserved.
3		//
4		// This software is licensed under the Apache License, Version 2.0 (the "License").
5		// You may not use this file except in compliance with the License. You may
6		// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		//
8		// Unless required by applicable law or agreed to in writing, software
9		// distributed under the License is distributed on an "AS IS" BASIS,
10		// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		// See the License for the specific language governing permissions and
12		// limitations under the License.
13
14		#ifndef EXACTEXTRACT_RASTER_AREA_H
15		#define EXACTEXTRACT_RASTER_AREA_H
16
17		#include "raster.h"
18
19		namespace exactextract {
20		template<typename T>
21		class CartesianAreaRaster : public AbstractRaster<T> {
22		public:
23	14x	explicit CartesianAreaRaster(const Grid<bounded_extent> & ex) :
24		AbstractRaster<T>(ex),
25	14x	m_area(grid_cell(AbstractRaster<T>::grid(), 0, 0).area())
26		{}
27
28	7116x	T operator()(size_t row, size_t column) const override {
29		(void) row;
30		(void) column;
31	7116x	return m_area;
32		}
33		private:
34		double m_area;
35		};
36
37		template<typename T>
38		class SphericalAreaRaster : public AbstractRaster<T> {
39		public:
40	25x	explicit SphericalAreaRaster(const Grid<bounded_extent> &ex) :
41		AbstractRaster<T>(ex),
42	25x	m_areas(ex.rows())
43		{
44	25x	const auto& g = AbstractRaster<T>::grid();
45
46	25x	double dlon = g.dx();
47	25x	double dlat = g.dy();
48
49	99x	for (size_t i = 0; i < g.rows(); i++) {
50	74x	double y = g.y_for_row(i);
51	74x	double ymin = y - 0.5 * dlat;
52	74x	double ymax = y + 0.5 * dlat;
53
54	74x	m_areas[i] = EARTH_RADIUS_SQ * PI_180 * std::abs(std::sin(ymin * PI_180) - std::sin(ymax * PI_180)) * dlon;
55		}
56		}
57
58	276x	T operator()(size_t row, size_t column) const override {
59		(void) column;
60	276x	return m_areas[row];
61		}
62
63		private:
64		static constexpr double EARTH_RADIUS = 6378137;
65		static constexpr double EARTH_RADIUS_SQ = EARTH_RADIUS * EARTH_RADIUS;
66		static constexpr double PI_180 = 3.141592653589793238462643383279502884197169399375105 / 180;
67
68		std::vector<double> m_areas;
69		};
70		}
71
72		#endif //EXACTEXTRACT_RASTER_AREA_H