exactextractr coverage - 90.75%

Files
Source

# 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("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`
#'
#' 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     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     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
                           ) {
  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')
  .validateNumericScalarOrNA(default_value, 'default_value')
  .validateNumericScalarOrNA(default_weight, 'default_weight')

  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 (progress && length(geoms) > 1) {
    n <- length(geoms)
    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() {}
  }

  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 <- sapply(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, quantiles)
        update_progress()
        return(ret)
      })

      if (!is.matrix(results) && !force_df) {
        # Single stat? Return a vector unless asked otherwise via force_df.
        return(results)
      } else {
        # Return a data frame with a column for each stat
        colnames <- .resultColNames(names(x), names(weights), fun, full_colnames, quantiles)

        if (is.matrix(results)) {
          results <- t(results)
        } else {
          results <- matrix(results, nrow=length(results))
        }

        dimnames(results) <- list(NULL, colnames)
        ret <- as.data.frame(results)

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

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

        return(ret)
      }
    } 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)
        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,
                                      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)
        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) 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.

#' 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
.resultColNames <- function(value_names, weight_names, fun, full_colnames, quantiles=numeric()) {
  if (inherits(fun, 'standardGeneric')) {
    stat_names <- fun@generic[1]
  } else if (is.function(fun)) {
    stat_names <- 'fun'
  } else {
    stat_names <- fun
  }

  quantile_index = which(stat_names == 'quantile')
  if (length(quantile_index) != 0) {
    stat_names <- c(stat_names[seq_along(stat_names) < quantile_index],
                    sprintf('q%02d', as.integer(100 * quantiles)),
                    stat_names[seq_along(stat_names) > quantile_index])
  }

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

  if (length(vn) > 1 || full_colnames) {
    # determine all combinations of index and stat
    z <- expand.grid(index=seq_along(vn),
                     stat=stat_names, stringsAsFactors=FALSE)
    z$value <- vn[z$index]
    if (is.null(wn)) {
      z$weights <- NA
    } else {
      z$weights <- wn[z$index]
    }

    # construct column names for each index, stat
    # add weight layer name only if layer is ambiguously weighted
    mapply(function(stat, value, weight) {
      ret <- stat
      if (full_colnames || length(value_names) > 1) {
        ret <- paste(ret, value, sep='.')
      }
      if (.includeWeightInColName(stat) && ((full_colnames & length(weight_names) > 0)
                                            || length(weight_names) > 1)) {
        ret <- paste(ret, weight, sep='.')
      }

      return(ret)
    }, z$stat, z$value, z$weight, USE.NAMES = FALSE)
  } else {
    stat_names
  }
}

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

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

#' 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')
  }
}

# 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))
  }
}

# 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.

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)

# Copyright (c) 2020-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("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 to be used for the resampling
#' @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) {
  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.")
    } 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 (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.")
  }

  x <- .startReading(x)
  tryCatch({
    ret <- CPP_resample(x, y, fun)
    if (class(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-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.

#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 raster = Rcpp::Environment::namespace_env("raster");
    Rcpp::Environment xx = Rcpp::Environment::namespace_env("exactextractr");
    Rcpp::Function rasterFn = raster["raster"];
    Rcpp::Function crsFn = xx[".crs"];

    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;
        }
      }
    }

    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::floor(disaggregation_factor_x)) > 1e-6 ||
                    std::abs(disaggregation_factor_y - std::floor(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>(disaggregation_factor_x);
                m_ry = static_cast<size_t>(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-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.

// [[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> & 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) {
  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);

      if (!weights_grid.compatible_with(grid)) {
        Rcpp::stop("Incompatible extents.");
      }

      common_grid = grid.common_grid(weights_grid);

      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, common_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() != common_grid.dx() || grid.dy() != common_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, common_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() != common_grid.dx() || weights_grid.dy() != common_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) {
      cols["cell"] = get_cell_numbers(rast, cov_grid)[covered];
    }

    if (include_area) {
      auto area_rast = get_area_raster(area_method, common_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
};

// 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,
                              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()) : rsrc.grid();

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

    bool store_values = false;
    int stat_result_size = 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("quantile")) {
        int num_quantiles = 0;

        if (quantiles.isNotNull()) {
          Rcpp::NumericVector qvec = quantiles.get();
          num_quantiles = qvec.size();
        }

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

        stat_result_size += num_quantiles;
      } else {
        stat_result_size += 1;
      }
    }

    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);
    }

    Rcpp::NumericMatrix stat_results = Rcpp::no_init(nresults, stat_result_size);

    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);
              }
            }
          }
        }
      }
    }

    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 Rcpp::stop("Unknown stat: " + stat);
      }
    }

    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-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_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:
        /**
         * Compute raster statistics from a Raster representing intersection percentages,
         * a Raster representing data values, and (optionally) a Raster representing weights.
         * and a set of raster values.
         */
        explicit RasterStats(bool store_values = false) :
                m_min{std::numeric_limits<T>::max()},
                m_max{std::numeric_limits<T>::lowest()},
                m_sum_ciwi{0},
                m_sum_ci{0},
                m_sum_xici{0},
                m_sum_xiciwi{0},
                m_store_values{store_values} {}

        void process(const Raster<float> & intersection_percentages, const AbstractRaster<T> & rast) {
            std::unique_ptr<AbstractRaster<T>> rvp;

            if (rast.grid() != intersection_percentages.grid()) {
                rvp = std::make_unique<RasterView<T>>(rast, intersection_percentages.grid());
            }

            const AbstractRaster<T>& rv = rvp ? *rvp : rast;

            for (size_t i = 0; i < rv.rows(); i++) {
                for (size_t j = 0; j < rv.cols(); j++) {
                    float pct_cov = intersection_percentages(i, j);
                    T val;
                    if (pct_cov > 0 && rv.get(i, j, val)) {
                        process_value(val, pct_cov, 1.0);
                    }
                }
            }
        }

        void process(const Raster<float> & intersection_percentages, const AbstractRaster<T> & rast, const AbstractRaster<T> & weights) {
            // Process the entire intersection_percentages grid, even though it may
            // be outside the extent of the weighting raster. Although we've been
            // provided a weighting raster, we still need to calculate correct values
            // for unweighted stats.
            auto& common = intersection_percentages.grid();

            if (common.empty())
                return;

            // If the value or weights grids do not correspond to the intersection_percentages grid,
            // construct a RasterView to perform the transformation. Even a no-op RasterView can be
            // expensive, so we avoid doing this unless necessary.
            std::unique_ptr<AbstractRaster<T>> rvp;
            std::unique_ptr<AbstractRaster<T>> wvp;

            if (rast.grid() != common) {
                rvp = std::make_unique<RasterView<T>>(rast, common);
            }

            if (weights.grid() != common) {
                wvp = std::make_unique<RasterView<T>>(weights, common);
            }

            const AbstractRaster<T>& rv = rvp ? *rvp : rast;
            const AbstractRaster<T>& wv = wvp ? *wvp : weights;

            for (size_t i = 0; i < rv.rows(); i++) {
                for (size_t j = 0; j < rv.cols(); j++) {
                    float pct_cov = intersection_percentages(i, j);
                    T weight;
                    T val;

                    if (pct_cov > 0 && rv.get(i, j, val)) {
                        if (wv.get(i, j, weight)) {
                            process_value(val, pct_cov, weight);
                        } else {
                            // Weight is NODATA, convert to NAN
                            process_value(val, pct_cov, std::numeric_limits<double>::quiet_NaN());
                        }
                    }
                }
            }
        }

        /**
         * The mean value of cells covered by this polygon, weighted
         * by the percent of the cell that is covered.
         */
        float mean() const {
            return sum() / count();
        }

        /**
         * The mean value of cells covered by this polygon, weighted
         * by the percent of the cell that is covered and a secondary
         * weighting raster.
         *
         * If any weights are undefined, will return NAN. If this is undesirable,
         * caller should replace undefined weights with a suitable default
         * before computing statistics.
         */
         float weighted_mean() const {
             return weighted_sum() / weighted_count();
         }

         /** The fraction of weighted cells to unweighted cells.
          *  Meaningful only when the values of the weighting
          *  raster are between 0 and 1.
          */
         float weighted_fraction() const {
             return weighted_sum() / sum();
         }

        /**
         * The raster value occupying the greatest number of cells
         * or partial cells within the polygon. When multiple values
         * cover the same number of cells, the greatest value will
         * be returned. Weights are not taken into account.
         */
        nonstd::optional<T> mode() const {
            if (variety() == 0) {
                return nonstd::nullopt;
            }

            return std::max_element(m_freq.cbegin(),
                                    m_freq.cend(),
                                    [](const auto &a, const auto &b) {
                                        return a.second < b.second || (a.second == b.second && a.first < b.first);
                                    })->first;
        }

        /**
         * The minimum value in any raster cell wholly or partially covered
         * by the polygon. Weights are not taken into account.
         */
        nonstd::optional<T> min() const {
            if (m_sum_ci == 0) {
                return nonstd::nullopt;
            }
            return m_min;
        }

        /**
         * The maximum value in any raster cell wholly or partially covered
         * by the polygon. Weights are not taken into account.
         */
        nonstd::optional<T> max() const {
            if (m_sum_ci == 0) {
                return nonstd::nullopt;
            }
            return m_max;
        }

        /**
         * The given quantile (0-1) of raster cell values. Coverage fractions
         * are taken into account but weights are not.
         */
        nonstd::optional<T> quantile(double q) const {
            if (m_sum_ci == 0) {
                return nonstd::nullopt;
            }

            // The weighted quantile computation is not processed incrementally.
            // Create it on demand and retain it in case we want multiple quantiles.
            if (!m_quantiles) {
                m_quantiles = std::make_unique<WeightedQuantiles>();

                for (const auto& entry : m_freq) {
                    m_quantiles->process(entry.first, entry.second);
                }
            }

            return m_quantiles->quantile(q);
        }

        /**
         * The sum of raster cells covered by the polygon, with each raster
         * value weighted by its coverage fraction.
         */
        float sum() const {
            return (float) m_sum_xici;
        }

        /**
         * The sum of raster cells covered by the polygon, with each raster
         * value weighted by its coverage fraction and weighting raster value.
         *
         * If any weights are undefined, will return NAN. If this is undesirable,
         * caller should replace undefined weights with a suitable default
         * before computing statistics.
         */
        float weighted_sum() const {
            return (float) m_sum_xiciwi;
        }

        /**
         * The number of raster cells with a defined value
         * covered by the polygon. Weights are not taken
         * into account.
         */
        float count() const {
            return (float) m_sum_ci;
        }

        /**
         * The population variance of raster cells touched
         * by the polygon. Cell coverage fractions are taken
         * into account; values of a weighting raster are not.
         */
        float variance() const {
            return static_cast<float>(m_variance.variance());
        }

        /**
         * The population standard deviation of raster cells
         * touched by the polygon. Cell coverage fractions
         * are taken into account; values of a weighting
         * raster are not.
         */
        float stdev() const {
            return static_cast<float>(m_variance.stdev());
        }

        /**
         * The population coefficient of variation of raster
         * cells touched by the polygon. Cell coverage fractions
         * are taken into account; values of a weighting
         * raster are not.
         */
        float coefficient_of_variation() const {
            return static_cast<float>(m_variance.coefficent_of_variation());
        }

        /**
         * The sum of weights for each cell covered by the
         * polygon, with each weight multiplied by the coverage
         * coverage fraction of each cell.
         *
         * If any weights are undefined, will return NAN. If this is undesirable,
         * caller should replace undefined weights with a suitable default
         * before computing statistics.
         */
        float weighted_count() const {
            return (float) m_sum_ciwi;
        }

        /**
         * The raster value occupying the least number of cells
         * or partial cells within the polygon. When multiple values
         * cover the same number of cells, the lowest value will
         * be returned.
         *
         * Cell weights are not taken into account.
         */
        nonstd::optional<T> minority() const {
            if (variety() == 0) {
                return nonstd::nullopt;
            }

            return std::min_element(m_freq.cbegin(),
                                    m_freq.cend(),
                                    [](const auto &a, const auto &b) {
                                        return a.second < b.second || (a.second == b.second && a.first < b.first);
                                    })->first;
        }

        /**
         * The number of distinct defined raster values in cells wholly
         * or partially covered by the polygon.
         */
        size_t variety() const {
            return m_freq.size();
        }

        bool stores_values() const {
            return m_store_values;
        }

    private:
        T m_min;
        T m_max;

        // ci: coverage fraction of pixel i
        // wi: weight of pixel i
        // xi: value of pixel i
        double m_sum_ciwi;
        double m_sum_ci;
        double m_sum_xici;
        double m_sum_xiciwi;
        WestVariance m_variance;

        mutable std::unique_ptr<WeightedQuantiles> m_quantiles;

        std::unordered_map<T, float> m_freq;

        bool m_store_values;

        void process_value(const T& val, float coverage, double weight) {
            m_sum_ci += static_cast<double>(coverage);
            m_sum_xici += val*static_cast<double>(coverage);

            m_variance.process(val, coverage);

            double ciwi = static_cast<double>(coverage)*weight;
            m_sum_ciwi += ciwi;
            m_sum_xiciwi += val * ciwi;

            if (val < m_min) {
                m_min = val;
            }

            if (val > m_max) {
                m_max = val;
            }

            // TODO should weights factor in here?
            if (m_store_values) {
                m_freq[val] += coverage;
                m_quantiles.reset();
            }
        }
    };

    template<typename T>
    std::ostream& operator<<(std::ostream& os, const RasterStats<T> & stats) {
        os << "{" << std::endl;
        os << "  \"count\" : " << stats.count() << "," << std::endl;

        os << "  \"min\" : ";
        if (stats.min().has_value()) {
            os << stats.min().value();
        } else {
            os << "null";
        }
        os << "," << std::endl;

        os << "  \"max\" : ";
        if (stats.max().has_value()) {
            os << stats.max().value();
        } else {
            os << "null";
        }
        os << "," << std::endl;

        os << "  \"mean\" : " << stats.mean() << "," << std::endl;
        os << "  \"sum\" : " << stats.sum() << "," << std::endl;
        os << "  \"weighted_mean\" : " << stats.weighted_mean() << "," << std::endl;
        os << "  \"weighted_sum\" : " << stats.weighted_sum();
        if (stats.stores_values()) {
            os << "," << std::endl;
            os << "  \"mode\" : ";
            if (stats.mode().has_value()) {
                os << stats.mode().value();
            } else {
                os << "null";
            }
            os << "," << std::endl;

            os << "  \"minority\" : ";
            if (stats.minority().has_value()) {
                os << stats.minority().value();
            } else {
                os << "null";
            }
            os << "," << std::endl;

            os << "  \"variety\" : " << stats.variety() << std::endl;
        } else {
            os << std::endl;
        }
        os << "}" << std::endl;
        return os;
    }


}

#endif

// Copyright (c) 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_VARIANCE_H
#define EXACTEXTRACT_VARIANCE_H

#include <cmath>

namespace exactextract {

class WestVariance {
    /** \brief Implements an incremental algorithm for weighted standard
     * deviation, variance, and coefficient of variation, as described in
     * formula WV2 of West, D.H.D. (1979) "Updating Mean and Variance
     * Estimates: An Improved Method". Communications of the ACM 22(9).
     */

private:
    double sum_w = 0;
    double mean = 0;

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		#' Return column names to be used for summary operations
15		#'
16		#' @param value_names names of value raster layers
17		#' @param weight_names names of weighting raster layers
18		#' @param fun functions or names of summary operations
19		#' @param full_colnames return a complete column name even when there is no
20		#' ambiguity?
21		#' @param quantiles quantiles to use when \code{stat_names} contains \code{quantile}
22		#' @return character vector of column names
23		#' @keywords internal
24		.resultColNames <- function(value_names, weight_names, fun, full_colnames, quantiles=numeric()) {
25	120x	if (inherits(fun, 'standardGeneric')) {
26	22x	stat_names <- fun@generic[1]
27	98x	} else if (is.function(fun)) {
28	32x	stat_names <- 'fun'
29		} else {
30	66x	stat_names <- fun
31		}
32
33	120x	quantile_index = which(stat_names == 'quantile')
34	120x	if (length(quantile_index) != 0) {
35	1x	stat_names <- c(stat_names[seq_along(stat_names) < quantile_index],
36	1x	sprintf('q%02d', as.integer(100 * quantiles)),
37	1x	stat_names[seq_along(stat_names) > quantile_index])
38		}
39
40	120x	ind <- .valueWeightIndexes(length(value_names), length(weight_names))
41	120x	vn <- value_names[ind$values]
42	120x	wn <- weight_names[ind$weights]
43
44	120x	if (length(vn) > 1 \|\| full_colnames) {
45		# determine all combinations of index and stat
46	37x	z <- expand.grid(index=seq_along(vn),
47	37x	stat=stat_names, stringsAsFactors=FALSE)
48	37x	z$value <- vn[z$index]
49	37x	if (is.null(wn)) {
50	17x	z$weights <- NA
51		} else {
52	20x	z$weights <- wn[z$index]
53		}
54
55		# construct column names for each index, stat
56		# add weight layer name only if layer is ambiguously weighted
57	37x	mapply(function(stat, value, weight) {
58	132x	ret <- stat
59	132x	if (full_colnames \|\| length(value_names) > 1) {
60	117x	ret <- paste(ret, value, sep='.')
61		}
62	132x	if (.includeWeightInColName(stat) && ((full_colnames & length(weight_names) > 0)
63	132x	\|\| length(weight_names) > 1)) {
64	48x	ret <- paste(ret, weight, sep='.')
65		}
66
67	132x	return(ret)
68	37x	}, z$stat, z$value, z$weight, USE.NAMES = FALSE)
69		} else {
70	83x	stat_names
71		}
72		}
73
74		.includeWeightInColName <- function(fun) {
75	132x	.isWeighted(fun) \| fun == 'fun'
76		}
77
78		.isWeighted <- function(stat_name) {
79	211x	stat_name %in% c('weighted_mean', 'weighted_sum')
80		}
81
82		#' Compute indexes for the value and weight layers that should be
83		#' processed together
84		#'
85		#' @param num_values number of layers in value raster
86		#' @param num_weights number of layers in weighting raster
87		#' @return list with \code{values} and \code{weights} elements
88		#' providing layer indexes
89		#' @keywords internal
90		.valueWeightIndexes <- function(num_values, num_weights) {
91	196x	if (num_weights == 0) {
92	141x	vi <- seq_len(num_values)
93	141x	wi <- NA
94	55x	} else if (num_values == num_weights) {
95		# process in parallel
96	40x	vi <- seq_len(num_values)
97	40x	wi <- seq_len(num_weights)
98	15x	} else if (num_values == 1 && num_weights > 1) {
99		# recycle values
100	7x	vi <- rep.int(1, num_weights)
101	7x	wi <- seq_len(num_weights)
102	8x	} else if (num_values > 1 && num_weights == 1) {
103		# recycle weights
104	8x	vi <- seq_len(num_values)
105	8x	wi <- rep.int(1, num_values)
106		}
107
108	196x	list(values = vi, weights = wi)
109		}
110
111		.areaMethod <- function(crs_obj) {
112	21x	if (!(is.na(crs_obj)) && crs_obj$units_gdal == 'degree') {
113	7x	return('spherical')
114		} else {
115	14x	return('cartesian')
116		}
117		}
118
119		.validateFlag <- function(value, name) {
120	2574x	if(!(is.logical(value) && length(value) == 1 && !is.na(value))) {
121	27x	stop(name, ' must be TRUE or FALSE')
122		}
123		}
124
125		.validateNumericScalar <- function(value, name) {
126	271x	if (!(is.numeric(value) && length(value) == 1 && !is.na(value))) {
127	4x	stop(name, ' must be a single numeric value')
128		}
129		}
130
131		.validateNumericScalarOrNA <- function(value, name) {
132	530x	if (!(is.numeric(value) && length(value) == 1)) {
133	4x	stop(name, ' must be a single numeric value')
134		}
135		}
136
137		# faster replacement for as.data.frame when input is a named list
138		# with equal-length columns
139		# from Advanced R, sec. 24.4.2
140		.quickDf <- function(lst) {
141	139x	class(lst) <- 'data.frame'
142	139x	attr(lst, 'row.names') <- .set_row_names(length(lst[[1]]))
143	139x	lst
144		}
145
146		.singleColumnToVector <- function(df) {
147	14x	if (ncol(df) == 1) {
148	3x	df[, 1]
149		} else {
150	11x	df
151		}
152		}
153
154		# Return the number of standard (non-...) arguments in a supplied function that
155		# do not have a default value. This is used to fail if the summary function
156		# provided by the user cannot accept arguments of values and weights.
157		.num_expected_args <- function(fun) {
158	63x	a <- formals(args(fun))
159	63x	a <- a[names(a) != '...']
160	63x	sum(sapply(a, nchar) == 0)
161		}
162
163		.startReading <- function(r) {
164	475x	if(inherits(r, 'BasicRaster')) {
165	306x	return(raster::readStart(r))
166	169x	} else if (inherits(r, 'SpatRaster')) {
167	14x	terra::readStart(r)
168		}
169
170	169x	return(r)
171		}
172
173		.stopReading <- function(r) {
174	475x	if(inherits(r, 'BasicRaster')) {
175	306x	return(raster::readStop(r))
176	169x	} else if (inherits(r, 'SpatRaster')) {
177	14x	terra::readStop(r)
178		}
179
180	169x	return(r)
181		}
182
183		.crs <- function(r) {
184	12x	if(inherits(r, 'BasicRaster')) {
185	12x	if (utils::packageVersion('raster') < numeric_version('3.5')) {
186	!	return(raster::crs(r))
187		} else {
188	12x	return(terra::crs(r))
189		}
190	!	} else if (inherits(r, 'SpatRaster')) {
191	!	return(terra::crs(r))
192		} else {
193	!	stop('Unknown type: ', class(r))
194		}
195		}
196
197		.setValues <- function(r, x) {
198	7x	if(inherits(r, 'BasicRaster')) {
199	7x	if (utils::packageVersion('raster') < numeric_version('3.5')) {
200	!	raster::values(r) <- x
201		} else {
202	7x	terra::values(r) <- x
203		}
204	!	} else if (inherits(r, 'SpatRaster')) {
205	!	raster::values(r) <- x
206		} else {
207	!	stop('Unknown type: ', class(r))
208		}
209
210	7x	return(r)
211		}
212
213		.numLayers <- function(r) {
214	501x	if(inherits(r, 'BasicRaster')) {
215	481x	return(raster::nlayers(r))
216	20x	} else if (inherits(r, 'SpatRaster')) {
217	20x	return(terra::nlyr(r))
218		} else {
219	!	stop('Unknown type: ', class(r))
220		}
221		}
222
223		.isRaster <- function(r) {
224	179x	inherits(r, 'BasicRaster') \| inherits(r, 'SpatRaster')
225		}
226
227		.xFromCol <- function(r, col) {
228	17x	if (inherits(r, 'BasicRaster')) {
229	16x	raster::xFromCol(r, col)
230	1x	} else if (inherits(r, 'SpatRaster')) {
231	1x	terra::xFromCol(r, col)
232		} else {
233	!	stop('Unknown type: ', class(r))
234		}
235		}
236
237		.colFromX <- function(r, x) {
238	28x	if (inherits(r, 'BasicRaster')) {
239	26x	raster::colFromX(r, x)
240	2x	} else if (inherits(r, 'SpatRaster')) {
241	2x	terra::colFromX(r, x)
242		} else {
243	!	stop('Unknown type: ', class(r))
244		}
245		}
246
247		.yFromRow <- function(r, row) {
248	17x	if (inherits(r, 'BasicRaster')) {
249	16x	raster::yFromRow(r, row)
250	1x	} else if (inherits(r, 'SpatRaster')) {
251	1x	terra::yFromRow(r, row)
252		} else {
253	!	stop('Unknown type: ', class(r))
254		}
255		}
256
257		.rowFromY <- function(r, y) {
258	28x	if (inherits(r, 'BasicRaster')) {
259	26x	raster::rowFromY(r, y)
260	2x	} else if (inherits(r, 'SpatRaster')) {
261	2x	terra::rowFromY(r, y)
262		} else {
263	!	stop('Unknown type: ', class(r))
264		}
265		}
266
267		.cellFromRowCol <- function(r, row, col) {
268	11x	if (inherits(r, 'BasicRaster')) {
269	10x	raster::cellFromRowCol(r, row, col)
270	1x	} else if (inherits(r, 'SpatRaster')) {
271	1x	terra::cellFromRowCol(r, row, col)
272		} else {
273	!	stop('Unknown type: ', class(r))
274		}
275		}
276
277		.extent <- function(r) {
278	569x	if (inherits(r, 'BasicRaster')) {
279	546x	ex <- r@extent
280	546x	c(ex@xmin, ex@ymin, ex@xmax, ex@ymax)
281	23x	} else if (inherits(r, 'SpatRaster')) {
282	23x	ex <- terra::ext(r)
283	23x	c(ex$xmin, ex$ymin, ex$xmax, ex$ymax)
284		} else {
285	!	stop('Unknown type: ', class(r))
286		}
287		}
288
289		.res <- function(r) {
290	569x	if (inherits(r, 'BasicRaster')) {
291	546x	raster::res(r)
292	23x	} else if (inherits(r, 'SpatRaster')) {
293	23x	terra::res(r)
294		} else {
295	!	stop('Unknown type: ', class(r))
296		}
297		}
298
299		.getValuesBlock <- function(r, row, nrows, col, ncols) {
300	946x	if (inherits(r, 'BasicRaster')) {
301	928x	raster::getValuesBlock(r, row, nrows, col, ncols, format = 'm')
302	18x	} else if (inherits(r, 'SpatRaster')) {
303	18x	terra::readValues(r, row, nrows, col, ncols, mat = TRUE)
304		} else {
305	!	stop('Unknown type: ', class(r))
306		}
307		}

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	12x	setGeneric('coverage_fraction', function(x, y, crop=FALSE, ...)
15	12x	standardGeneric('coverage_fraction'))
16
17		.coverage_fraction <- function(x, y, crop) {
18	12x	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	11x	} 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	10x	} 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	12x	lapply(sf::st_as_binary(y, EWKB=TRUE), function(wkb) {
28	12x	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)

1		# Copyright (c) 2020-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	10x	setGeneric("exact_resample", function(x, y, ...)
15	10x	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 to be used for the resampling
23		#' @return a resampled version of \code{x}, returned as a \code{RasterLayer} or
24		#' \code{SpatRaster}, depending on the values of \code{x} and \code{y}
25		#'
26		#' @name exact_resample
27		NULL
28
29		.exact_resample <- function(x, y, fun) {
30	10x	if (sf::st_crs(x) != sf::st_crs(y)) {
31	3x	if (is.na(sf::st_crs(x))) {
32	1x	warning("No CRS specified for source raster; assuming it has the same CRS as destination raster.")
33	2x	} else if (is.na(sf::st_crs(y))) {
34	1x	warning("No CRS specified for destination raster; assuming it has the same CRS as source raster.")
35		} else {
36	1x	stop('Destination raster must have same CRS as source.')
37		}
38		}
39
40	9x	if (length(fun) != 1) {
41	1x	stop("Only a single operation may be used for resampling.")
42		}
43
44	8x	if (startsWith(fun, 'weighted')) {
45	1x	stop("Weighted operations cannot be used for resampling.")
46		}
47
48	7x	x <- .startReading(x)
49	7x	tryCatch({
50	7x	ret <- CPP_resample(x, y, fun)
51	7x	if (class(x) == 'SpatRaster') {
52	1x	terra::rast(ret)
53		} else {
54	6x	ret
55		}
56	7x	}, finally={
57	7x	.stopReading(x)
58		})
59		}
60
61		#' @import sf
62		#' @import raster
63		#' @useDynLib exactextractr
64		#' @rdname exact_resample
65		#' @export
66		setMethod('exact_resample', signature(x='RasterLayer', y='RasterLayer'), .exact_resample)
67
68		#' @useDynLib exactextractr
69		#' @rdname exact_resample
70		#' @export
71		setMethod('exact_resample', signature(x='SpatRaster', y='SpatRaster'), .exact_resample)

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		#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	12x	Rcpp::S4 CPP_coverage_fraction(Rcpp::S4 & rast, const Rcpp::RawVector & wkb, bool crop)
28		{
29		try {
30	24x	GEOSAutoHandle geos;
31	36x	Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster");
32	36x	Rcpp::Environment xx = Rcpp::Environment::namespace_env("exactextractr");
33	36x	Rcpp::Function rasterFn = raster["raster"];
34	36x	Rcpp::Function crsFn = xx[".crs"];
35
36	12x	auto grid = make_grid(rast);
37	24x	auto coverage_fraction = raster_cell_intersection(grid, geos.handle, read_wkb(geos.handle, wkb).get());
38
39	12x	if (crop) {
40	2x	grid = coverage_fraction.grid();
41		}
42	24x	RasterView<float> coverage_view(coverage_fraction, grid);
43
44	12x	Rcpp::NumericMatrix weights{static_cast<int>(grid.rows()),
45	12x	static_cast<int>(grid.cols())};
46
47	629x	for (size_t i = 0; i < grid.rows(); i++) {
48	159018x	for (size_t j = 0; j < grid.cols(); j++) {
49	158401x	weights(i, j) = coverage_view(i, j);
50	158401x	if (!crop && std::isnan(weights(i, j))) {
51	139692x	weights(i, j) = 0;
52		}
53		}
54		}
55
56		return rasterFn(weights,
57	24x	Rcpp::Named("xmn")=grid.xmin(),
58	24x	Rcpp::Named("xmx")=grid.xmax(),
59	24x	Rcpp::Named("ymn")=grid.ymin(),
60	24x	Rcpp::Named("ymx")=grid.ymax(),
61	48x	Rcpp::Named("crs")=crsFn(rast));
62		} catch (std::exception & e) {
63		// throw predictable exception class
64		#ifdef __SUNPRO_CC
65		// Rcpp::stop crashes CRAN Solaris build
66		// https://github.com/RcppCore/Rcpp/issues/1159
67		Rf_error(e.what());
68		return R_NilValue;
69		#else
70		Rcpp::stop(e.what());
71		#endif
72		}
73		}

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	411x	static inline bool is_integral(double d, double tol) {
33	411x	return std::abs(d - std::round(d)) <= tol;
34		}
35
36		template<typename extent_tag>
37		class Grid {
38
39		public:
40	1835108x	Grid(const Box & extent, double dx, double dy) :
41		m_extent{extent},
42		m_dx{dx},
43		m_dy{dy},
44	1835108x	m_num_rows{2*extent_tag::padding + (extent.ymax > extent.ymin ? static_cast<size_t>(std::round((extent.ymax - extent.ymin) / dy)) : 0)},
45	3670216x	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	118250x	static Grid make_empty() {
49	118250x	return Grid({0, 0, 0, 0}, 0, 0);
50		}
51
52	4792155x	size_t get_column(double x) const {
53		if (extent_tag::padding) {
54	6382908x	if (x < m_extent.xmin)
55	878x	return 0;
56	6382030x	if (x > m_extent.xmax)
57	1180x	return m_num_cols - 1;
58	6380850x	if (x == m_extent.xmax) // special case, returning the cell for which xmax is the right
59	3191314x	return m_num_cols - 2;
60		} else {
61	1600701x	if (x < m_extent.xmin \|\| x > m_extent.xmax)
62	!	throw std::out_of_range("x");
63
64	1600701x	if (x == m_extent.xmax)
65	801062x	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	7183221x	return std::min(
72	2394407x	extent_tag::padding + static_cast<size_t>(std::floor((x - m_extent.xmin) / m_dx)),
73	4788814x	get_column(m_extent.xmax));
74		}
75
76	4790894x	size_t get_row(double y) const {
77		if (extent_tag::padding) {
78	6380254x	if (y > m_extent.ymax)
79	2246x	return 0;
80	6378008x	if (y < m_extent.ymin)
81	1888x	return m_num_rows - 1;
82	6376120x	if (y == m_extent.ymin) // special case, returning the cell for which ymin is the bottom
83	3189238x	return m_num_rows - 2;
84		} else {
85	1600767x	if (y < m_extent.ymin \|\| y > m_extent.ymax)
86	!	throw std::out_of_range("y");
87
88	1600767x	if (y == m_extent.ymin)
89	801062x	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	7179438x	return std::min(
96	2393146x	extent_tag::padding + static_cast<size_t>(std::floor((m_extent.ymax - y) / m_dy)),
97	4786292x	get_row(m_extent.ymin));
98		}
99
100	2217181x	bool empty() const { return m_num_rows <= 2extent_tag::padding && m_num_cols <= 2extent_tag::padding; }
101
102	10645683x	size_t rows() const { return m_num_rows; }
103
104	12130535x	size_t cols() const { return m_num_cols; }
105
106	159x	size_t size() const { return rows()*cols(); }
107
108	5318573x	double xmin() const { return m_extent.xmin; }
109
110	3810145x	double xmax() const { return m_extent.xmax; }
111
112	3668372x	double ymin() const { return m_extent.ymin; }
113
114	5553910x	double ymax() const { return m_extent.ymax; }
115
116	5834397x	double dx() const { return m_dx; }
117
118	6071555x	double dy() const { return m_dy; }
119
120	2235128x	const Box& extent() const { return m_extent; }
121
122	400188x	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	400188x	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	409x	double x_for_col(size_t col) const { return m_extent.xmin + ((col - extent_tag::padding) + 0.5) * m_dx; }
127
128	1051x	double y_for_row(size_t row) const { return m_extent.ymax - ((row - extent_tag::padding) + 0.5) * m_dy; }
129
130	1377x	Grid<extent_tag> crop(const Box & b) const {
131	1377x	if (extent().intersects(b)) {
132	1294x	return shrink_to_fit(b.intersection(extent()));
133		} else {
134	83x	return make_empty();
135		}
136		}
137
138	400531x	Grid<extent_tag> shrink_to_fit(const Box & b) const {
139	400531x	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	400531x	size_t col0 = get_column(b.xmin);
144	400531x	size_t row1 = get_row(b.ymax);
145
146		// Shrink xmin and ymax to fit the upper-left corner of the supplied extent
147	400531x	double snapped_xmin = m_extent.xmin + (col0 - extent_tag::padding) * m_dx;
148	400531x	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	400531x	if (b.xmin < snapped_xmin) {
153	!	snapped_xmin -= m_dx;
154	!	col0--;
155		}
156
157	400531x	if (b.ymax > snapped_ymax) {
158	!	snapped_ymax += m_dy;
159	!	row1--;
160		}
161
162	400531x	size_t col1 = get_column(b.xmax);
163	400531x	size_t row0 = get_row(b.ymin);
164
165	400531x	size_t num_rows = 1 + (row0 - row1);
166	400531x	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	400531x	if (num_rows > 2 && (snapped_ymax - (num_rows-1)*m_dy <= b.ymin)) {
174	429x	num_rows--;
175		}
176	400531x	if (num_cols > 2 && (snapped_xmin + (num_cols-1)*m_dx >= b.xmax)) {
177	423x	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	400531x	Box reduced_box = {
185		snapped_xmin,
186	400531x	std::min(snapped_ymax - num_rows * m_dy, b.ymin),
187	400531x	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	400531x	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	400531x	if (reduced_box.ymin < m_extent.ymin) {
202	!	if (std::round((reduced_box.ymax - reduced_box.ymin)/m_dy) ==
203	!	std::round((reduced_box.ymax - m_extent.ymin)/m_dy)) {
204	!	reduced_box.ymin = m_extent.ymin;
205		} else {
206	!	throw std::runtime_error("Shrink operation failed.");
207		}
208		}
209
210	400531x	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	801062x	return reduced;
217		}
218
219		template<typename extent_tag2>
220	103x	bool compatible_with(const Grid<extent_tag2> &b) 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		// Perhaps it's possible to remove this hardcoded tolerance and express something in terms of
226		// std::numeric_limits<double>::epsilon(), but I wasn't able to make that work.
227	103x	constexpr double compatability_tol = 1e-6;
228
229		if (empty() \|\| b.empty()) {
230	!	return true;
231		}
232
233		// Check x-resolution compatibility
234	103x	if (!is_integral(std::max(m_dx, b.m_dx) / std::min(m_dx, b.m_dx), std::min(m_dx, b.m_dx)*compatability_tol)) {
235	!	return false;
236		}
237
238		// Check y-resolution compatibility
239	103x	if (!is_integral(std::max(m_dy, b.m_dy) / std::min(m_dy, b.m_dy), std::min(m_dy, b.m_dy)*compatability_tol)) {
240	!	return false;
241		}
242
243		// Check left-hand boundary compatibility
244	103x	if (!is_integral(std::abs(b.m_extent.xmin - m_extent.xmin) / std::min(m_dx, b.m_dx), std::min(m_dx, b.m_dx)*compatability_tol)) {
245	1x	return false;
246		}
247
248		// Check upper boundary compatibility
249	102x	if (!is_integral(std::abs(b.m_extent.ymax - m_extent.ymax) / std::min(m_dy, b.m_dy), std::min(m_dy, b.m_dy)*compatability_tol)) {
250	!	return false;
251		}
252
253	102x	return true;
254		}
255
256		template<typename extent_tag2>
257	78x	Grid<extent_tag> common_grid(const Grid<extent_tag2> &b) const {
258	78x	if (!compatible_with(b)) {
259	1x	throw std::runtime_error("Incompatible extents.");
260		}
261
262	77x	if (b.empty()) {
263	!	return *this;
264		}
265
266	77x	const double common_dx = std::min(m_dx, b.m_dx);
267	77x	const double common_dy = std::min(m_dy, b.m_dy);
268
269	77x	const double common_xmin = std::min(m_extent.xmin, b.m_extent.xmin);
270	77x	const double common_ymax = std::max(m_extent.ymax, b.m_extent.ymax);
271
272	77x	double common_xmax = std::max(m_extent.xmax, b.m_extent.xmax);
273	77x	double common_ymin = std::min(m_extent.ymin, b.m_extent.ymin);
274
275	77x	const long nx = static_cast<long>(std::round((common_xmax - common_xmin) / common_dx));
276	77x	const long ny = static_cast<long>(std::round((common_ymax - common_ymin) / common_dy));
277
278	77x	common_xmax = std::max(common_xmax, common_xmin + nx*common_dx);
279	77x	common_ymin = std::min(common_ymin, common_ymax - ny*common_dy);
280
281	77x	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) const {
286		if (!compatible_with(b)) {
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	399180x	bool operator==(const Grid<extent_tag> &b) const {
313		return
314	399180x	m_extent == b.m_extent &&
315	399406x	m_dx == b.m_dx &&
316	399406x	m_dy == b.m_dy;
317		}
318
319	399180x	bool operator!=(const Grid<extent_tag> &b) const {
320	399180x	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-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_H
15		#define EXACTEXTRACT_RASTER_H
16
17		#include <array>
18		#include <cmath>
19		#include <limits>
20
21		#include "grid.h"
22		#include "matrix.h"
23
24		namespace exactextract {
25		template<typename T>
26		class AbstractRaster {
27		public:
28		using value_type = T;
29
30	916982x	explicit AbstractRaster(const Grid<bounded_extent> & ex) :
31		m_grid{ex},
32	916982x	m_nodata{std::is_floating_point<T>::value ? std::numeric_limits<T>::quiet_NaN() : std::numeric_limits<T>::min()},
33	916982x	m_has_nodata{false}
34		{}
35
36		AbstractRaster(const Grid<bounded_extent> & ex, const T& nodata_val) : m_grid{ex}, m_nodata{nodata_val}, m_has_nodata{true} {}
37
38	916982x	virtual ~AbstractRaster() = default;
39
40	3275154x	size_t rows() const {
41	3275154x	return m_grid.rows();
42		}
43
44	4899362x	size_t cols() const {
45	4899362x	return m_grid.cols();
46		}
47
48		size_t size() const {
49		return rows() * cols();
50		}
51
52	398975x	double xres() const {
53	398975x	return m_grid.dx();
54		}
55
56	398975x	double yres() const {
57	398975x	return m_grid.dy();
58		}
59
60	398975x	double xmin() const {
61	398975x	return m_grid.xmin();
62		}
63
64		double ymin() const {
65		return m_grid.ymin();
66		}
67
68		double xmax() const {
69		return m_grid.xmax();
70		}
71
72	398975x	double ymax() const {
73	398975x	return m_grid.ymax();
74		}
75
76	5790452x	const Grid<bounded_extent>& grid() const {
77	5790452x	return m_grid;
78		}
79
80		virtual T operator()(size_t row, size_t col) const = 0;
81
82	398987x	bool has_nodata() const { return m_has_nodata; }
83
84	139706x	T nodata() const { return m_nodata; }
85
86	639963x	bool get(size_t row, size_t col, T & val) const {
87	639963x	val = operator()(row, col);
88
89	!	if (m_has_nodata && val == m_nodata) {
90	!	return false;
91		}
92	639963x	if (std::is_floating_point<T>::value && std::isnan(val)) {
93	23x	return false;
94		}
95
96	639940x	return true;
97		}
98
99	!	void set_nodata(const T & val) {
100	!	m_has_nodata = true;
101	!	m_nodata = val;
102		}
103
104		bool operator==(const AbstractRaster<T> & other) const {
105		if (rows() != other.rows())
106		return false;
107
108		if (cols() != other.cols())
109		return false;
110
111		if (xres() != other.xres())
112		return false;
113
114		if (yres() != other.yres())
115		return false;
116
117		if (xmin() != other.xmin())
118		return false;
119
120		if (ymin() != other.ymin())
121		return false;
122
123		// Do the rasters differ in their definition of NODATA? If so, we need to do some more detailed
124		// checking below.
125		bool nodata_differs = has_nodata() != other.has_nodata() \|\|
126		(has_nodata() && other.has_nodata() && nodata() != other.nodata());
127
128		for (size_t i = 0; i < rows(); i++) {
129		for (size_t j = 0; j < cols(); j++) {
130		if(operator()(i, j) != other(i, j)) {
131		// Override default behavior of NAN != NAN
132		if (!std::isnan(operator()(i, j)) \|\| !std::isnan(other(i, j))) {
133		return false;
134		}
135		} else if (nodata_differs && (operator()(i, j) == nodata() \|\| other(i, j) == other.nodata())) {
136		// For data types that do not have NAN, or even for floating point types where the user has
137		// selected some other value to represent NODATA, we need to reverse a positive equality test
138		// where the value is considered to be NODATA in one raster but not the other.
139		return false;
140		}
141		}
142		}
143
144		return true;
145		}
146
147		bool operator!=(const AbstractRaster<T> & other) const {
148		return !(operator==(other));
149		}
150
151		class Iterator : public std::iterator<std::forward_iterator_tag, T> {
152		public:
153		Iterator(const AbstractRaster<T>* r, size_t i, size_t j) :
154		m_rast(r), m_row(i), m_col(j) {}
155
156		const T& operator*() const {
157		m_val = m_rast->operator()(m_row, m_col);
158		return m_val;
159		}
160
161		Iterator& operator++() {
162		m_col++;
163		if (m_col == m_rast->cols()) {
164		m_col = 0;
165		m_row++;
166		}
167		return *this;
168		}
169
170		friend bool operator==(const Iterator& a,
171		const Iterator& b) {
172		return a.m_rast == b.m_rast && a.m_row == b.m_row && a.m_col == b.m_col;
173		}
174
175		friend bool operator!=(const Iterator& a,
176		const Iterator& b) {
177		return a.m_row != b.m_row \|\| a.m_col != b.m_col \|\| a.m_rast != b.m_rast;
178		}
179
180		private:
181		const AbstractRaster<T>* m_rast;
182		mutable T m_val;
183		size_t m_row;
184		size_t m_col;
185		};
186
187		Iterator begin() const {
188		return Iterator(this, 0, 0);
189		}
190
191		Iterator end() const {
192		return Iterator(this, rows(), 0);
193		}
194		private:
195		Grid<bounded_extent> m_grid;
196		T m_nodata;
197		bool m_has_nodata;
198		};
199
200		template<typename T>
201		class Raster : public AbstractRaster<T> {
202		public:
203		Raster(Matrix<T>&& values, const Box & box) : AbstractRaster<T>(
204		Grid<bounded_extent>(
205		box,
206		(box.xmax - box.xmin) / values.cols(),
207		(box.ymax - box.ymin) / values.rows())),
208		m_values{std::move(values)} {}
209
210	516886x	Raster(Matrix<T>&& values, const Grid<bounded_extent> & g) :
211		AbstractRaster<T>(g),
212	516886x	m_values{std::move(values)} {}
213
214		Raster(const Box & box, size_t nrow, size_t ncol) :
215		AbstractRaster<T>(Grid<bounded_extent>(box, (box.xmax-box.xmin) / ncol, (box.ymax-box.ymin) / nrow)),
216		m_values{nrow, ncol}
217		{}
218
219		explicit Raster(const Grid<bounded_extent> & ex) :
220		AbstractRaster<T>(ex),
221		m_values{ex.rows(), ex.cols()}
222		{}
223
224		Matrix<T>& data() {
225		return m_values;
226		}
227
228	16x	T& operator()(size_t row, size_t col) {
229	16x	return m_values(row, col);
230		}
231
232	656707x	T operator()(size_t row, size_t col) const override {
233	656707x	return m_values(row, col);
234		}
235
236		private:
237		Matrix<T> m_values;
238		};
239
240		template<typename T>
241		class RasterView : public AbstractRaster<T>{
242		public:
243		// Construct a view of a raster r at an extent ex that is larger
244		// and/or of finer resolution than r
245	398987x	RasterView(const AbstractRaster<T> & r, Grid<bounded_extent> ex) :
246		AbstractRaster<T>(ex),
247		m_raster{r},
248		m_x_off{0},
249		m_y_off{0},
250		m_rx{1},
251	398987x	m_ry{1} {
252	398987x	if (!this->grid().empty()) {
253	398975x	double disaggregation_factor_x = r.xres() / ex.dx();
254	398975x	double disaggregation_factor_y = r.yres() / ex.dy();
255
256	!	if (std::abs(disaggregation_factor_x - std::floor(disaggregation_factor_x)) > 1e-6 \|\|
257	398975x	std::abs(disaggregation_factor_y - std::floor(disaggregation_factor_y)) > 1e-6) {
258	!	throw std::runtime_error("Must construct view at resolution that is an integer multiple of original.");
259		}
260
261	398975x	if (disaggregation_factor_x < 0 \|\| disaggregation_factor_y < 0) {
262	!	throw std::runtime_error("Must construct view at equal or higher resolution than original.");
263		}
264
265	398975x	m_x_off = static_cast<long>(std::round((ex.xmin() - r.xmin()) / ex.dx()));
266	398975x	m_y_off = static_cast<long>(std::round((r.ymax() - ex.ymax()) / ex.dy()));
267	398975x	m_rx = static_cast<size_t>(disaggregation_factor_x);
268	398975x	m_ry = static_cast<size_t>(disaggregation_factor_y);
269		}
270
271	398987x	if (r.has_nodata()) {
272	!	this->set_nodata(r.nodata());
273		}
274		}
275
276	782104x	T operator()(size_t row, size_t col) const override {
277	782104x	if (m_raster.grid().empty()) {
278	8x	return this->nodata();
279		}
280
281	158401x	if (m_x_off < 0 && static_cast<size_t>(-m_x_off) > col) {
282	33718x	return this->nodata();
283		}
284	124683x	if (m_y_off < 0 && static_cast<size_t>(-m_y_off) > row) {
285	19136x	return this->nodata();
286		}
287
288	729242x	size_t i0 = (row + m_y_off) / m_ry;
289	729242x	size_t j0 = (col + m_x_off) / m_rx;
290
291	729242x	if (i0 > m_raster.rows() - 1 \|\| j0 > m_raster.cols() - 1) {
292	86844x	return this->nodata();
293		}
294
295	642398x	return m_raster(i0, j0);
296		}
297
298		private:
299		const AbstractRaster<T>& m_raster;
300
301		long m_x_off;
302		long m_y_off;
303		size_t m_rx;
304		size_t m_ry;
305		};
306
307
308		template<typename T>
309		std::ostream& operator<<(std::ostream & os, const AbstractRaster<T> & m) {
310		for (size_t i = 0; i < m.rows(); i++) {
311		for (size_t j = 0; j < m.cols(); j++) {
312		if (m(i, j) != 0) {
313		os << std::right << std::fixed << std::setw(10) << std::setprecision(6) <<
314		m(i, j) << " ";
315		} else {
316		os << " ";
317		}
318		}
319		os << std::endl;
320		}
321
322		return os;
323		}
324		}
325
326		#endif //EXACTEXTRACT_RASTER_H

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	8726129x	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	8726129x	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	336x	static Box make_empty() {
47	336x	return {0, 0, 0, 0};
48		}
49
50	6377052x	double width() const {
51	6377052x	return xmax - xmin;
52		}
53
54	6377052x	double height() const {
55	6377052x	return ymax - ymin;
56		}
57
58	6369286x	double area() const {
59	6369286x	return width() * height();
60		}
61
62	3883x	double perimeter() const {
63	3883x	return 2 * width() + 2 * height();
64		}
65
66	401176x	bool intersects(const Box & other) const {
67	401176x	if (other.ymin > ymax)
68	42x	return false;
69	401134x	if (other.ymax < ymin)
70	50x	return false;
71	401084x	if (other.xmin > xmax)
72	!	return false;
73	401084x	if (other.xmax < xmin)
74	!	return false;
75
76	401084x	return true;
77		}
78
79	4100088x	Box intersection(const Box & other) const {
80		return {
81	4100088x	std::max(xmin, other.xmin),
82	4100088x	std::max(ymin, other.ymin),
83	4100088x	std::min(xmax, other.xmax),
84	4100088x	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	517232x	bool empty() const {
118	517232x	return xmin >= xmax \|\| ymin >= ymax;
119		}
120
121	5x	Box expand_to_include(const Box & other) const {
122	5x	if (empty()) {
123	!	return other;
124		}
125
126	5x	if (other.empty()) {
127	!	return *this;
128		}
129
130	5x	return { std::min(xmin, other.xmin),
131	5x	std::min(ymin, other.ymin),
132	5x	std::max(xmax, other.xmax),
133	5x	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	400610x	bool operator==(const Box& other) const {
143	400610x	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	13x	explicit CartesianAreaRaster(const Grid<bounded_extent> & ex) :
24		AbstractRaster<T>(ex),
25	13x	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	8x	explicit SphericalAreaRaster(const Grid<bounded_extent> &ex) :
41		AbstractRaster<T>(ex),
42	8x	m_areas(ex.rows())
43		{
44	8x	const auto& g = AbstractRaster<T>::grid();
45
46	8x	double dlon = g.dx();
47	8x	double dlat = g.dy();
48
49	24x	for (size_t i = 0; i < g.rows(); i++) {
50	16x	double y = g.y_for_row(i);
51	16x	double ymin = y - 0.5 * dlat;
52	16x	double ymax = y + 0.5 * dlat;
53
54	16x	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	32x	T operator()(size_t row, size_t column) const override {
59		(void) column;
60	32x	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

1		// Copyright (c) 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_VARIANCE_H
15		#define EXACTEXTRACT_VARIANCE_H
16
17		#include <cmath>
18
19		namespace exactextract {
20
21		class WestVariance {
22		/** \brief Implements an incremental algorithm for weighted standard
23		* deviation, variance, and coefficient of variation, as described in
24		* formula WV2 of West, D.H.D. (1979) "Updating Mean and Variance
25		* Estimates: An Improved Method". Communications of the ACM 22(9).
26		*/
27
28		private:
29		double sum_w = 0;
30		double mean = 0;