exactextractr coverage - 88.74%

Files
Source

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

if (!isGeneric("exact_extract")) {
	setGeneric("exact_extract", function(x, y, ...)
		standardGeneric("exact_extract"))
}

#' Extract or summarize values from Raster* objects
#'
#' Extracts the values of cells in a Raster* that are covered by a
#' simple feature collection containing polygonal geometries, as well as the
#' fraction of each cell that is covered by the polygon. Returns either
#' the result of a summary operation or function applied to the values
#' and coverage fractions (if \code{fun} is specified), or a data frame
#' containing the values and coverage fractions themselves (if \code{fun}
#' is \code{NULL}.)
#'
#' The value of \code{fun} may be set to a string (or vector of strings)
#' representing summary operations supported by the exactextract library.
#' If the input raster has a single layer and a single summary operation
#' is specified, \code{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,
#' \code{exact_extract} will return a data frame with a row for each feature
#' and a column for each summary operation / layer combination. (The
#' \code{force_df} can be used to always return a data frame instead of a vector.)
#' In all of the summary operations, \code{NA} values in the raster are ignored
#' (i.e., \code{na.rm = TRUE}.)
#'
#' The following summary operations are supported:
#'
#' \itemize{
#'  \item{\code{min} - the minimum defined value in any raster cell wholly or
#'                     partially covered by the polygon}
#'  \item{\code{max} - the maximum defined value in any raster cell wholly or
#'                     partially covered by the polygon}
#'  \item{\code{count} - the sum of fractions of raster cells with defined values
#'                       covered by the polygon}
#'  \item{\code{sum}   - the sum of defined raster cell values, multiplied by
#'                       the fraction of the cell that is covered by the polygon}
#'  \item{\code{mean} - the mean cell value, weighted by the fraction of each cell
#'                      that is covered by the polygon}
#'  \item{\code{median} - the median cell value, weighted by the fraction of each
#'                        cell that is covered by the polygon}
#'  \item{\code{quantile} - arbitrary quantile(s) of cell values, specified in
#'                          \code{quantiles}, weighted by the fraction of each
#'                          cell that is covered by the polygon}
#'  \item{\code{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.}
#'  \item{\code{majority} - synonym for \code{mode}}
#'  \item{\code{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.}
#'  \item{\code{variety} - the number of distinct values in cells that are wholly
#'                         or partially covered by the polygon.}
#'  \item{\code{variance} - the population variance of cell values, weighted by the
#'                          fraction of each cell that is covered by the polygon.}
#'  \item{\code{stdev} - the population standard deviation of cell values, weighted
#'                       by the fraction of each cell that is covered by the polygon.}
#'  \item{\code{coefficient_of_variation} - the population coefficient of variation of
#'                       cell values, weighted by the fraction of each cell that is
#'                       covered by the polygon.}
#'  \item{\code{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 \code{weights}}
#'  \item{\code{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 \code{weights}}
#' }
#'
#' Alternatively, an R function may be provided as \code{fun}. The function will be
#' called for each feature with with vectors of cell values and weights as arguments.
#' \code{exact_extract} will then return a vector of the return values of \code{fun}.
#'
#' If \code{fun} is not specified, \code{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 values from each layer in the input `Raster*`,
#' and a final column indicating the fraction of the cell that is covered by the
#' polygon.
#'
#' @param     x a \code{RasterLayer}, \code{RasterStack}, or \code{RasterBrick}
#' @param     y a sf 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 \code{weighted_mean}
#'                     and \code{weighted_sum} summary operations.
#' @param     quantiles   quantiles to be computed when \code{fun == 'quantile'}
#' @param     append_cols when \code{fun} is not \code{NULL}, an optional
#'                        character vector of columns from \code{y} to be
#'                        included in returned data frame.
#' @param     force_df always return a data frame instead of a vector, even if
#'                     \code{x} has only one layer and \code{fun} has length 1
#' @param     full_colnames include the names of \code{x} in the names of the
#'                          returned data frame, even if \code{x} has only one
#'                          layer. This is useful when the results of multiple
#'                          calls to \code{exact_extract} are combined with
#'                          \code{cbind}.
#' @param     include_cell if \code{TRUE}, and \code{fun} is \code{NULL}, augment
#'                       the returned data frame for each feature with a column
#'                       for the cell index (\code{cell}). If \code{TRUE} and
#'                       \code{fun} is not \code{NULL}, add \code{cell} to the
#'                       data frame passed to \code{fun} for each feature.
#' @param     include_cols an optional character vector of column names in
#'                         \code{y} to be added to the data frame for each
#'                         feature that is either returned (when \code{fun} is
#'                         \code{NULL}) or passed to \code{fun}.
#' @param     include_xy if \code{TRUE}, and \code{fun} is \code{NULL}, augment
#'                       the returned data frame for each feature with columns
#'                       for cell center coordinates (\code{x} and \code{y}). If
#'                       \code{TRUE} and \code{fun} is not \code{NULL}, add
#'                       \code{x} and {y} to the data frame passed to \code{fun}
#'                       for each feature.
#' @param     stack_apply   if \code{TRUE}, apply \code{fun} to each layer of
#'                          \code{x} independently. If \code{FALSE}, apply \code{fun}
#'                          to all layers of \code{x} 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 \code{fun} (as opposed to a function defined using
#'                                R code). If a polygon covers more than \code{max_cells_in_memory}
#'                                raster cells, it will be processed in multiple chunks.
#' @param     progress if \code{TRUE}, display a progress bar during processing
#' @param     ... additional arguments to pass to \code{fun}
#' @return a vector or list of data frames, depending on the type of \code{x} and the
#'         value of \code{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
NULL

#' @import sf
#' @import raster
#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='Raster', y='sf'),
          function(x, y, fun=NULL, ...,
                   include_xy=FALSE,
                   progress=TRUE,
                   max_cells_in_memory=30000000,
                   include_cell=FALSE,
                   force_df=FALSE,
                   full_colnames=FALSE,
                   stack_apply=FALSE,
                   append_cols=NULL,
                   include_cols=NULL,
                   quantiles=NULL) {
  .exact_extract(x, y, fun=fun, ...,
                include_xy=include_xy,
                progress=progress,
                max_cells_in_memory=max_cells_in_memory,
                include_cell=include_cell,
                force_df=force_df,
                full_colnames=full_colnames,
                stack_apply=stack_apply,
                append_cols=append_cols,
                include_cols=include_cols,
                quantiles=quantiles)
})

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

emptyVector <- function(rast) {
  switch(substr(raster::dataType(rast), 1, 3),
         LOG=logical(),
         INT=integer(),
         numeric())
}

.exact_extract <- function(x, y, fun=NULL, ...,
                           weights=NULL,
                           include_xy=FALSE,
                           progress=TRUE,
                           max_cells_in_memory=30000000,
                           include_cell=FALSE,
                           force_df=FALSE,
                           full_colnames=FALSE,
                           stack_apply=FALSE,
                           append_cols=NULL,
                           include_cols=NULL,
                           quantiles=NULL) {
  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(sf::st_geometry_type(y, by_geometry = FALSE) == 'GEOMETRY') {
    if (!all(sf::st_dimension(y) == 2)) {
      stop("Features in sfc_GEOMETRY must be polygonal")
    }
    y <- sf::st_cast(y, 'MULTIPOLYGON')
  }

  if(!is.null(weights)) {
    if (!startsWith(class(weights), 'Raster')) {
      stop("Weights must be a Raster object.")
    }

    if (!is.character(fun)) {
      stop("Weighting raster can only be used with named summary operations.")
    }

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

    if (!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.")
      }
    }
  }

  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, ")")
  }

  if (!is.null(fun) && !is.character(fun) && .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 (is.character(fun)) {
    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 (!is.null(include_cols)) {
      stop("include_cols not supported for named_summary operations (see argument append_cols)")
    }
  }

  if (progress && length(y) > 1) {
    n <- length(y)
    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 <- raster::readStart(x)
    if (!is.null(weights)) {
      weights <- raster::readStart(weights)
    }

    if (is.character(fun)) {
      # Compute all stats in C++
      results <- sapply(sf::st_as_binary(sf::st_geometry(y), EWKB=TRUE), function(wkb) {
        ret <- CPP_stats(x, weights, wkb, fun, max_cells_in_memory, quantiles)
        update_progress()
        return(ret)
      })

      if (length(fun) == 1 &&
          raster::nlayers(x) == 1 &&
          (is.null(quantiles) || length(quantiles) == 1) &&
          !force_df) {
        # Just return a vector of stat results
        return(as.vector(results))
      } else {
        # Return a data frame with a column for each stat
        colnames <- .cppStatColNames(x, 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)

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

        return(ret)
      }
    } else {
      geoms <- sf::st_geometry(y)

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

        ret <- CPP_exact_extract(x, wkb)

        if (length(ret$weights) > 0) {
          vals <- raster::getValuesBlock(x,
                                         row=ret$row,
                                         col=ret$col,
                                         nrow=nrow(ret$weights),
                                         ncol=ncol(ret$weights))

          if(is.matrix(vals)) {
            vals <- as.data.frame(vals)
          } else {
            vals <- data.frame(value=vals)
          }
        } else {
          # Polygon does not intersect raster.
          # Construct a zero-row data frame with correct column names/types.
          vals <- do.call(data.frame, lapply(seq_len(raster::nlayers(x)),
                                                     function(i) emptyVector(x[[i]])))
          if (raster::nlayers(x) == 1) {
            names(vals) <- 'value'
          } else {
            names(vals) <- names(x)
          }
        }

        if (include_xy) {
           vals <- .appendXY(vals, x, ret$row, nrow(ret$weights), ret$col, ncol(ret$weights))
        }

        if (include_cell) {
          vals <- .appendCell(vals, x, ret$row, nrow(ret$weights), ret$col, ncol(ret$weights))
        }

        if (!is.null(include_cols)) {
          # use vals as first argument to cbind, then rearrange names so that
          # include_cols come first
          vals <- cbind(vals,
                        sf::st_drop_geometry(y[feature_num, include_cols]),
                        row.names = NULL)[, c(include_cols, names(vals))]
        }

        cov_fracs <- as.vector(t(ret$weights))
        vals <- vals[cov_fracs > 0, , drop=FALSE]
        cov_fracs <- cov_fracs[cov_fracs > 0]

        update_progress()

        if (is.null(fun)) {
          vals$coverage_fraction <- cov_fracs
          return(vals)
        } else {
          if (ncol(vals) == 1) {
            # Only one layer, nothing appended (cells or XY)
            return(fun(vals[,1], cov_fracs, ...))
          } else {
            if (stack_apply) {
              # Pass each layer in stack to callback individually
              nlay <- raster::nlayers(x)
              appended_cols <- seq_len(ncol(vals))[-seq_len(nlay)]

              if (length(appended_cols) == 0) {
                result <- lapply(seq_len(nlay), function(z)
                  fun(vals[, z], cov_fracs, ...))
              } else {
                result <- lapply(seq_len(nlay), function(z)
                  fun(cbind(data.frame(value=vals[, z]),
                            vals[, c(appended_cols)]), cov_fracs, ...))
              }

              names(result) <- paste('fun', names(x), sep='.')
              return(do.call(data.frame, result))
            } else {
              # Pass all layers to callback, to be handled together
              return(fun(vals, cov_fracs, ...))
            }
          }
        }
      })

      if (!is.null(fun)) {
        if (class(ret[[1]]) == 'data.frame') {
          if (requireNamespace('dplyr', quietly = TRUE)) {
            ret <- dplyr::bind_rows(ret) # handle column name mismatches
          } else {
            ret <- do.call(rbind, ret)
          }
        } else {
          ret <- simplify2array(ret)

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

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

      return(ret)
    }
  }, finally={
    raster::readStop(x)
    if (!is.null(weights)) {
      raster::readStop(weights)
    }
  })
}

.appendXY <- function(vals_df, rast, first_row, nrow, first_col, ncol) {
  if (nrow(vals_df) == 0) {
    vals_df$x <- numeric()
    vals_df$y <- numeric()
  } else {
    x_coords <- raster::xFromCol(rast, col=seq(first_col, first_col + ncol - 1))
    y_coords <- raster::yFromRow(rast, row=seq(first_row, first_row + nrow - 1))

    vals_df$x <- rep(x_coords, times=nrow)
    vals_df$y <- rep(y_coords, each=ncol)
  }

  return(vals_df)
}

.appendCell <- function(vals_df, rast, first_row, nrow, first_col, ncol) {
  if (nrow(vals_df) == 0) {
    vals_df$cell <- numeric()
  } else {
    rows <- rep(seq(first_row, first_row + nrow - 1), each = ncol)
    cols <- rep(seq(first_col, first_col + ncol - 1), times = nrow)

    vals_df$cell <- raster::cellFromRowCol(rast, row=rows, col=cols)
  }

  return(vals_df)
}

.cppStatColNames <- function(rast, stat_names, full_colnames, quantiles) {
  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])
  }

  if (raster::nlayers(rast) > 1 || full_colnames) {
    z <- expand.grid(names(rast), stat_names, stringsAsFactors=TRUE)
    mapply(paste, z[[2]], z[[1]], MoreArgs=list(sep='.'))
  } else {
    stat_names
  }
}

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

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

if (!isGeneric("exact_resample")) {
	setGeneric("exact_resample", function(x, y, ...)
		standardGeneric("exact_resample"))
}

#' Resample a raster to a new grid
#'
#' @param x a \code{RasterLayer} to be resampled
#' @param y a \code{RasterLayer} 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}
#'
#' @name exact_resample
NULL

#' @import sf
#' @import raster
#' @useDynLib exactextractr
#' @rdname exact_resample
#' @export
setMethod('exact_resample',
          signature(x='RasterLayer', y='RasterLayer'),
          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.')
              }
            }

            x <- raster::readStart(x)
            tryCatch({
              CPP_resample(x, y, fun)
            }, finally={
              raster::readStop(x)
            })
          })

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

if (!isGeneric('coverage_fraction')) {
	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) 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) 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-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) {
            RasterView<T> rv{rast, intersection_percentages.grid()};

            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;

            RasterView<float> iv{intersection_percentages, common};
            RasterView<T> rv{rast,    common};
            RasterView<T> wv{weights, common};

            for (size_t i = 0; i < rv.rows(); i++) {
                for (size_t j = 0; j < rv.cols(); j++) {
                    float pct_cov = iv(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) 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_H
#define EXACTEXTRACT_RASTER_H

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

        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) {
            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));
        }
    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} {
            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_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) 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;
    double t = 0;

public:
    /** \brief Update variance estimate with another value
     *
     * @param x value to add
     * @param w weight of `x`
     */
    void process(double x, double w) {
        double mean_old = mean;

        sum_w += w;
        mean += (w / sum_w) * (x - mean_old);
        t += w * (x - mean_old) * (x - mean);
    }

    /** \brief Return the population variance.
     */
    constexpr double variance() const {
        return t / sum_w;
    }

    /** \brief Return the population standard deviation
     */
    double stdev() const {
        return std::sqrt(variance());
    }

    /** \brief Return the population coefficient of variation
     */
    double coefficent_of_variation() const {
        return stdev() / mean;
    }

};

}

#endif //EXACTEXTRACT_VARIANCE_H

// Copyright (c) 2019-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_WEIGHTED_QUANTILES_H
#define EXACTEXTRACT_WEIGHTED_QUANTILES_H

#include <algorithm>
#include <cmath>
#include <stdexcept>
#include <vector>

namespace exactextract {

    class WeightedQuantiles {
        // Compute weighted quantiles based on https://stats.stackexchange.com/a/13223

    public:

        void process(double x, double w) {
            if (w < 0) {
                throw std::runtime_error("Weighted quantile calculation does not support negative weights.");
            }

            if (!std::isfinite(w)) {
                throw std::runtime_error("Weighted quantile does not support non-finite weights.");
            }

            m_ready_to_query = false;

            m_elems.emplace_back(x, w);
        }

        double quantile(double q) const;

    private:
        struct elem_t {
            elem_t(double _x, double _w) : x(_x), w(_w), cumsum(0), s(0) {}

            double x;
            double w;
            double cumsum;
            double s;
        };

        void prepare() const;

        mutable std::vector<elem_t> m_elems;
        mutable double m_sum_w;
        mutable bool m_ready_to_query;
    };

}

#endif //EXACTEXTRACT_WEIGHTED_QUANTILES_H

//
// Copyright (c) 2014-2018 Martin Moene
//
// https://github.com/martinmoene/optional-lite
//
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#pragma once

#ifndef NONSTD_OPTIONAL_LITE_HPP
#define NONSTD_OPTIONAL_LITE_HPP

#define optional_lite_MAJOR  3
#define optional_lite_MINOR  2
#define optional_lite_PATCH  0

#define optional_lite_VERSION  optional_STRINGIFY(optional_lite_MAJOR) "." optional_STRINGIFY(optional_lite_MINOR) "." optional_STRINGIFY(optional_lite_PATCH)

#define optional_STRINGIFY(  x )  optional_STRINGIFY_( x )
#define optional_STRINGIFY_( x )  #x

// optional-lite configuration:

#define optional_OPTIONAL_DEFAULT  0
#define optional_OPTIONAL_NONSTD   1
#define optional_OPTIONAL_STD      2

#if !defined( optional_CONFIG_SELECT_OPTIONAL )
# define optional_CONFIG_SELECT_OPTIONAL  ( optional_HAVE_STD_OPTIONAL ? optional_OPTIONAL_STD : optional_OPTIONAL_NONSTD )
#endif

// Control presence of exception handling (try and auto discover):

#ifndef optional_CONFIG_NO_EXCEPTIONS
# if defined(__cpp_exceptions) || defined(__EXCEPTIONS) || defined(_CPPUNWIND)
#  define optional_CONFIG_NO_EXCEPTIONS  0
# else
#  define optional_CONFIG_NO_EXCEPTIONS  1
# endif
#endif

// C++ language version detection (C++20 is speculative):
// Note: VC14.0/1900 (VS2015) lacks too much from C++14.

#ifndef   optional_CPLUSPLUS
# if defined(_MSVC_LANG ) && !defined(__clang__)
#  define optional_CPLUSPLUS  (_MSC_VER == 1900 ? 201103L : _MSVC_LANG )
# else
#  define optional_CPLUSPLUS  __cplusplus
# endif
#endif

#define optional_CPP98_OR_GREATER  ( optional_CPLUSPLUS >= 199711L )
#define optional_CPP11_OR_GREATER  ( optional_CPLUSPLUS >= 201103L )
#define optional_CPP11_OR_GREATER_ ( optional_CPLUSPLUS >= 201103L )
#define optional_CPP14_OR_GREATER  ( optional_CPLUSPLUS >= 201402L )
#define optional_CPP17_OR_GREATER  ( optional_CPLUSPLUS >= 201703L )
#define optional_CPP20_OR_GREATER  ( optional_CPLUSPLUS >= 202000L )

// C++ language version (represent 98 as 3):

#define optional_CPLUSPLUS_V  ( optional_CPLUSPLUS / 100 - (optional_CPLUSPLUS > 200000 ? 2000 : 1994) )

// Use C++17 std::optional if available and requested:

#if optional_CPP17_OR_GREATER && defined(__has_include )
# if __has_include( <optional> )
#  define optional_HAVE_STD_OPTIONAL  1
# else
#  define optional_HAVE_STD_OPTIONAL  0
# endif
#else
# define  optional_HAVE_STD_OPTIONAL  0
#endif

#define optional_USES_STD_OPTIONAL  ( (optional_CONFIG_SELECT_OPTIONAL == optional_OPTIONAL_STD) || ((optional_CONFIG_SELECT_OPTIONAL == optional_OPTIONAL_DEFAULT) && optional_HAVE_STD_OPTIONAL) )

//
// in_place: code duplicated in any-lite, expected-lite, optional-lite, value-ptr-lite, variant-lite:
//

#ifndef nonstd_lite_HAVE_IN_PLACE_TYPES
#define nonstd_lite_HAVE_IN_PLACE_TYPES  1

// C++17 std::in_place in <utility>:

#if optional_CPP17_OR_GREATER

#include <utility>

namespace nonstd {

using std::in_place;
using std::in_place_type;
using std::in_place_index;
using std::in_place_t;
using std::in_place_type_t;
using std::in_place_index_t;

#define nonstd_lite_in_place_t(      T)  std::in_place_t
#define nonstd_lite_in_place_type_t( T)  std::in_place_type_t<T>
#define nonstd_lite_in_place_index_t(K)  std::in_place_index_t<K>

#define nonstd_lite_in_place(      T)    std::in_place_t{}
#define nonstd_lite_in_place_type( T)    std::in_place_type_t<T>{}
#define nonstd_lite_in_place_index(K)    std::in_place_index_t<K>{}

} // namespace nonstd

#else // optional_CPP17_OR_GREATER

#include <cstddef>

namespace nonstd {
namespace detail {

template< class T >
struct in_place_type_tag {};

template< std::size_t K >
struct in_place_index_tag {};

} // namespace detail

struct in_place_t {};

template< class T >
inline in_place_t in_place( detail::in_place_type_tag<T> /*unused*/ = detail::in_place_type_tag<T>() )
{
    return in_place_t();
}

template< std::size_t K >
inline in_place_t in_place( detail::in_place_index_tag<K> /*unused*/ = detail::in_place_index_tag<K>() )
{
    return in_place_t();
}

template< class T >
inline in_place_t in_place_type( detail::in_place_type_tag<T> /*unused*/ = detail::in_place_type_tag<T>() )
{
    return in_place_t();
}

template< std::size_t K >
inline in_place_t in_place_index( detail::in_place_index_tag<K> /*unused*/ = detail::in_place_index_tag<K>() )
{
    return in_place_t();
}

// mimic templated typedef:

#define nonstd_lite_in_place_t(      T)  nonstd::in_place_t(&)( nonstd::detail::in_place_type_tag<T>  )
#define nonstd_lite_in_place_type_t( T)  nonstd::in_place_t(&)( nonstd::detail::in_place_type_tag<T>  )
#define nonstd_lite_in_place_index_t(K)  nonstd::in_place_t(&)( nonstd::detail::in_place_index_tag<K> )

#define nonstd_lite_in_place(      T)    nonstd::in_place_type<T>
#define nonstd_lite_in_place_type( T)    nonstd::in_place_type<T>
#define nonstd_lite_in_place_index(K)    nonstd::in_place_index<K>

} // namespace nonstd

#endif // optional_CPP17_OR_GREATER
#endif // nonstd_lite_HAVE_IN_PLACE_TYPES

//
// Using std::optional:
//

#if optional_USES_STD_OPTIONAL

#include <optional>

namespace nonstd {

    using std::optional;
    using std::bad_optional_access;
    using std::hash;

    using std::nullopt;
    using std::nullopt_t;

    using std::operator==;
    using std::operator!=;
    using std::operator<;
    using std::operator<=;
    using std::operator>;
    using std::operator>=;
    using std::make_optional;
    using std::swap;
}

#else // optional_USES_STD_OPTIONAL

#include <cassert>
#include <utility>

// optional-lite alignment configuration:

#ifndef  optional_CONFIG_MAX_ALIGN_HACK
# define optional_CONFIG_MAX_ALIGN_HACK  0
#endif

#ifndef  optional_CONFIG_ALIGN_AS
// no default, used in #if defined()
#endif

#ifndef  optional_CONFIG_ALIGN_AS_FALLBACK
# define optional_CONFIG_ALIGN_AS_FALLBACK  double
#endif

// Compiler warning suppression:

#if defined(__clang__)
# pragma clang diagnostic push
# pragma clang diagnostic ignored "-Wundef"
#elif defined(__GNUC__)
# pragma GCC   diagnostic push
# pragma GCC   diagnostic ignored "-Wundef"
#elif defined(_MSC_VER )
# pragma warning( push )
#endif

// half-open range [lo..hi):
#define optional_BETWEEN( v, lo, hi ) ( (lo) <= (v) && (v) < (hi) )

// Compiler versions:
//
// MSVC++ 6.0  _MSC_VER == 1200 (Visual Studio 6.0)
// MSVC++ 7.0  _MSC_VER == 1300 (Visual Studio .NET 2002)
// MSVC++ 7.1  _MSC_VER == 1310 (Visual Studio .NET 2003)
// MSVC++ 8.0  _MSC_VER == 1400 (Visual Studio 2005)
// MSVC++ 9.0  _MSC_VER == 1500 (Visual Studio 2008)
// MSVC++ 10.0 _MSC_VER == 1600 (Visual Studio 2010)
// MSVC++ 11.0 _MSC_VER == 1700 (Visual Studio 2012)
// MSVC++ 12.0 _MSC_VER == 1800 (Visual Studio 2013)
// MSVC++ 14.0 _MSC_VER == 1900 (Visual Studio 2015)
// MSVC++ 14.1 _MSC_VER >= 1910 (Visual Studio 2017)

#if defined(_MSC_VER ) && !defined(__clang__)
# define optional_COMPILER_MSVC_VER      (_MSC_VER )
# define optional_COMPILER_MSVC_VERSION  (_MSC_VER / 10 - 10 * ( 5 + (_MSC_VER < 1900 ) ) )
#else
# define optional_COMPILER_MSVC_VER      0
# define optional_COMPILER_MSVC_VERSION  0
#endif

#define optional_COMPILER_VERSION( major, minor, patch )  ( 10 * (10 * (major) + (minor) ) + (patch) )

#if defined(__GNUC__) && !defined(__clang__)
# define optional_COMPILER_GNUC_VERSION   optional_COMPILER_VERSION(__GNUC__, __GNUC_MINOR__, __GNUC_PATCHLEVEL__)
#else
# define optional_COMPILER_GNUC_VERSION   0
#endif

#if defined(__clang__)
# define optional_COMPILER_CLANG_VERSION  optional_COMPILER_VERSION(__clang_major__, __clang_minor__, __clang_patchlevel__)
#else
# define optional_COMPILER_CLANG_VERSION  0
#endif

#if optional_BETWEEN(optional_COMPILER_MSVC_VERSION, 70, 140 )
# pragma warning( disable: 4345 )   // initialization behavior changed
#endif

#if optional_BETWEEN(optional_COMPILER_MSVC_VERSION, 70, 150 )
# pragma warning( disable: 4814 )   // in C++14 'constexpr' will not imply 'const'
#endif

// Presence of language and library features:

#define optional_HAVE(FEATURE) ( optional_HAVE_##FEATURE )

#ifdef _HAS_CPP0X
# define optional_HAS_CPP0X  _HAS_CPP0X
#else
# define optional_HAS_CPP0X  0
#endif

// Unless defined otherwise below, consider VC14 as C++11 for optional-lite:

#if optional_COMPILER_MSVC_VER >= 1900
# undef  optional_CPP11_OR_GREATER
# define optional_CPP11_OR_GREATER  1
#endif

#define optional_CPP11_90   (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1500)
#define optional_CPP11_100  (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1600)
#define optional_CPP11_110  (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1700)
#define optional_CPP11_120  (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1800)
#define optional_CPP11_140  (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1900)
#define optional_CPP11_141  (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1910)

#define optional_CPP14_000  (optional_CPP14_OR_GREATER)
#define optional_CPP17_000  (optional_CPP17_OR_GREATER)

// Presence of C++11 language features:

#define optional_HAVE_CONSTEXPR_11      optional_CPP11_140
#define optional_HAVE_IS_DEFAULT        optional_CPP11_140
#define optional_HAVE_NOEXCEPT          optional_CPP11_140
#define optional_HAVE_NULLPTR           optional_CPP11_100
#define optional_HAVE_REF_QUALIFIER     optional_CPP11_140

// Presence of C++14 language features:

#define optional_HAVE_CONSTEXPR_14      optional_CPP14_000

// Presence of C++17 language features:

#define optional_HAVE_NODISCARD         optional_CPP17_000

// Presence of C++ library features:

#define optional_HAVE_CONDITIONAL       optional_CPP11_120
#define optional_HAVE_REMOVE_CV         optional_CPP11_120
#define optional_HAVE_TYPE_TRAITS       optional_CPP11_90

#define optional_HAVE_TR1_TYPE_TRAITS   (!! optional_COMPILER_GNUC_VERSION )
#define optional_HAVE_TR1_ADD_POINTER   (!! optional_COMPILER_GNUC_VERSION )

// C++ feature usage:

#if optional_HAVE( CONSTEXPR_11 )
# define optional_constexpr  constexpr
#else
# define optional_constexpr  /*constexpr*/
#endif

#if optional_HAVE( IS_DEFAULT )
# define optional_is_default  = default;
#else
# define optional_is_default  {}
#endif

#if optional_HAVE( CONSTEXPR_14 )
# define optional_constexpr14  constexpr
#else
# define optional_constexpr14  /*constexpr*/
#endif

#if optional_HAVE( NODISCARD )
# define optional_nodiscard  [[nodiscard]]
#else
# define optional_nodiscard  /*[[nodiscard]]*/
#endif

#if optional_HAVE( NOEXCEPT )
# define optional_noexcept  noexcept
#else
# define optional_noexcept  /*noexcept*/
#endif

#if optional_HAVE( NULLPTR )
# define optional_nullptr  nullptr
#else
# define optional_nullptr  NULL
#endif

#if optional_HAVE( REF_QUALIFIER )
// NOLINTNEXTLINE( bugprone-macro-parentheses )
# define optional_ref_qual  &
# define optional_refref_qual  &&
#else
# define optional_ref_qual  /*&*/
# define optional_refref_qual  /*&&*/
#endif

// additional includes:

#if optional_CONFIG_NO_EXCEPTIONS
// already included: <cassert>
#else
# include <stdexcept>
#endif

#if optional_CPP11_OR_GREATER
# include <functional>
#endif

#if optional_HAVE( INITIALIZER_LIST )
# include <initializer_list>
#endif

#if optional_HAVE( TYPE_TRAITS )
# include <type_traits>
#elif optional_HAVE( TR1_TYPE_TRAITS )
# include <tr1/type_traits>
#endif

// Method enabling

#if optional_CPP11_OR_GREATER

#define optional_REQUIRES_0(...) \
    template< bool B = (__VA_ARGS__), typename std::enable_if<B, int>::type = 0 >

#define optional_REQUIRES_T(...) \
    , typename = typename std::enable_if< (__VA_ARGS__), nonstd::optional_lite::detail::enabler >::type

#define optional_REQUIRES_R(R, ...) \
    typename std::enable_if< (__VA_ARGS__), R>::type

#define optional_REQUIRES_A(...) \
    , typename std::enable_if< (__VA_ARGS__), void*>::type = nullptr

#endif

//
// optional:
//

namespace nonstd { namespace optional_lite {

namespace std11 {

#if optional_CPP11_OR_GREATER
    using std::move;
#else
    template< typename T > T & move( T & t ) { return t; }
#endif

#if optional_HAVE( CONDITIONAL )
    using std::conditional;
#else
    template< bool B, typename T, typename F > struct conditional              { typedef T type; };
    template<         typename T, typename F > struct conditional<false, T, F> { typedef F type; };
#endif // optional_HAVE_CONDITIONAL

} // namespace std11

#if optional_CPP11_OR_GREATER

/// type traits C++17:

namespace std17 {

#if optional_CPP17_OR_GREATER

using std::is_swappable;
using std::is_nothrow_swappable;

#elif optional_CPP11_OR_GREATER

namespace detail {

using std::swap;

struct is_swappable
{
    template< typename T, typename = decltype( swap( std::declval<T&>(), std::declval<T&>() ) ) >
    static std::true_type test( int /*unused*/ );

    template< typename >
    static std::false_type test(...);
};

struct is_nothrow_swappable
{
    // wrap noexcept(expr) in separate function as work-around for VC140 (VS2015):

    template< typename T >
    static constexpr bool satisfies()
    {
        return noexcept( swap( std::declval<T&>(), std::declval<T&>() ) );
    }

    template< typename T >
    static auto test( int /*unused*/ ) -> std::integral_constant<bool, satisfies<T>()>{}

    template< typename >
    static auto test(...) -> std::false_type;
};

} // namespace detail

// is [nothow] swappable:

template< typename T >
struct is_swappable : decltype( detail::is_swappable::test<T>(0) ){};

template< typename T >
struct is_nothrow_swappable : decltype( detail::is_nothrow_swappable::test<T>(0) ){};

#endif // optional_CPP17_OR_GREATER

} // namespace std17

/// type traits C++20:

namespace std20 {

template< typename T >
struct remove_cvref
{
    typedef typename std::remove_cv< typename std::remove_reference<T>::type >::type type;
};

} // namespace std20

#endif // optional_CPP11_OR_GREATER

/// class optional

template< typename T >
class optional;

namespace detail {

// for optional_REQUIRES_T

#if optional_CPP11_OR_GREATER
enum class enabler{};
#endif

// C++11 emulation:

struct nulltype{};

template< typename Head, typename Tail >
struct typelist
{
    typedef Head head;
    typedef Tail tail;
};

#if optional_CONFIG_MAX_ALIGN_HACK

// Max align, use most restricted type for alignment:

#define optional_UNIQUE(  name )       optional_UNIQUE2( name, __LINE__ )
#define optional_UNIQUE2( name, line ) optional_UNIQUE3( name, line )
#define optional_UNIQUE3( name, line ) name ## line

#define optional_ALIGN_TYPE( type ) \
    type optional_UNIQUE( _t ); struct_t< type > optional_UNIQUE( _st )

template< typename T >
struct struct_t { T _; };

union max_align_t
{
    optional_ALIGN_TYPE( char );
    optional_ALIGN_TYPE( short int );
    optional_ALIGN_TYPE( int );
    optional_ALIGN_TYPE( long int  );
    optional_ALIGN_TYPE( float  );
    optional_ALIGN_TYPE( double );
    optional_ALIGN_TYPE( long double );
    optional_ALIGN_TYPE( char * );
    optional_ALIGN_TYPE( short int * );
    optional_ALIGN_TYPE( int *  );
    optional_ALIGN_TYPE( long int * );
    optional_ALIGN_TYPE( float * );
    optional_ALIGN_TYPE( double * );
    optional_ALIGN_TYPE( long double * );
    optional_ALIGN_TYPE( void * );

#ifdef HAVE_LONG_LONG
    optional_ALIGN_TYPE( long long );
#endif

    struct Unknown;

    Unknown ( * optional_UNIQUE(_) )( Unknown );
    Unknown * Unknown::* optional_UNIQUE(_);
    Unknown ( Unknown::* optional_UNIQUE(_) )( Unknown );

    struct_t< Unknown ( * )( Unknown)         > optional_UNIQUE(_);
    struct_t< Unknown * Unknown::*            > optional_UNIQUE(_);
    struct_t< Unknown ( Unknown::* )(Unknown) > optional_UNIQUE(_);
};

#undef optional_UNIQUE
#undef optional_UNIQUE2
#undef optional_UNIQUE3

#undef optional_ALIGN_TYPE

#elif defined( optional_CONFIG_ALIGN_AS ) // optional_CONFIG_MAX_ALIGN_HACK

// Use user-specified type for alignment:

#define optional_ALIGN_AS( unused ) \
    optional_CONFIG_ALIGN_AS

#else // optional_CONFIG_MAX_ALIGN_HACK

// Determine POD type to use for alignment:

#define optional_ALIGN_AS( to_align ) \
    typename type_of_size< alignment_types, alignment_of< to_align >::value >::type

template< typename T >
struct alignment_of;

template< typename T >
struct alignment_of_hack
{
    char c;
    T t;
    alignment_of_hack();
};

template< size_t A, size_t S >
struct alignment_logic
{
    enum { value = A < S ? A : S };
};

template< typename T >
struct alignment_of
{
    enum { value = alignment_logic<
        sizeof( alignment_of_hack<T> ) - sizeof(T), sizeof(T) >::value };
};

template< typename List, size_t N >
struct type_of_size
{
    typedef typename std11::conditional<
        N == sizeof( typename List::head ),
            typename List::head,
            typename type_of_size<typename List::tail, N >::type >::type type;
};

template< size_t N >
struct type_of_size< nulltype, N >
{
    typedef optional_CONFIG_ALIGN_AS_FALLBACK type;
};

template< typename T>
struct struct_t { T _; };

#define optional_ALIGN_TYPE( type ) \
    typelist< type , typelist< struct_t< type >

struct Unknown;

typedef
    optional_ALIGN_TYPE( char ),
    optional_ALIGN_TYPE( short ),
    optional_ALIGN_TYPE( int ),
    optional_ALIGN_TYPE( long ),
    optional_ALIGN_TYPE( float ),
    optional_ALIGN_TYPE( double ),
    optional_ALIGN_TYPE( long double ),

    optional_ALIGN_TYPE( char *),
    optional_ALIGN_TYPE( short * ),
    optional_ALIGN_TYPE( int * ),
    optional_ALIGN_TYPE( long * ),
    optional_ALIGN_TYPE( float * ),
    optional_ALIGN_TYPE( double * ),
    optional_ALIGN_TYPE( long double * ),

    optional_ALIGN_TYPE( Unknown ( * )( Unknown ) ),
    optional_ALIGN_TYPE( Unknown * Unknown::*     ),
    optional_ALIGN_TYPE( Unknown ( Unknown::* )( Unknown ) ),

    nulltype
    > > > > > > >    > > > > > > >
    > > > > > > >    > > > > > > >
    > > > > > >
    alignment_types;

#undef optional_ALIGN_TYPE

#endif // optional_CONFIG_MAX_ALIGN_HACK

/// C++03 constructed union to hold value.

template< typename T >
union storage_t
{
//private:
//    template< typename > friend class optional;

    typedef T value_type;

    storage_t() optional_is_default

    explicit storage_t( value_type const & v )
    {
        construct_value( v );
    }

    void construct_value( value_type const & v )
    {
        ::new( value_ptr() ) value_type( v );
    }

#if optional_CPP11_OR_GREATER

    explicit storage_t( value_type && v )
    {
        construct_value( std::move( v ) );
    }

    void construct_value( value_type && v )
    {
        ::new( value_ptr() ) value_type( std::move( v ) );
    }

    template< class... Args >
    void emplace( Args&&... args )
    {
        ::new( value_ptr() ) value_type( std::forward<Args>(args)... );
    }

    template< class U, class... Args >
    void emplace( std::initializer_list<U> il, Args&&... args )
    {
        ::new( value_ptr() ) value_type( il, std::forward<Args>(args)... );
    }

#endif

    void destruct_value()
    {
        value_ptr()->~T();
    }

    optional_nodiscard value_type const * value_ptr() const
    {
        return as<value_type>();
    }

    value_type * value_ptr()
    {
        return as<value_type>();
    }

    optional_nodiscard value_type const & value() const optional_ref_qual
    {
        return * value_ptr();
    }

    value_type & value() optional_ref_qual
    {
        return * value_ptr();
    }

#if optional_CPP11_OR_GREATER

    optional_nodiscard value_type const && value() const optional_refref_qual
    {
        return std::move( value() );
    }

    value_type && value() optional_refref_qual
    {
        return std::move( value() );
    }

#endif

#if optional_CPP11_OR_GREATER

    using aligned_storage_t = typename std::aligned_storage< sizeof(value_type), alignof(value_type) >::type;
    aligned_storage_t data;

#elif optional_CONFIG_MAX_ALIGN_HACK

    typedef struct { unsigned char data[ sizeof(value_type) ]; } aligned_storage_t;

    max_align_t hack;
    aligned_storage_t data;

#else
    typedef optional_ALIGN_AS(value_type) align_as_type;

    typedef struct { align_as_type data[ 1 + ( sizeof(value_type) - 1 ) / sizeof(align_as_type) ]; } aligned_storage_t;
    aligned_storage_t data;

#   undef optional_ALIGN_AS

#endif // optional_CONFIG_MAX_ALIGN_HACK

    optional_nodiscard void * ptr() optional_noexcept
    {
        return &data;
    }

    optional_nodiscard void const * ptr() const optional_noexcept
    {
        return &data;
    }

    template <typename U>
    optional_nodiscard U * as()
    {
        return reinterpret_cast<U*>( ptr() );
    }

    template <typename U>
    optional_nodiscard U const * as() const
    {
        return reinterpret_cast<U const *>( ptr() );
    }
};

} // namespace detail

/// disengaged state tag

struct nullopt_t
{
    struct init{};
    explicit optional_constexpr nullopt_t( init /*unused*/ ) optional_noexcept {}
};

#if optional_HAVE( CONSTEXPR_11 )
constexpr nullopt_t nullopt{ nullopt_t::init{} };
#else
// extra parenthesis to prevent the most vexing parse:
const nullopt_t nullopt(( nullopt_t::init() ));
#endif

/// optional access error

#if ! optional_CONFIG_NO_EXCEPTIONS

class bad_optional_access : public std::logic_error
{
public:
  explicit bad_optional_access()
  : logic_error( "bad optional access" ) {}
};

#endif //optional_CONFIG_NO_EXCEPTIONS

/// optional

template< typename T>
class optional
{
private:
    template< typename > friend class optional;

    typedef void (optional::*safe_bool)() const;

public:
    typedef T value_type;

    // x.x.3.1, constructors

    // 1a - default construct
    optional_constexpr optional() optional_noexcept
    : has_value_( false )
    , contained()
    {}

    // 1b - construct explicitly empty
    // NOLINTNEXTLINE( google-explicit-constructor, hicpp-explicit-conversions )
    optional_constexpr optional( nullopt_t /*unused*/ ) optional_noexcept
    : has_value_( false )
    , contained()
    {}

    // 2 - copy-construct
    optional_constexpr14 optional( optional const & other
#if optional_CPP11_OR_GREATER
        optional_REQUIRES_A(
            true || std::is_copy_constructible<T>::value
        )
#endif
    )
    : has_value_( other.has_value() )
    {
        if ( other.has_value() )
        {
            contained.construct_value( other.contained.value() );
        }
    }

#if optional_CPP11_OR_GREATER

    // 3 (C++11) - move-construct from optional
    optional_constexpr14 optional( optional && other
        optional_REQUIRES_A(
            true || std::is_move_constructible<T>::value
        )
        // NOLINTNEXTLINE( performance-noexcept-move-constructor )
    ) noexcept( std::is_nothrow_move_constructible<T>::value )
    : has_value_( other.has_value() )
    {
        if ( other.has_value() )
        {
            contained.construct_value( std::move( other.contained.value() ) );
        }
    }

    // 4a (C++11) - explicit converting copy-construct from optional
    template< typename U >
    explicit optional( optional<U> const & other
        optional_REQUIRES_A(
            std::is_constructible<T, U const &>::value
            && !std::is_constructible<T, optional<U> &          >::value
            && !std::is_constructible<T, optional<U> &&         >::value
            && !std::is_constructible<T, optional<U> const &    >::value
            && !std::is_constructible<T, optional<U> const &&   >::value
            && !std::is_convertible<     optional<U> &       , T>::value
            && !std::is_convertible<     optional<U> &&      , T>::value
            && !std::is_convertible<     optional<U> const & , T>::value
            && !std::is_convertible<     optional<U> const &&, T>::value
            && !std::is_convertible<               U const & , T>::value /*=> explicit */
        )
    )
    : has_value_( other.has_value() )
    {
        if ( other.has_value() )
        {
            contained.construct_value( T{ other.contained.value() } );
        }
    }
#endif // optional_CPP11_OR_GREATER

    // 4b (C++98 and later) - non-explicit converting copy-construct from optional
    template< typename U >
    // NOLINTNEXTLINE( google-explicit-constructor, hicpp-explicit-conversions )
    optional( optional<U> const & other
#if optional_CPP11_OR_GREATER
        optional_REQUIRES_A(
            std::is_constructible<T, U const &>::value
            && !std::is_constructible<T, optional<U> &          >::value
            && !std::is_constructible<T, optional<U> &&         >::value
            && !std::is_constructible<T, optional<U> const &    >::value
            && !std::is_constructible<T, optional<U> const &&   >::value
            && !std::is_convertible<     optional<U> &       , T>::value
            && !std::is_convertible<     optional<U> &&      , T>::value
            && !std::is_convertible<     optional<U> const & , T>::value
            && !std::is_convertible<     optional<U> const &&, T>::value
            &&  std::is_convertible<               U const & , T>::value /*=> non-explicit */
        )
#endif // optional_CPP11_OR_GREATER
    )
    : has_value_( other.has_value() )
    {
        if ( other.has_value() )
        {
            contained.construct_value( other.contained.value() );
        }
    }

#if optional_CPP11_OR_GREATER

    // 5a (C++11) - explicit converting move-construct from optional
    template< typename U >
    explicit optional( optional<U> && other
        optional_REQUIRES_A(
            std::is_constructible<T, U &&>::value
            && !std::is_constructible<T, optional<U> &          >::value
            && !std::is_constructible<T, optional<U> &&         >::value
            && !std::is_constructible<T, optional<U> const &    >::value
            && !std::is_constructible<T, optional<U> const &&   >::value
            && !std::is_convertible<     optional<U> &       , T>::value
            && !std::is_convertible<     optional<U> &&      , T>::value
            && !std::is_convertible<     optional<U> const & , T>::value
            && !std::is_convertible<     optional<U> const &&, T>::value
            && !std::is_convertible<                     U &&, T>::value /*=> explicit */
        )
    )
    : has_value_( other.has_value() )
    {
        if ( other.has_value() )
        {
            contained.construct_value( T{ std::move( other.contained.value() ) } );
        }
    }

    // 5a (C++11) - non-explicit converting move-construct from optional
    template< typename U >
    // NOLINTNEXTLINE( google-explicit-constructor, hicpp-explicit-conversions )
    optional( optional<U> && other
        optional_REQUIRES_A(
            std::is_constructible<T, U &&>::value
            && !std::is_constructible<T, optional<U> &          >::value
            && !std::is_constructible<T, optional<U> &&         >::value
            && !std::is_constructible<T, optional<U> const &    >::value
            && !std::is_constructible<T, optional<U> const &&   >::value
            && !std::is_convertible<     optional<U> &       , T>::value
            && !std::is_convertible<     optional<U> &&      , T>::value
            && !std::is_convertible<     optional<U> const & , T>::value
            && !std::is_convertible<     optional<U> const &&, T>::value
            &&  std::is_convertible<                     U &&, T>::value /*=> non-explicit */
        )
    )
    : has_value_( other.has_value() )
    {
        if ( other.has_value() )
        {
            contained.construct_value( std::move( other.contained.value() ) );
        }
    }

    // 6 (C++11) - in-place construct
    template< typename... Args
        optional_REQUIRES_T(
            std::is_constructible<T, Args&&...>::value
        )
    >
    optional_constexpr explicit optional( nonstd_lite_in_place_t(T), Args&&... args )
    : has_value_( true )
    , contained( T( std::forward<Args>(args)...) )
    {}

    // 7 (C++11) - in-place construct,  initializer-list
    template< typename U, typename... Args
        optional_REQUIRES_T(
            std::is_constructible<T, std::initializer_list<U>&, Args&&...>::value
        )
    >
    optional_constexpr explicit optional( nonstd_lite_in_place_t(T), std::initializer_list<U> il, Args&&... args )
    : has_value_( true )
    , contained( T( il, std::forward<Args>(args)...) )
    {}

    // 8a (C++11) - explicit move construct from value
    template< typename U = value_type >
    optional_constexpr explicit optional( U && value
        optional_REQUIRES_A(
            std::is_constructible<T, U&&>::value
            && !std::is_same<typename std20::remove_cvref<U>::type, nonstd_lite_in_place_t(U)>::value
            && !std::is_same<typename std20::remove_cvref<U>::type, optional<T>>::value
            && !std::is_convertible<U&&, T>::value /*=> explicit */
        )
    )
    : has_value_( true )
    , contained( T{ std::forward<U>( value ) } )
    {}

    // 8b (C++11) - non-explicit move construct from value
    template< typename U = value_type >
    // NOLINTNEXTLINE( google-explicit-constructor, hicpp-explicit-conversions )
    optional_constexpr optional( U && value
        optional_REQUIRES_A(
            std::is_constructible<T, U&&>::value
            && !std::is_same<typename std20::remove_cvref<U>::type, nonstd_lite_in_place_t(U)>::value
            && !std::is_same<typename std20::remove_cvref<U>::type, optional<T>>::value
            && std::is_convertible<U&&, T>::value /*=> non-explicit */
        )
    )
    : has_value_( true )
    , contained( std::forward<U>( value ) )
    {}

#else // optional_CPP11_OR_GREATER

    // 8 (C++98)
    optional( value_type const & value )
    : has_value_( true )
    , contained( value )
    {}

#endif // optional_CPP11_OR_GREATER

    // x.x.3.2, destructor

    ~optional()
    {
        if ( has_value() )
        {
            contained.destruct_value();
        }
    }

    // x.x.3.3, assignment

    // 1 (C++98and later) -  assign explicitly empty
    optional & operator=( nullopt_t /*unused*/) optional_noexcept
    {
        reset();
        return *this;
    }

    // 2 (C++98and later) - copy-assign from optional
#if optional_CPP11_OR_GREATER
    // NOLINTNEXTLINE( cppcoreguidelines-c-copy-assignment-signature, misc-unconventional-assign-operator )
    optional_REQUIRES_R(
        optional &,
        true
//      std::is_copy_constructible<T>::value
//      && std::is_copy_assignable<T>::value
    )
    operator=( optional const & other )
        noexcept(
            std::is_nothrow_move_assignable<T>::value
            && std::is_nothrow_move_constructible<T>::value
        )
#else
    optional & operator=( optional const & other )
#endif
    {
        if      ( (has_value() == true ) && (other.has_value() == false) ) { reset(); }
        else if ( (has_value() == false) && (other.has_value() == true ) ) { initialize( *other ); }
        else if ( (has_value() == true ) && (other.has_value() == true ) ) { contained.value() = *other; }
        return *this;
    }

#if optional_CPP11_OR_GREATER

    // 3 (C++11) - move-assign from optional
    // NOLINTNEXTLINE( cppcoreguidelines-c-copy-assignment-signature, misc-unconventional-assign-operator )
    optional_REQUIRES_R(
        optional &,
        true
//      std::is_move_constructible<T>::value
//      && std::is_move_assignable<T>::value
    )
    operator=( optional && other ) noexcept
    {
        if      ( (has_value() == true ) && (other.has_value() == false) ) { reset(); }
        else if ( (has_value() == false) && (other.has_value() == true ) ) { initialize( std::move( *other ) ); }
        else if ( (has_value() == true ) && (other.has_value() == true ) ) { contained.value() = std::move( *other ); }
        return *this;
    }

    // 4 (C++11) - move-assign from value
    template< typename U = T >
        // NOLINTNEXTLINE( cppcoreguidelines-c-copy-assignment-signature, misc-unconventional-assign-operator )
        optional_REQUIRES_R(
            optional &,
            std::is_constructible<T , U>::value
            && std::is_assignable<T&, U>::value
            && !std::is_same<typename std20::remove_cvref<U>::type, nonstd_lite_in_place_t(U)>::value
            && !std::is_same<typename std20::remove_cvref<U>::type, optional<T>>::value
            && !(std::is_scalar<T>::value && std::is_same<T, typename std::decay<U>::type>::value)
        )
    operator=( U && value )
    {
        if ( has_value() )
        {
            contained.value() = std::forward<U>( value );
        }
        else
        {
            initialize( T( std::forward<U>( value ) ) );
        }
        return *this;
    }

#else // optional_CPP11_OR_GREATER

    // 4 (C++98) - copy-assign from value
    template< typename U /*= T*/ >
    optional & operator=( U const & value )
    {
        if ( has_value() ) contained.value() = value;
        else               initialize( T( value ) );
        return *this;
    }

#endif // optional_CPP11_OR_GREATER

    // 5 (C++98 and later) - converting copy-assign from optional
    template< typename U >
#if optional_CPP11_OR_GREATER
        // NOLINTNEXTLINE( cppcoreguidelines-c-copy-assignment-signature, misc-unconventional-assign-operator )
        optional_REQUIRES_R(
            optional&,
            std::is_constructible<  T , U const &>::value
            &&  std::is_assignable< T&, U const &>::value
            && !std::is_constructible<T, optional<U> &          >::value
            && !std::is_constructible<T, optional<U> &&         >::value
            && !std::is_constructible<T, optional<U> const &    >::value
            && !std::is_constructible<T, optional<U> const &&   >::value
            && !std::is_convertible<     optional<U> &       , T>::value
            && !std::is_convertible<     optional<U> &&      , T>::value
            && !std::is_convertible<     optional<U> const & , T>::value
            && !std::is_convertible<     optional<U> const &&, T>::value
            && !std::is_assignable<  T&, optional<U> &          >::value
            && !std::is_assignable<  T&, optional<U> &&         >::value
            && !std::is_assignable<  T&, optional<U> const &    >::value
            && !std::is_assignable<  T&, optional<U> const &&   >::value
        )
#else
    optional&
#endif // optional_CPP11_OR_GREATER
    operator=( optional<U> const & other )
    {
        return *this = optional( other );
    }

#if optional_CPP11_OR_GREATER

    // 6 (C++11) -  converting move-assign from optional
    template< typename U >
        // NOLINTNEXTLINE( cppcoreguidelines-c-copy-assignment-signature, misc-unconventional-assign-operator )
        optional_REQUIRES_R(
            optional&,
            std::is_constructible<  T , U>::value
            &&  std::is_assignable< T&, U>::value
            && !std::is_constructible<T, optional<U> &          >::value
            && !std::is_constructible<T, optional<U> &&         >::value
            && !std::is_constructible<T, optional<U> const &    >::value
            && !std::is_constructible<T, optional<U> const &&   >::value
            && !std::is_convertible<     optional<U> &       , T>::value
            && !std::is_convertible<     optional<U> &&      , T>::value
            && !std::is_convertible<     optional<U> const & , T>::value
            && !std::is_convertible<     optional<U> const &&, T>::value
            && !std::is_assignable<  T&, optional<U> &          >::value
            && !std::is_assignable<  T&, optional<U> &&         >::value
            && !std::is_assignable<  T&, optional<U> const &    >::value
            && !std::is_assignable<  T&, optional<U> const &&   >::value
        )
    operator=( optional<U> && other )
    {
        return *this = optional( std::move( other ) );
    }

    // 7 (C++11) - emplace
    template< typename... Args
        optional_REQUIRES_T(
            std::is_constructible<T, Args&&...>::value
        )
    >
    T& emplace( Args&&... args )
    {
        *this = nullopt;
        contained.emplace( std::forward<Args>(args)...  );
        has_value_ = true;
        return contained.value();
    }

    // 8 (C++11) - emplace, initializer-list
    template< typename U, typename... Args
        optional_REQUIRES_T(
            std::is_constructible<T, std::initializer_list<U>&, Args&&...>::value
        )
    >
    T& emplace( std::initializer_list<U> il, Args&&... args )
    {
        *this = nullopt;
        contained.emplace( il, std::forward<Args>(args)...  );
        has_value_ = true;
        return contained.value();
    }

#endif // optional_CPP11_OR_GREATER

    // x.x.3.4, swap

    void swap( optional & other )
#if optional_CPP11_OR_GREATER
        noexcept(
            std::is_nothrow_move_constructible<T>::value
            && std17::is_nothrow_swappable<T>::value
        )
#endif
    {
        using std::swap;
        if      ( (has_value() == true ) && (other.has_value() == true ) ) { swap( **this, *other ); }
        else if ( (has_value() == false) && (other.has_value() == true ) ) { initialize( std11::move(*other) ); other.reset(); }
        else if ( (has_value() == true ) && (other.has_value() == false) ) { other.initialize( std11::move(**this) ); reset(); }
    }

    // x.x.3.5, observers

    optional_constexpr value_type const * operator ->() const
    {
        return assert( has_value() ),
            contained.value_ptr();
    }

    optional_constexpr14 value_type * operator ->()
    {
        return assert( has_value() ),
            contained.value_ptr();
    }

    optional_constexpr value_type const & operator *() const optional_ref_qual
    {
        return assert( has_value() ),
            contained.value();
    }

    optional_constexpr14 value_type & operator *() optional_ref_qual
    {
        return assert( has_value() ),
            contained.value();
    }

#if optional_HAVE( REF_QUALIFIER )  &&  ( !optional_COMPILER_GNUC_VERSION || optional_COMPILER_GNUC_VERSION >= 490 )

    optional_constexpr value_type const && operator *() const optional_refref_qual
    {
        return std::move( **this );
    }

    optional_constexpr14 value_type && operator *() optional_refref_qual
    {
        return std::move( **this );
    }

#endif

#if optional_CPP11_OR_GREATER
    optional_constexpr explicit operator bool() const optional_noexcept
    {
        return has_value();
    }
#else
    optional_constexpr operator safe_bool() const optional_noexcept
    {
        return has_value() ? &optional::this_type_does_not_support_comparisons : 0;
    }
#endif

    // NOLINTNEXTLINE( modernize-use-nodiscard )
    /*optional_nodiscard*/ optional_constexpr bool has_value() const optional_noexcept
    {
        return has_value_;
    }

    // NOLINTNEXTLINE( modernize-use-nodiscard )
    /*optional_nodiscard*/ optional_constexpr14 value_type const & value() const optional_ref_qual
    {
#if optional_CONFIG_NO_EXCEPTIONS
        assert( has_value() );
#else
        if ( ! has_value() )
        {
            throw bad_optional_access();
        }
#endif
        return contained.value();
    }

    optional_constexpr14 value_type & value() optional_ref_qual
    {
#if optional_CONFIG_NO_EXCEPTIONS
        assert( has_value() );
#else
        if ( ! has_value() )
        {
            throw bad_optional_access();
        }
#endif
        return contained.value();
    }

#if optional_HAVE( REF_QUALIFIER )  &&  ( !optional_COMPILER_GNUC_VERSION || optional_COMPILER_GNUC_VERSION >= 490 )

    // NOLINTNEXTLINE( modernize-use-nodiscard )
    /*optional_nodiscard*/ optional_constexpr value_type const && value() const optional_refref_qual
    {
        return std::move( value() );
    }

    optional_constexpr14 value_type && value() optional_refref_qual
    {
        return std::move( value() );
    }

#endif

#if optional_CPP11_OR_GREATER

    template< typename U >
    optional_constexpr value_type value_or( U && v ) const optional_ref_qual
    {
        return has_value() ? contained.value() : static_cast<T>(std::forward<U>( v ) );
    }

    template< typename U >
    optional_constexpr14 value_type value_or( U && v ) optional_refref_qual
    {
        return has_value() ? std::move( contained.value() ) : static_cast<T>(std::forward<U>( v ) );
    }

#else

    template< typename U >
    optional_constexpr value_type value_or( U const & v ) const
    {
        return has_value() ? contained.value() : static_cast<value_type>( v );
    }

#endif // optional_CPP11_OR_GREATER

    // x.x.3.6, modifiers

    void reset() optional_noexcept
    {
        if ( has_value() )
        {
            contained.destruct_value();
        }

        has_value_ = false;
    }

private:
    void this_type_does_not_support_comparisons() const {}

    template< typename V >
    void initialize( V const & value )
    {
        assert( ! has_value()  );
        contained.construct_value( value );
        has_value_ = true;
    }

#if optional_CPP11_OR_GREATER
    template< typename V >
    void initialize( V && value )
    {
        assert( ! has_value()  );
        contained.construct_value( std::move( value ) );
        has_value_ = true;
    }

#endif

private:
    bool has_value_;
    detail::storage_t< value_type > contained;

};

// Relational operators

template< typename T, typename U >
inline optional_constexpr bool operator==( optional<T> const & x, optional<U> const & y )
{
    return bool(x) != bool(y) ? false : !bool( x ) ? true : *x == *y;
}

template< typename T, typename U >
inline optional_constexpr bool operator!=( optional<T> const & x, optional<U> const & y )
{
    return !(x == y);
}

template< typename T, typename U >
inline optional_constexpr bool operator<( optional<T> const & x, optional<U> const & y )
{
    return (!y) ? false : (!x) ? true : *x < *y;
}

template< typename T, typename U >
inline optional_constexpr bool operator>( optional<T> const & x, optional<U> const & y )
{
    return (y < x);
}

template< typename T, typename U >
inline optional_constexpr bool operator<=( optional<T> const & x, optional<U> const & y )
{
    return !(y < x);
}

template< typename T, typename U >
inline optional_constexpr bool operator>=( optional<T> const & x, optional<U> const & y )
{
    return !(x < y);
}

// Comparison with nullopt

template< typename T >
inline optional_constexpr bool operator==( optional<T> const & x, nullopt_t /*unused*/ ) optional_noexcept
{
    return (!x);
}

template< typename T >
inline optional_constexpr bool operator==( nullopt_t /*unused*/, optional<T> const & x ) optional_noexcept
{
    return (!x);
}

template< typename T >
inline optional_constexpr bool operator!=( optional<T> const & x, nullopt_t /*unused*/ ) optional_noexcept
{
    return bool(x);
}

template< typename T >
inline optional_constexpr bool operator!=( nullopt_t /*unused*/, optional<T> const & x ) optional_noexcept
{
    return bool(x);
}

template< typename T >
inline optional_constexpr bool operator<( optional<T> const & /*unused*/, nullopt_t /*unused*/ ) optional_noexcept
{
    return false;
}

template< typename T >
inline optional_constexpr bool operator<( nullopt_t /*unused*/, optional<T> const & x ) optional_noexcept
{
    return bool(x);
}

template< typename T >
inline optional_constexpr bool operator<=( optional<T> const & x, nullopt_t /*unused*/ ) optional_noexcept
{
    return (!x);
}

template< typename T >
inline optional_constexpr bool operator<=( nullopt_t /*unused*/, optional<T> const & /*unused*/ ) optional_noexcept
{
    return true;
}

template< typename T >
inline optional_constexpr bool operator>( optional<T> const & x, nullopt_t /*unused*/ ) optional_noexcept
{
    return bool(x);
}

template< typename T >
inline optional_constexpr bool operator>( nullopt_t /*unused*/, optional<T> const & /*unused*/ ) optional_noexcept
{
    return false;
}

template< typename T >
inline optional_constexpr bool operator>=( optional<T> const & /*unused*/, nullopt_t /*unused*/ ) optional_noexcept
{
    return true;
}

template< typename T >
inline optional_constexpr bool operator>=( nullopt_t /*unused*/, optional<T> const & x ) optional_noexcept
{
    return (!x);
}

// Comparison with T

template< typename T, typename U >
inline optional_constexpr bool operator==( optional<T> const & x, U const & v )
{
    return bool(x) ? *x == v : false;
}

template< typename T, typename U >
inline optional_constexpr bool operator==( U const & v, optional<T> const & x )
{
    return bool(x) ? v == *x : false;
}

template< typename T, typename U >
inline optional_constexpr bool operator!=( optional<T> const & x, U const & v )
{
    return bool(x) ? *x != v : true;
}

template< typename T, typename U >
inline optional_constexpr bool operator!=( U const & v, optional<T> const & x )
{
    return bool(x) ? v != *x : true;
}

template< typename T, typename U >
inline optional_constexpr bool operator<( optional<T> const & x, U const & v )
{
    return bool(x) ? *x < v : true;
}

template< typename T, typename U >
inline optional_constexpr bool operator<( U const & v, optional<T> const & x )
{
    return bool(x) ? v < *x : false;
}

template< typename T, typename U >
inline optional_constexpr bool operator<=( optional<T> const & x, U const & v )
{
    return bool(x) ? *x <= v : true;
}

template< typename T, typename U >
inline optional_constexpr bool operator<=( U const & v, optional<T> const & x )
{
    return bool(x) ? v <= *x : false;
}

template< typename T, typename U >
inline optional_constexpr bool operator>( optional<T> const & x, U const & v )
{
    return bool(x) ? *x > v : false;
}

template< typename T, typename U >
inline optional_constexpr bool operator>( U const & v, optional<T> const & x )
{
    return bool(x) ? v > *x : true;
}

template< typename T, typename U >
inline optional_constexpr bool operator>=( optional<T> const & x, U const & v )
{
    return bool(x) ? *x >= v : false;
}

template< typename T, typename U >
inline optional_constexpr bool operator>=( U const & v, optional<T> const & x )
{
    return bool(x) ? v >= *x : true;
}

// Specialized algorithms

template< typename T
#if optional_CPP11_OR_GREATER
    optional_REQUIRES_T(
        std::is_move_constructible<T>::value
        && std17::is_swappable<T>::value )
#endif
>
void swap( optional<T> & x, optional<T> & y )
#if optional_CPP11_OR_GREATER
    noexcept( noexcept( x.swap(y) ) )
#endif
{
    x.swap( y );
}

#if optional_CPP11_OR_GREATER

template< typename T >
optional_constexpr optional< typename std::decay<T>::type > make_optional( T && value )
{
    return optional< typename std::decay<T>::type >( std::forward<T>( value ) );
}

template< typename T, typename...Args >
optional_constexpr optional<T> make_optional( Args&&... args )
{
    return optional<T>( nonstd_lite_in_place(T), std::forward<Args>(args)...);
}

template< typename T, typename U, typename... Args >
optional_constexpr optional<T> make_optional( std::initializer_list<U> il, Args&&... args )
{
    return optional<T>( nonstd_lite_in_place(T), il, std::forward<Args>(args)...);
}

#else

template< typename T >
optional<T> make_optional( T const & value )
{
    return optional<T>( value );
}

#endif // optional_CPP11_OR_GREATER

} // namespace optional_lite

using optional_lite::optional;
using optional_lite::nullopt_t;
using optional_lite::nullopt;
using optional_lite::bad_optional_access;

using optional_lite::make_optional;

} // namespace nonstd

#if optional_CPP11_OR_GREATER

// specialize the std::hash algorithm:

namespace std {

template< class T >
struct hash< nonstd::optional<T> >
{
public:
    std::size_t operator()( nonstd::optional<T> const & v ) const optional_noexcept
    {
        return bool( v ) ? std::hash<T>{}( *v ) : 0;
    }
};

} //namespace std

#endif // optional_CPP11_OR_GREATER

#if defined(__clang__)
# pragma clang diagnostic pop
#elif defined(__GNUC__)
# pragma GCC   diagnostic pop
#elif defined(_MSC_VER )
# pragma warning( pop )
#endif

#endif // optional_USES_STD_OPTIONAL

#endif // NONSTD_OPTIONAL_LITE_HPP

// Copyright (c) 2018 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 "area.h"

#include <cmath>
#include <cstddef>

namespace exactextract {

    double area_signed(const std::vector<Coordinate> &ring) {
        if (ring.size() < 3) {
            return 0;
        }

        double sum{0};

        double x0{ring[0].x};
        for (size_t i = 1; i < ring.size() - 1; i++) {
            double x = ring[i].x - x0;
            double y1 = ring[i + 1].y;
            double y2 = ring[i - 1].y;
            sum += x * (y2 - y1);
        }

        return sum / 2.0;
    }

    double area(const std::vector<Coordinate> &ring) {
        return std::abs(area_signed(ring));
    }

}

// Copyright (c) 2018 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 <cmath>
#include <stdexcept>

#include "box.h"

namespace exactextract {

    static inline double clamp(double x, double low, double high) {
        return std::min(std::max(x, low), high);
    }

    Side Box::side(const Coordinate &c) const {
        if (c.x == xmin) {
            return Side::LEFT;
        } else if (c.x == xmax) {
            return Side::RIGHT;
        } else if (c.y == ymin) {
            return Side::BOTTOM;
        } else if (c.y == ymax) {
            return Side::TOP;
        }

        return Side::NONE;
    }

    Crossing Box::crossing(const Coordinate &c1, const Coordinate &c2) const {
        // vertical line
        if (c1.x == c2.x) {
            if (c2.y >= ymax) {
                return Crossing{Side::TOP, c1.x, ymax};
            } else if (c2.y <= ymin) {
                return Crossing{Side::BOTTOM, c1.x, ymin};
            } else {
                throw std::runtime_error("Never get here.");
            }
        }

        // horizontal line
        if (c1.y == c2.y) {
            if (c2.x >= xmax) {
                return Crossing{Side::RIGHT, xmax, c1.y};
            } else if (c2.x <= xmin) {
                return Crossing{Side::LEFT, xmin, c1.y};
            } else {
                throw std::runtime_error("Never get here");
            }
        }

        double m = std::abs((c2.y - c1.y) / (c2.x - c1.x));

        bool up = c2.y > c1.y;
        bool right = c2.x > c1.x;

        if (up) {
            if (right) {
                // 1st quadrant
                double y2 = c1.y + m * (xmax - c1.x);

                if (y2 < ymax) {
                    return Crossing{Side::RIGHT, xmax, clamp(y2, ymin, ymax)};
                } else {
                    double x2 = c1.x + (ymax - c1.y) / m;
                    return Crossing{Side::TOP, clamp(x2, xmin, xmax), ymax};
                }
            } else {
                // 2nd quadrant
                double y2 = c1.y + m * (c1.x - xmin);

                if (y2 < ymax) {
                    return Crossing{Side::LEFT, xmin, clamp(y2, ymin, ymax)};
                } else {
                    double x2 = c1.x - (ymax - c1.y) / m;
                    return Crossing{Side::TOP, clamp(x2, xmin, xmax), ymax};
                }
            }
        } else {
            if (right) {
                // 4th quadrant
                double y2 = c1.y - m * (xmax - c1.x);

                if (y2 > ymin) {
                    return Crossing{Side::RIGHT, xmax, clamp(y2, ymin, ymax)};
                } else {
                    double x2 = c1.x + (c1.y - ymin) / m;
                    return Crossing{Side::BOTTOM, clamp(x2, xmin, xmax), ymin};
                }
            } else {
                // 3rd quadrant
                double y2 = c1.y - m * (c1.x - xmin);

                if (y2 > ymin) {
                    return Crossing{Side::LEFT, xmin, clamp(y2, ymin, ymax)};
                } else {
                    double x2 = c1.x - (c1.y - ymin) / m;
                    return Crossing{Side::BOTTOM, clamp(x2, xmin, xmax), ymin};
                }
            }
        }

    }

    bool Box::contains(const Box& b) const {
        return b.xmin >= xmin && b.xmax <= xmax && b.ymin >= ymin && b.ymax <= ymax;
    }

    bool Box::contains(const Coordinate &c) const {
        return c.x >= xmin && c.x <= xmax && c.y >= ymin && c.y <= ymax;
    }

    bool Box::strictly_contains(const Coordinate &c) const {
        return c.x > xmin && c.x < xmax && c.y > ymin && c.y < ymax;
    }

    std::ostream &operator<<(std::ostream &os, const Box &b) {
        os << "POLYGON ((";
        os << b.xmin << " " << b.ymin << ", ";
        os << b.xmax << " " << b.ymin << ", ";
        os << b.xmax << " " << b.ymax << ", ";
        os << b.xmin << " " << b.ymax << ", ";
        os << b.xmin << " " << b.ymin << "))";
        return os;
    }

}

// Copyright (c) 2018 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_COORDINATE_H
#define EXACTEXTRACT_COORDINATE_H

#include <cmath>
#include <iostream>

namespace exactextract {

    struct Coordinate {
        double x;
        double y;

        Coordinate() = default;

        Coordinate(double p_x, double p_y) : x{p_x}, y{p_y} {}

        bool equals(const Coordinate &other, double tol) const {
            return std::abs(other.x - x) < tol && std::abs(other.y - y) < tol;
        }

        bool operator==(const Coordinate &other) const {
            return x == other.x && y == other.y;
        }

        bool operator!=(const Coordinate &other) const {
            return !(*this == other);
        }
    };

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

}

#endif

// Copyright (c) 2018 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_CELL_CROSSING_H
#define EXACTEXTRACT_CELL_CROSSING_H

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

namespace exactextract {

    class Crossing {
    public:
        Crossing(Side s, double x, double y) : m_side{s}, m_coord{x, y} {}

        Crossing(Side s, const Coordinate &c) : m_side{s}, m_coord{c} {}

        const Side &side() const {
            return m_side;
        }

        const Coordinate &coord() const {
            return m_coord;
        }

    private:
        Side m_side;
        Coordinate m_coord;

    };

}

#endif

// Copyright (c) 2018 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 <stdexcept>

#include "area.h"
#include "cell.h"
#include "crossing.h"
#include "traversal_areas.h"
#include "geos_utils.h"

namespace exactextract {

    double Cell::height() const {
        return m_box.height();
    }

    double Cell::width() const {
        return m_box.width();
    }

    double Cell::area() const {
        return m_box.area();
    }

    Side Cell::side(const Coordinate &c) const {
        if (c.x == m_box.xmin) {
            return Side::LEFT;
        } else if (c.x == m_box.xmax) {
            return Side::RIGHT;
        } else if (c.y == m_box.ymin) {
            return Side::BOTTOM;
        } else if (c.y == m_box.ymax) {
            return Side::TOP;
        }

        return Side::NONE;
    }

    void Cell::force_exit() {
        if (last_traversal().exited()) {
            return;
        }

        const Coordinate &last = last_traversal().last_coordinate();

        if (location(last) == Location::BOUNDARY) {
            last_traversal().force_exit(side(last));
        }
    }

    Cell::Location Cell::location(const Coordinate &c) const {
        if (m_box.strictly_contains(c)) {
            return Cell::Location::INSIDE;
        }

        if (m_box.contains(c)) {
            return Cell::Location::BOUNDARY;
        }

        return Cell::Location::OUTSIDE;
    }

    Traversal &Cell::traversal_in_progress() {
        if (m_traversals.empty() || m_traversals[m_traversals.size() - 1].exited()) {
            m_traversals.emplace_back();
        }

        return m_traversals[m_traversals.size() - 1];
    }

    Traversal &Cell::last_traversal() {
        return m_traversals.at(m_traversals.size() - 1);
    }

    bool Cell::take(const Coordinate &c) {
        Traversal &t = traversal_in_progress();

        if (t.empty()) {
            //std::cout << "Entering " << m_box << " from " << side(c) << " at " << c << std::endl;

            t.enter(c, side(c));
            return true;
        }

        if (location(c) != Cell::Location::OUTSIDE) {
            //std::cout << "Still in " << m_box << " with " << c << std::endl;

            t.add(c);
            return true;
        }

        Crossing x = m_box.crossing(t.last_coordinate(), c);
        t.exit(x.coord(), x.side());

        //std::cout << "Leaving " << m_box << " from " << x.side() << " at " << x.coord();
        //std::cout << " on the way to " << c << std::endl;

        return false;
    }

    double Cell::covered_fraction() const {
        // Handle the special case of a ring that is enclosed within a
        // single pixel of our raster
        if (m_traversals.size() == 1 && m_traversals[0].is_closed_ring()) {
            return exactextract::area(m_traversals[0].coords()) / area();
        }

        // TODO consider porting in simplified single-traversal area calculations
        // from Java code. Do they really make a performance difference?
        //if (m_traversals.size() == 1) {
        //    double a = area();

        //    return (a - area_right_of(m_traversals.at(m_traversals.size() - 1))) / a;
        //}

        std::vector<const std::vector<Coordinate> *> coord_lists;

        for (const auto &t : m_traversals) {
            if (!t.traversed() || !t.multiple_unique_coordinates()) {
                continue;
            }

            coord_lists.push_back(&t.coords());
        }

        return left_hand_area(m_box, coord_lists) / area();
    }

#if 0
    Crossing Cell::crossing(const Coordinate & c1, const Coordinate & c2) const {
        Coordinate result(0, 0);

        if (c2.y > c1.y && segment_intersection(c1, c2, m_box.upper_left(), m_box.upper_right(), result)) {
            return Crossing(Side::TOP, result);
        }

        if (c2.y < c1.y && segment_intersection(c1, c2, m_box.lower_right(), m_box.lower_left(), result)) {
            return Crossing(Side::BOTTOM, result);
        }

        if (c2.x < c1.x && segment_intersection(c1, c2, m_box.lower_left(), m_box.upper_left(), result)) {
            return Crossing(Side::LEFT, result);
        }

        if (c2.x > c1.x && segment_intersection(c1, c2, m_box.lower_right(), m_box.upper_right(), result)) {
            return Crossing(Side::RIGHT, result);
        }

        throw std::runtime_error("Never get here!");
    }
#endif


}

// Copyright (c) 2018 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_TRAVERSAL_H
#define EXACTEXTRACT_TRAVERSAL_H

#include <vector>

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

namespace exactextract {

    class Traversal {
    public:
        Traversal() : m_entry{Side::NONE}, m_exit{Side::NONE} {}

        bool is_closed_ring() const;

        bool empty() const;

        bool entered() const;

        bool exited() const;

        bool traversed() const;

        bool multiple_unique_coordinates() const;

        void enter(const Coordinate &c, Side s);

        void exit(const Coordinate &c, Side s);

        Side entry_side() const { return m_entry; }

        Side exit_side() const { return m_exit; }

        const Coordinate &last_coordinate() const;

        const Coordinate &exit_coordinate() const;

        void add(const Coordinate &c);

        void force_exit(Side s) { m_exit = s; }

        const std::vector<Coordinate> &coords() const { return m_coords; }

    private:
        std::vector<Coordinate> m_coords;
        Side m_entry;
        Side m_exit;
    };

}

#endif

// Copyright (c) 2018 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 "coordinate.h"

namespace exactextract {

    std::ostream &operator<<(std::ostream &os, const Coordinate &c) {
        os << "POINT (" << c.x << " " << c.y << ")";
        return os;
    }

}

// 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 <geos_c.h>

#include "grid.h"
#include "floodfill.h"
#include "geos_utils.h"

namespace exactextract {

    FloodFill::FloodFill(GEOSContextHandle_t context, const GEOSGeometry *g, const Grid<bounded_extent> &extent) :
            m_extent{extent},
            m_geos_context{context},
            m_g{nullptr},
            m_pg{nullptr} {
        geom_ptr_r ring_copy = geos_ptr(context, GEOSGeom_clone_r(context, g));
        m_g = geos_ptr(context, GEOSGeom_createPolygon_r(context, ring_copy.release(), nullptr, 0));
        m_pg = GEOSPrepare_ptr(context, m_g.get());
    }

    bool FloodFill::cell_is_inside(size_t i, size_t j) const {
        double x = m_extent.x_for_col(j);
        double y = m_extent.y_for_row(i);

        auto point = GEOSGeom_createPoint_ptr(m_geos_context, x, y);

        return GEOSPreparedContains_r(m_geos_context, m_pg.get(), point.get());
    }

}

// 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_GEOS_UTILS_H
#define EXACTEXTRACT_GEOS_UTILS_H

#include <algorithm>
#include <functional>
#include <limits>
#include <memory>
#include <stdexcept>
#include <vector>

#include <geos_c.h>

#define HAVE_370 (GEOS_VERSION_MAJOR >= 3 && GEOS_VERSION_MINOR >= 7)
#define HAVE_380 (GEOS_VERSION_MAJOR >= 3 && GEOS_VERSION_MINOR >= 8)

#include "box.h"
#include "coordinate.h"

namespace exactextract {

    using geom_ptr_r = std::unique_ptr<GEOSGeometry, std::function<void(GEOSGeometry*)>>;
    using tree_ptr_r = std::unique_ptr<GEOSSTRtree, std::function<void(GEOSSTRtree*)>>;
    using seq_ptr_r = std::unique_ptr<GEOSCoordSequence, std::function<void(GEOSCoordSequence*)>>;
    using prep_geom_ptr_r = std::unique_ptr<const GEOSPreparedGeometry, std::function<void(const GEOSPreparedGeometry*)>>;

    inline geom_ptr_r geos_ptr(GEOSContextHandle_t context, GEOSGeometry* geom) {
        auto deleter = [context](GEOSGeometry* g){ GEOSGeom_destroy_r(context, g); };
        return geom_ptr_r{geom, deleter};
    }

    inline tree_ptr_r geos_ptr(GEOSContextHandle_t context, GEOSSTRtree* tree) {
        auto deleter = [context](GEOSSTRtree_t* t){ GEOSSTRtree_destroy_r(context, t); };
        return tree_ptr_r{tree, deleter};
    }

    inline seq_ptr_r geos_ptr(GEOSContextHandle_t context, GEOSCoordSequence* seq) {
        auto deleter = [context](GEOSCoordSequence* s){ GEOSCoordSeq_destroy_r(context, s); };
        return seq_ptr_r{seq, deleter};
    }

    inline prep_geom_ptr_r
    GEOSPrepare_ptr(GEOSContextHandle_t context, const GEOSGeometry *g) {
        auto deleter = [context](const GEOSPreparedGeometry* pg) { GEOSPreparedGeom_destroy_r(context, pg); };
        return prep_geom_ptr_r{GEOSPrepare_r(context, g), deleter};
    }

    inline seq_ptr_r
    GEOSCoordSeq_create_ptr(GEOSContextHandle_t context, unsigned int size, unsigned int dims) {
        return geos_ptr(context, GEOSCoordSeq_create_r(context, size, dims));
    }

    inline geom_ptr_r
    GEOSGeom_createPoint_ptr(GEOSContextHandle_t context, double x, double y) {
#if HAVE_380
        return geos_ptr(context, GEOSGeom_createPointFromXY_r(context, x, y));
#else
        auto seq = GEOSCoordSeq_create_ptr(context, 1, 2);
        GEOSCoordSeq_setX_r(context, seq.get(), 0, x);
        GEOSCoordSeq_setY_r(context, seq.get(), 0, y);
        return geos_ptr(context, GEOSGeom_createPoint_r(context, seq.release()));
#endif
    }

    inline geom_ptr_r
    GEOSGeom_createLineString_ptr(GEOSContextHandle_t context, GEOSCoordSequence *seq) {
        return geos_ptr(context, GEOSGeom_createLineString_r(context, seq));
    }

    inline unsigned int geos_get_num_points(GEOSContextHandle_t context, const GEOSCoordSequence *s) {
        unsigned int result;
        if (!GEOSCoordSeq_getSize_r(context, s, &result)) {
            throw std::runtime_error("Error calling GEOSCoordSeq_getSize_r.");
        }
        return result;
    }

    inline geom_ptr_r
    GEOSGeom_read_r(GEOSContextHandle_t context, const std::string &s) {
        return geos_ptr(context, GEOSGeomFromWKT_r(context, s.c_str()));
    }
    geom_ptr_r geos_make_box_polygon(GEOSContextHandle_t context, const Box & b);

    Box geos_get_box(GEOSContextHandle_t context, const GEOSGeometry* g);

    std::vector<Box> geos_get_component_boxes(GEOSContextHandle_t context, const GEOSGeometry* g);

    bool segment_intersection(GEOSContextHandle_t context, const Coordinate &a0, const Coordinate &a1, const Coordinate &b0, const Coordinate &b1,
                              Coordinate &result);

    bool geos_is_ccw(GEOSContextHandle_t context, const GEOSCoordSequence *s);

    std::vector<Coordinate> read(GEOSContextHandle_t context, const GEOSCoordSequence *s);

}

#endif //RASTER_OVERLAY_CPP_GEOS_UTILS_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.

#include "geos_utils.h"

#if !HAVE_380
static inline int GEOSCoordSeq_setXY_r(GEOSContextHandle_t context,
        GEOSCoordSequence* seq,
        unsigned int idx,
        double x,
        double y) {
    return GEOSCoordSeq_setX_r(context, seq, idx, x) && GEOSCoordSeq_setY_r(context, seq, idx, y);
}
#endif

namespace exactextract {

    geom_ptr_r geos_make_box_polygon(GEOSContextHandle_t context, const Box & b) {
        auto seq = geos_ptr(context, GEOSCoordSeq_create_r(context, 5, 2));

        GEOSCoordSeq_setXY_r(context, seq.get(), 0, b.xmin, b.ymin);
        GEOSCoordSeq_setXY_r(context, seq.get(), 1, b.xmax, b.ymin);
        GEOSCoordSeq_setXY_r(context, seq.get(), 2, b.xmax, b.ymax);
        GEOSCoordSeq_setXY_r(context, seq.get(), 3, b.xmin, b.ymax);
        GEOSCoordSeq_setXY_r(context, seq.get(), 4, b.xmin, b.ymin);

        auto shell = geos_ptr(context, GEOSGeom_createLinearRing_r(context, seq.release()));

        return geos_ptr(context, GEOSGeom_createPolygon_r(context, shell.release(), nullptr, 0));
    }

    bool segment_intersection(
            GEOSContextHandle_t context,
            const Coordinate &a0,
            const Coordinate &a1,
            const Coordinate &b0,
            const Coordinate &b1,
            Coordinate &result) {
#if HAVE_370
        int code = GEOSSegmentIntersection_r(context,
                                             a0.x, a0.y,
                                             a1.x, a1.y,
                                             b0.x, b0.y,
                                             b1.x, b1.y,
                                             &result.x, &result.y);
        if (!code) {
            throw std::runtime_error("Error in GEOSSegmentIntersection_r");
        }

        return code == 1;
#else
        auto seqa = GEOSCoordSeq_create_ptr(context, 2, 2);
        auto seqb = GEOSCoordSeq_create_ptr(context, 2, 2);

        GEOSCoordSeq_setX_r(context, seqa.get(), 0, a0.x);
        GEOSCoordSeq_setY_r(context, seqa.get(), 0, a0.y);
        GEOSCoordSeq_setX_r(context, seqa.get(), 1, a1.x);
        GEOSCoordSeq_setY_r(context, seqa.get(), 1, a1.y);

        GEOSCoordSeq_setX_r(context, seqb.get(), 0, b0.x);
        GEOSCoordSeq_setY_r(context, seqb.get(), 0, b0.y);
        GEOSCoordSeq_setX_r(context, seqb.get(), 1, b1.x);
        GEOSCoordSeq_setY_r(context, seqb.get(), 1, b1.y);

        auto geom_a = GEOSGeom_createLineString_ptr(context, seqa.release());
        auto geom_b = GEOSGeom_createLineString_ptr(context, seqb.release());

        geom_ptr_r intersection = geos_ptr(context, GEOSIntersection_r(context, geom_a.get(), geom_b.get()));

        if (GEOSisEmpty_r(context, intersection.get())) {
            return false;
        }

        if (GEOSGeomTypeId_r(context, intersection.get()) != GEOS_POINT) {
            return false;
        }

        GEOSGeomGetX_r(context, intersection.get(), &result.x);
        GEOSGeomGetY_r(context, intersection.get(), &result.y);

        return true;
#endif
    }

    Box geos_get_box(GEOSContextHandle_t context, const GEOSGeometry* g) {
        double xmin, ymin, xmax, ymax;

#if HAVE_370
        if (!(GEOSGeom_getXMin_r(context, g, &xmin) &&
              GEOSGeom_getYMin_r(context, g, &ymin) &&
              GEOSGeom_getXMax_r(context, g, &xmax) &&
              GEOSGeom_getYMax_r(context, g, &ymax))) {
            throw std::runtime_error("Error getting geometry extent.");
        }
#else
        geom_ptr_r env = geos_ptr(context, GEOSEnvelope_r(context, g));

        const GEOSGeometry* ring = GEOSGetExteriorRing_r(context, env.get());
        const GEOSCoordSequence* seq = GEOSGeom_getCoordSeq_r(context, ring);

        xmin = std::numeric_limits<double>::max();
        ymin = std::numeric_limits<double>::max();
        xmax = std::numeric_limits<double>::lowest();
        ymax = std::numeric_limits<double>::lowest();

        for (unsigned int i = 0; i < 4; i++) {
            double x, y;

            if (!GEOSCoordSeq_getX_r(context, seq, i, &x) || !GEOSCoordSeq_getY_r(context, seq, i, &y)) {
                throw std::runtime_error("Error reading coordinates.");
            }

            xmin = std::min(xmin, x);
            ymin = std::min(ymin, y);
            xmax = std::max(xmax, x);
            ymax = std::max(ymax, y);
        }
#endif
        return {xmin, ymin, xmax, ymax};
    }

    std::vector<Box> geos_get_component_boxes(GEOSContextHandle_t context, const GEOSGeometry* g) {
        size_t n = static_cast<size_t>(GEOSGetNumGeometries_r(context, g));
        std::vector<Box> boxes;
        boxes.reserve(n);

        for (size_t i = 0; i < n; i++) {
            boxes.push_back(geos_get_box(context, GEOSGetGeometryN_r(context, g, i)));
        }

        return boxes;
    }

    bool geos_is_ccw(GEOSContextHandle_t context, const GEOSCoordSequence *s) {
#if HAVE_370
        char result;
        if (!GEOSCoordSeq_isCCW_r(context, s, &result)) {
            throw std::runtime_error("Error calling GEOSCoordSeq_isCCW_r.");
        }
        return result;
#else
        std::vector<Coordinate> coords = read(context, s);

        if (coords.size() < 4) {
            throw std::runtime_error("Ring has fewer than 4 points, so orientation cannot be determined.");
        }

        // find highest point
        size_t hi_index = (size_t) std::distance(
                coords.begin(),
                std::max_element(coords.begin(), coords.end(), [](const auto& a, const auto&b) {
                    return a.y < b.y;
                })
        );

        // find distinct point before highest point
        size_t i_prev = hi_index;
        do {
            if (i_prev == 0) {
                i_prev = coords.size() - 1;
            } else {
                i_prev--;
            }
        } while (i_prev != hi_index && coords[i_prev] == coords[hi_index]);

        // find distinct point after highest point
        size_t i_next = hi_index;
        do {
            i_next = (i_next + 1) % coords.size();
        } while (i_next != hi_index && coords[i_next] == coords[hi_index]);

        Coordinate& a = coords[i_prev];
        Coordinate& b = coords[hi_index];
        Coordinate& c = coords[i_next];

        if (a == b || b == c || a == c) {
            return false;
        }

        int disc = GEOSOrientationIndex_r(context, a.x, a.y, b.x, b.y, c.x, c.y);

        if (disc == 0) {
            // poly is CCW if prev x is right of next x
            return (a.x > b.x);
        } else {
            // if area is positive, points are ordered CCW
            return disc > 0;
        }
#endif
    }

    std::vector<Coordinate> read(GEOSContextHandle_t context, const GEOSCoordSequence *s) {
        unsigned int size;

        if (!GEOSCoordSeq_getSize_r(context, s, &size)) {
            throw std::runtime_error("Error calling GEOSCoordSeq_getSize.");
        }

        std::vector<Coordinate> coords{size};

        for (unsigned int i = 0; i < size; i++) {
            if (!GEOSCoordSeq_getX_r(context, s, i, &(coords[i].x)) || !GEOSCoordSeq_getY_r(context, s, i, &(coords[i].y))) {
                throw std::runtime_error("Error reading coordinates.");
            }
        }

        return coords;
    }

}

// Copyright (c) 2018 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 "grid.h"

namespace exactextract {
    Box grid_cell(const Grid<bounded_extent> & grid, size_t row, size_t col) {
        // The ternary clauses below are used to make sure that the cells along
        // the right and bottom edges of our grid are slightly larger than m_dx,dy
        // if needed to make sure that we capture our whole extent. This is necessary
        // because xmin + nx*m_dx may be less than xmax because of floating point
        // errors.
        return {
                grid.xmin() + col * grid.dx(),
                row == (grid.rows() - 1) ? grid.ymin() : (grid.ymax() - (row + 1) * grid.dy()),
                col == (grid.cols() - 1) ? grid.xmax() : (grid.xmin() + (col + 1) * grid.dx()),
                grid.ymax() - row * grid.dy()

        };
    }

    Box grid_cell(const Grid<infinite_extent> & grid, size_t row, size_t col) {
        double xmin, xmax, ymin, ymax;

        if (col == 0) {
            xmin = std::numeric_limits<double>::lowest();
        } else if (col == grid.cols() - 1) {
            xmin = grid.xmax(); // because rightmost col of regular grid may have different width from others
        } else  {
            xmin = grid.xmin() + (col - 1) * grid.dx();
        }

        switch(grid.cols() - col) {
            case 1: xmax = std::numeric_limits<double>::max(); break;
            case 2: xmax = grid.xmax(); break;
            default: xmax = grid.xmin() + col*grid.dx();
        }

        if (row == 0) {
            ymax = std::numeric_limits<double>::max();
        } else if (row == grid.rows() - 1) {
            ymax = grid.ymin(); // because bottom row of regular grid may have different height from others
        } else {
            ymax = grid.ymax() - (row - 1) * grid.dy();
        }

        switch(grid.rows() - row) {
            case 1: ymin = std::numeric_limits<double>::lowest(); break;
            case 2: ymin = grid.ymin(); break;
            default: ymin = grid.ymax() - row*grid.dy();
        }

        return { xmin, ymin, xmax, ymax };
    }

    Grid<infinite_extent> make_infinite(const Grid<bounded_extent> & grid) {
        return { grid.extent(), grid.dx(), grid.dy() };
    }

    Grid<bounded_extent> make_finite(const Grid<infinite_extent> & grid) {
        return { grid.extent(), grid.dx(), grid.dy() };
    }

    std::vector<Grid<bounded_extent>> subdivide(const Grid<bounded_extent> & grid, size_t max_size) {
        if (grid.size() < max_size) {
            return { grid };
        }

        size_t cols_per_block = std::min(max_size, grid.cols());
        size_t rows_per_block = max_size / cols_per_block;

        size_t col_blocks = (grid.cols() - 1) / cols_per_block + 1;
        size_t row_blocks = (grid.rows() - 1) / rows_per_block + 1;

        std::vector<Grid<bounded_extent>> subgrids;
        for (size_t i = 0; i < row_blocks; i++) {
            for (size_t j = 0; j < col_blocks; j++) {
                double xmin = grid.xmin() + grid.dx()*cols_per_block*j;
                double xmax = j == (col_blocks - 1) ? grid.xmax() : (grid.xmin() + grid.dx()*cols_per_block*(j+1));
                double ymax = grid.ymax() - grid.dy()*rows_per_block*i;
                double ymin = i == (row_blocks - 1) ? grid.ymin() : (grid.ymax() - grid.dy()*rows_per_block*(i+1));

                Box reduced = {xmin, ymin, xmax, ymax};
                subgrids.emplace_back(reduced, grid.dx(), grid.dy());
            }
        }

        return subgrids;
    }
}

// Copyright (c) 2018 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 <stdexcept>

#include "perimeter_distance.h"

namespace exactextract {

    double perimeter_distance(double xmin, double ymin, double xmax, double ymax, double x, double y) {
        if (x == xmin) {
            // Left
            return y - ymin;
        }

        if (y == ymax) {
            // Top
            return (ymax - ymin) + x - xmin;
        }

        if (x == xmax) {
            // Right
            return (xmax - xmin) + (ymax - ymin) + ymax - y;
        }

        if (y == ymin) {
            // Bottom
            return (xmax - xmin) + 2 * (ymax - ymin) + (xmax - x);
        }

        throw std::runtime_error("Never get here");
    }

    double perimeter_distance(double xmin, double ymin, double xmax, double ymax, const Coordinate &c) {
        return perimeter_distance(xmin, ymin, xmax, ymax, c.x, c.y);
    }

    double perimeter_distance(const Box &b, const Coordinate &c) {
        return perimeter_distance(b.xmin, b.ymin, b.xmax, b.ymax, c.x, c.y);
    }

    double perimeter_distance_ccw(double measure1, double measure2, double perimeter) {
        if (measure2 <= measure1) {
            return measure1 - measure2;
        }
        return perimeter + measure1 - measure2;
    }

}

// Copyright (c) 2018 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_CELL_H
#define EXACTEXTRACT_CELL_H

#include <memory>

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

namespace exactextract {

    class Cell {

    public:

        Cell(double xmin, double ymin, double xmax, double ymax) :
                m_box{xmin, ymin, xmax, ymax} {}

        explicit Cell(const Box & b) : m_box{b} {}

        void force_exit();

        double width() const;

        double height() const;

        double area() const;

        double covered_fraction() const;

        Traversal &last_traversal();

        const Box &box() const { return m_box; }

        bool take(const Coordinate &c);

    private:
        enum class Location {
            INSIDE, OUTSIDE, BOUNDARY
        };

        Box m_box;

        std::vector<Traversal> m_traversals;

        Side side(const Coordinate &c) const;

        Location location(const Coordinate &c) const;

        Traversal &traversal_in_progress();
    };

}

#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_FLOODFILL_H
#define EXACTEXTRACT_FLOODFILL_H

#include <queue>

#include <geos_c.h>

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

namespace exactextract {

    template<typename T>
    struct fill_values {
        static T EXTERIOR;
    };

    template<>
    struct fill_values<float> {
        static constexpr float EXTERIOR{0.0f};  // Cell is known to be entirely outside the polygon
        static constexpr float INTERIOR{1.0f};  // Cell is known to be entirely within the polygon
        static constexpr float FILLABLE{-1.0f}; // Cell location relative to polygon unknown, but
                                                // can be determined by fill.
        static constexpr float UNKNOWN{-2.0f};  // Cell location relative to polygon unknown
                                                // and cannot be determined from a flood fill
                                                // (must be explicitly tested)

    };

    class FloodFill {

    public:
        FloodFill(GEOSContextHandle_t context, const GEOSGeometry *g, const Grid<bounded_extent> &extent);

        template<typename T>
        void flood(Matrix<T> &arr) const;

        bool cell_is_inside(size_t i, size_t j) const;

    private:
        Grid<bounded_extent> m_extent;
        GEOSContextHandle_t m_geos_context;
        geom_ptr_r m_g;
        prep_geom_ptr_r m_pg;
    };

    template<typename T>
    void flood_from_pixel(Matrix<T> &arr, size_t i, size_t j, T fill_value) {
        std::queue<std::pair<size_t, size_t> > locations;

        locations.emplace(i, j);

        while (!locations.empty()) {
            i = locations.front().first;
            j = locations.front().second;
            locations.pop();

            if (arr(i, j) == fill_value) {
                continue;
            }

            // Left
            if (j > 0 && arr(i, j - 1) == fill_values<T>::FILLABLE) {
                locations.emplace(i, j - 1);
            }

            auto j0 = j;

            // Fill along this row until we hit something
            for (; j < arr.cols() && arr(i, j) == fill_values<T>::FILLABLE; j++) {
                arr(i, j) = fill_value;
            }

            auto j1 = j;

            // Initiate scanlines above our current row
            if (i > 0) {
                for (j = j0; j < j1; j++) {
                    // Up
                    if (arr(i - 1, j) == fill_values<T>::FILLABLE) {
                        locations.emplace(i - 1, j);
                    }
                }
            }

            // Initiate scanlines below our current row
            if (i < arr.rows() - 1) {
                for (j = j0; j < j1; j++) {
                    // Down
                    if (arr(i + 1, j) == fill_values<T>::FILLABLE) {
                        locations.emplace(i + 1, j);
                    }
                }
            }

        }
    }

    template<typename T>
    void FloodFill::flood(Matrix<T> &arr) const {

        for (size_t i = 0; i < arr.rows(); i++) {
            for (size_t j = 0; j < arr.cols(); j++) {
                if (arr(i, j) == fill_values<T>::UNKNOWN) {
                    throw std::runtime_error("Cell with unknown position encountered.");
                } else if (arr(i, j) == fill_values<T>::FILLABLE) {
                    // Cell position relative to polygon is unknown but can
                    // be determined from adjacent cells.
                    if (cell_is_inside(i, j)) {
                        flood_from_pixel(arr, i, j, fill_values<T>::INTERIOR);
                    } else {
                        flood_from_pixel(arr, i, j, fill_values<T>::EXTERIOR);
                    }
                }
            }
        }
    }

}

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

#include <stdexcept>

#include <geos_c.h>

#include "area.h"
#include "cell.h"
#include "floodfill.h"
#include "geos_utils.h"
#include "raster_cell_intersection.h"

namespace exactextract {

    Raster<float> raster_cell_intersection(const Grid<bounded_extent> & raster_grid, GEOSContextHandle_t context, const GEOSGeometry* g) {
        RasterCellIntersection rci(raster_grid, context, g);

        return { std::move(const_cast<Matrix<float>&>(rci.overlap_areas())),
                 make_finite(rci.m_geometry_grid) };
    }

    Raster<float> raster_cell_intersection(const Grid<bounded_extent> & raster_grid, const Box & box) {
        RasterCellIntersection rci(raster_grid, box);

        return { std::move(const_cast<Matrix<float>&>(rci.overlap_areas())),
                 make_finite(rci.m_geometry_grid) };
    }

    static Cell *get_cell(Matrix<std::unique_ptr<Cell>> &cells, const Grid<infinite_extent> &ex, size_t row, size_t col) {
        //std::cout << " getting cell " << row << ", " << col << std::endl;

        if (cells(row, col) == nullptr) {
            cells(row, col) = std::make_unique<Cell>(grid_cell(ex, row, col));
        }

        return cells(row, col).get();
    }

    Box processing_region(const Box & raster_extent, const std::vector<Box> & component_boxes) {
        Box ret = Box::make_empty();
        for (const auto& box : component_boxes) {
            if (ret == raster_extent) {
                // No more expansion is possible
                return ret;
            }

            if (!box.intersects(raster_extent)) {
                continue;
            }

            Box isect = raster_extent.intersection(box);
            if (ret.empty()) {
                ret = isect;
            } else if (!ret.contains(isect)) {
                ret = ret.expand_to_include(isect);
            }
        }

        return ret;
    }

    static Grid<infinite_extent> get_geometry_grid(const Grid<bounded_extent> &raster_grid, GEOSContextHandle_t context, const GEOSGeometry* g) {
        if (GEOSisEmpty_r(context, g)) {
            throw std::invalid_argument("Can't get statistics for empty geometry");
        }

        Box region = processing_region(raster_grid.extent(), geos_get_component_boxes(context, g));

        if (!region.empty()) {
            return make_infinite(raster_grid.shrink_to_fit(region));
        } else {
            return Grid<infinite_extent>::make_empty();
        }
    }

    static Grid<infinite_extent> get_geometry_grid(const Grid<bounded_extent> & raster_grid, const Box & box) {
        auto region = box.intersection(raster_grid.extent());

        if (!region.empty()) {
            return make_infinite(raster_grid.shrink_to_fit(region));
        } else {
            return Grid<infinite_extent>::make_empty();
        }
    }


    RasterCellIntersection::RasterCellIntersection(const Grid<bounded_extent> &raster_grid, GEOSContextHandle_t context, const GEOSGeometry *g)
        : m_geometry_grid{get_geometry_grid(raster_grid, context, g)},
          m_overlap_areas{std::make_unique<Matrix<float>>(m_geometry_grid.rows() - 2, m_geometry_grid.cols() - 2)}
    {
        if (!m_geometry_grid.empty())
            process(context, g);
    }

    RasterCellIntersection::RasterCellIntersection(const Grid<bounded_extent> & raster_grid, const Box & box)
        : m_geometry_grid{get_geometry_grid(raster_grid, box)},
          m_overlap_areas{std::make_unique<Matrix<float>>(m_geometry_grid.rows() - 2, m_geometry_grid.cols() - 2)} {
        if (!m_geometry_grid.empty()) {
            process_rectangular_ring(box, true);
        }
    }

    void RasterCellIntersection::process(GEOSContextHandle_t context, const GEOSGeometry *g) {
        if (GEOSGeomTypeId_r(context, g) == GEOS_POLYGON) {
            process_ring(context, GEOSGetExteriorRing_r(context, g), true);

            for (int i = 0; i < GEOSGetNumInteriorRings_r(context, g); i++) {
                process_ring(context, GEOSGetInteriorRingN_r(context, g, i), false);
            }
        } else if (GEOSGeomTypeId_r(context, g) == GEOS_MULTIPOLYGON) {
            for (int i = 0; i < GEOSGetNumGeometries_r(context, g); i++) {
                process(context, GEOSGetGeometryN_r(context, g, i));
            }
        } else {
            throw std::invalid_argument("Unsupported geometry type.");
        }
    }

    static Grid<infinite_extent> get_box_grid(const Box & box, const Grid<infinite_extent> & geometry_grid) {
        Box cropped_ring_extent = geometry_grid.extent().intersection(box);
        return geometry_grid.shrink_to_fit(cropped_ring_extent);
    }

    static Grid<infinite_extent> get_ring_grid(GEOSContextHandle_t context, const GEOSGeometry* ls, const Grid<infinite_extent> & geometry_grid) {
        return get_box_grid(geos_get_box(context, ls), geometry_grid);
    }

    void RasterCellIntersection::process_rectangular_ring(const Box& box, bool exterior_ring) {
        if (!box.intersects(m_geometry_grid.extent())) {
            return;
        }

        auto ring_grid = get_box_grid(box, m_geometry_grid);

        auto row_min = ring_grid.get_row(box.ymax);
        auto row_max = ring_grid.get_row(box.ymin);
        auto col_min = ring_grid.get_column(box.xmin);
        auto col_max = ring_grid.get_column(box.xmax);

        Matrix<float> areas(ring_grid.rows() - 2, ring_grid.cols() - 2);

        // upper-left
        if (row_min > 0 && col_min > 0) {
            auto ul = grid_cell(ring_grid, row_min, col_min);
            areas(row_min - 1, col_min - 1) = ul.intersection(box).area() / ul.area();
        }

        // upper-right
        if (row_min > 0 && col_max < ring_grid.cols() - 1) {
            auto ur = grid_cell(ring_grid, row_min, col_max);
            auto frac = ur.intersection(box).area() / ur.area();
            areas(row_min - 1, col_max - 1) = frac;
        }

        // lower-left
        if (row_max < ring_grid.rows() - 1 && col_min > 0) {
            auto ll = grid_cell(ring_grid, row_max, col_min);
            areas(row_max - 1, col_min - 1) = ll.intersection(box).area() / ll.area();
        }

        // lower-right
        if (row_max < ring_grid.rows() - 1 && col_max < ring_grid.cols() - 1) {
            auto lr = grid_cell(ring_grid, row_max, col_max);
            areas(row_max - 1, col_max - 1) = lr.intersection(box).area() / lr.area();
        }

        // left
        if (col_min > 0) {
            auto left = grid_cell(ring_grid, row_min + 1, col_min);
            auto frac = left.intersection(box).area() / left.area();
            for (size_t row = row_min + 1; row < row_max; row++) {
                areas(row - 1, col_min - 1) = frac;
            }
        }

        // right
        if (col_max < ring_grid.cols() - 1) {
            auto right = grid_cell(ring_grid, row_min + 1, col_max);
            auto frac = right.intersection(box).area() / right.area();
            for (size_t row = row_min + 1; row < row_max; row++) {
                areas(row - 1, col_max - 1) = frac;
            }
        }

        // top
        if (row_min > 0) {
            auto top = grid_cell(ring_grid, row_min, col_min + 1);
            auto frac = top.intersection(box).area() / top.area();
            for (size_t col = col_min + 1; col < col_max; col++) {
                areas(row_min - 1, col - 1) = frac;
            }
        }

        // bottom
        if (row_max < ring_grid.rows() - 1) {
            auto bottom = grid_cell(ring_grid, row_max, col_min + 1);
            auto frac = bottom.intersection(box).area() / bottom.area();
            for (size_t col = col_min + 1; col < col_max; col++) {
                areas(row_max - 1, col - 1) = frac;
            }
        }

        // interior
        for (size_t row = row_min + 1; row < row_max; row++) {
            for(size_t col = col_min + 1; col < col_max; col++) {
                areas(row - 1, col - 1) = 1.0f;
            }
        }

        // Transfer these areas to our global set
        size_t i0 = ring_grid.row_offset(m_geometry_grid);
        size_t j0 = ring_grid.col_offset(m_geometry_grid);

        add_ring_areas(i0, j0, areas, exterior_ring);
    }

    void RasterCellIntersection::process_ring(GEOSContextHandle_t context, const GEOSGeometry *ls, bool exterior_ring) {
        auto geom_box = geos_get_box(context, ls);

        if (!geom_box.intersects(m_geometry_grid.extent())) {
            return;
        }

        const GEOSCoordSequence *seq = GEOSGeom_getCoordSeq_r(context, ls);
        unsigned int npoints = geos_get_num_points(context, seq);

        if (npoints == 5) {
            auto coords = read(context, seq);
            if (area(coords) == geom_box.area()) {
                process_rectangular_ring(geom_box, exterior_ring);
                return;
            }
        }

        Grid<infinite_extent> ring_grid = get_ring_grid(context, ls, m_geometry_grid);

        size_t rows = ring_grid.rows();
        size_t cols = ring_grid.cols();

        // Short circuit for small rings that are entirely contained
        // within a single grid cell.
        if (rows == (1 + 2*infinite_extent::padding) &&
            cols == (1 + 2*infinite_extent::padding) &&
            grid_cell(ring_grid, 1, 1).contains(geom_box)) {

            auto coords = read(context, seq);
            auto ring_area = area(coords) / grid_cell(ring_grid, 1, 1).area();

            size_t i0 = ring_grid.row_offset(m_geometry_grid);
            size_t j0 = ring_grid.col_offset(m_geometry_grid);

            if (exterior_ring) {
                m_overlap_areas->increment(i0, j0, ring_area);
            } else {
                m_overlap_areas->increment(i0, j0, -1 * ring_area);
            }

            return;
        }

        bool is_ccw = geos_is_ccw(context, seq);
        Matrix<std::unique_ptr<Cell>> cells(rows, cols);

        std::deque<Coordinate> stk;
        {
            for (unsigned int i = 0; i < npoints; i++) {
                double x, y;
#if HAVE_380
                if (!GEOSCoordSeq_getXY_r(context, seq, i, &x, &y)) {
                    throw std::runtime_error("Error reading coordinates.");
                }
#else
                if (!GEOSCoordSeq_getX_r(context, seq, i, &x) || !GEOSCoordSeq_getY_r(context, seq, i, &y)) {
                    throw std::runtime_error("Error reading coordinates.");
                }
#endif

                if (is_ccw) {
                    stk.emplace_back(x, y);
                } else {
                    stk.emplace_front(x, y);
                }
            }
        }

        size_t row = ring_grid.get_row(stk.front().y);
        size_t col = ring_grid.get_column(stk.front().x);

        while (!stk.empty()) {
            Cell &cell = *get_cell(cells, ring_grid, row, col);

            while (!stk.empty()) {
                cell.take(stk.front());

                if (cell.last_traversal().exited()) {
                    // Only push our exit coordinate if it's not same as the
                    // coordinate we just took. This covers the case where
                    // the next coordinate in the stack falls exactly on
                    // the cell boundary.
                    const Coordinate &exc = cell.last_traversal().exit_coordinate();
                    if (exc != stk.front()) {
                        stk.emplace_front(exc.x, exc.y);
                    }
                    break;
                } else {
                    stk.pop_front();
                }
            }

            cell.force_exit();

            if (cell.last_traversal().exited()) {
                // When we start in the middle of a cell, we need to save the coordinates
                // from our incomplete traversal and reprocess them at the end of the line.
                // The effect is just to restructure the polygon so that the start/end
                // coordinate falls on a cell boundary.
                if (!cell.last_traversal().traversed()) {
                    for (const auto &coord : cell.last_traversal().coords()) {
                        stk.push_back(coord);
                    }
                }

                switch (cell.last_traversal().exit_side()) {
                    case Side::TOP:
                        row--;
                        break;
                    case Side::BOTTOM:
                        row++;
                        break;
                    case Side::LEFT:
                        col--;
                        break;
                    case Side::RIGHT:
                        col++;
                        break;
                    default:
                        throw std::runtime_error("Invalid traversal");
                }
            }
        }

        // Compute the fraction covered for all cells and assign it to
        // the area matrix
        // TODO avoid copying matrix when geometry has only one polygon, and polygon has only one ring
        Matrix<float> areas(rows - 2, cols - 2, fill_values<float>::FILLABLE);

        FloodFill ff(context, ls, make_finite(ring_grid));

        for (size_t i = 1; i <= areas.rows(); i++) {
            for (size_t j = 1; j <= areas.cols(); j++) {
                if (cells(i, j) != nullptr) {
                    // When we encounter a cell that has been processed (ie, it is not nullptr)
                    // but has zero covered fraction, we have no way to know if that cell is on
                    // the inside of the polygon. So we perform point-in-polygon test and set
                    // the covered fraction to 1.0 if needed.

                    auto frac = static_cast<float>(cells(i, j)->covered_fraction());
                    if (frac == 0) {
                        areas(i-1, j-1) = ff.cell_is_inside(i-1, j-1) ? fill_values<float>::INTERIOR : fill_values<float>::EXTERIOR;
                    } else {
                        areas(i-1, j-1) = frac;
                    }
                }
            }
        }

        ff.flood(areas);

        // Transfer these areas to our global set
        size_t i0 = ring_grid.row_offset(m_geometry_grid);
        size_t j0 = ring_grid.col_offset(m_geometry_grid);

        add_ring_areas(i0, j0, areas, exterior_ring);
    }

    void RasterCellIntersection::add_ring_areas(size_t i0, size_t j0, const Matrix<float> &areas, bool exterior_ring) {
        int factor = exterior_ring ? 1 : -1;

        for (size_t i = 0; i < areas.rows(); i++) {
            for (size_t j = 0; j < areas.cols(); j++) {
                m_overlap_areas->increment(i0 + i, j0 + j, factor * areas(i, j));
            }
        }
    }

}

// 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_RASTER_CELL_INTERSECTION_H
#define EXACTEXTRACT_RASTER_CELL_INTERSECTION_H

#include <memory>

#include <geos_c.h>

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

namespace exactextract {

    class RasterCellIntersection {

    public:
        RasterCellIntersection(const Grid<bounded_extent> &raster_grid, GEOSContextHandle_t context, const GEOSGeometry *g);

        RasterCellIntersection(const Grid<bounded_extent> &raster_grid, const Box & box);

        size_t rows() const { return m_overlap_areas->rows(); }

        size_t cols() const { return m_overlap_areas->cols(); }

        const Matrix<float> &overlap_areas() const { return *m_overlap_areas; }

        Grid<infinite_extent> m_geometry_grid;
    private:
        void process(GEOSContextHandle_t context, const GEOSGeometry *g);

        void process_ring(GEOSContextHandle_t context, const GEOSGeometry *ls, bool exterior_ring);

        void process_rectangular_ring(const Box & box, bool exterior_ring);

        void add_ring_areas(size_t i0, size_t j0, const Matrix<float> &areas, bool exterior_ring);

        std::unique_ptr<Matrix<float>> m_overlap_areas;

    };

    Raster<float> raster_cell_intersection(const Grid<bounded_extent> & raster_grid, GEOSContextHandle_t context, const GEOSGeometry* g);
    Raster<float> raster_cell_intersection(const Grid<bounded_extent> & raster_grid, const Box & box);
    Box processing_region(const Box & raster_extent, const std::vector<Box> & component_boxes);
}

#endif

// Copyright (c) 2018 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 <ostream>

#include "side.h"

namespace exactextract {

    std::ostream &operator<<(std::ostream &os, const Side &s) {
        switch (s) {
            case Side::NONE:
                os << "none";
                return os;
            case Side::LEFT:
                os << "left";
                return os;
            case Side::RIGHT:
                os << "right";
                return os;
            case Side::TOP:
                os << "top";
                return os;
            case Side::BOTTOM:
                os << "bottom";
                return os;
        }

        return os; // unreachable statement needed for -Werror=return-type
    }

}

// Copyright (c) 2018 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 <iostream>
#include <limits>
#include <vector>

#include "area.h"
#include "box.h"
#include "coordinate.h"
#include "perimeter_distance.h"

namespace exactextract {

    struct CoordinateChain {
        double start;
        double stop;
        const std::vector<Coordinate> *coordinates;
        bool visited;

        CoordinateChain(double p_start, double p_stop, const std::vector<Coordinate>* p_coords) :
                start{p_start},
                stop{p_stop},
                coordinates{p_coords},
                visited{false} {}
    };

    static double
    exit_to_entry_perimeter_distance_ccw(const CoordinateChain &c1, const CoordinateChain &c2, double perimeter) {
        return perimeter_distance_ccw(c1.stop, c2.start, perimeter);
    }

    static CoordinateChain *next_chain(std::vector<CoordinateChain> &chains,
                                       const CoordinateChain *chain,
                                       const CoordinateChain *kill,
                                       double perimeter) {

        CoordinateChain *min = nullptr;
        double min_distance = std::numeric_limits<double>::max();

        for (CoordinateChain &candidate : chains) {
            if (candidate.visited && std::addressof(candidate) != kill) {
                continue;
            }

            double distance = exit_to_entry_perimeter_distance_ccw(*chain, candidate, perimeter);
            if (distance < min_distance) {
                min_distance = distance;
                min = std::addressof(candidate);
            }
        }

        return min;
    }

    double left_hand_area(const Box &box, const std::vector<const std::vector<Coordinate> *> &coord_lists) {
        std::vector<CoordinateChain> chains;

        for (const auto &coords : coord_lists) {
            double start = perimeter_distance(box, (*coords)[0]);
            double stop = perimeter_distance(box, (*coords)[coords->size() - 1]);

            chains.emplace_back(start, stop, coords);
        }

        double height{box.height()};
        double width{box.width()};
        double perimeter{box.perimeter()};

        // create coordinate lists for corners
        std::vector<Coordinate> bottom_left = {Coordinate(box.xmin, box.ymin)};
        std::vector<Coordinate> top_left = {Coordinate(box.xmin, box.ymax)};
        std::vector<Coordinate> top_right = {Coordinate(box.xmax, box.ymax)};
        std::vector<Coordinate> bottom_right = {Coordinate(box.xmax, box.ymin)};

        // Add chains for corners
        chains.emplace_back(0.0, 0.0, &bottom_left);
        chains.emplace_back(height, height, &top_left);
        chains.emplace_back(height + width, height + width, &top_right);
        chains.emplace_back(2 * height + width, 2 * height + width, &bottom_right);

        double sum{0.0};
        for (auto &chain_ref : chains) {
            if (chain_ref.visited || chain_ref.coordinates->size() == 1) {
                continue;
            }

            std::vector<Coordinate> coords;
            CoordinateChain *chain = std::addressof(chain_ref);
            CoordinateChain *first_chain = chain;
            do {
                chain->visited = true;
                coords.insert(coords.end(), chain->coordinates->cbegin(), chain->coordinates->cend());
                chain = next_chain(chains, chain, first_chain, perimeter);
            } while (chain != first_chain);

            coords.push_back(coords[0]);

            sum += area(coords);
        }

        return sum;
    }

}

// Copyright (c) 2018 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 <cstddef>
#include <stdexcept>

#include "traversal.h"

namespace exactextract {

    void Traversal::add(const Coordinate &c) {
        m_coords.push_back(c);
    }

    bool Traversal::empty() const {
        return m_coords.empty();
    }

    void Traversal::enter(const Coordinate &c, Side s) {
        if (!m_coords.empty()) {
            throw std::runtime_error("Traversal already started");
        }

        add(c);
        m_entry = s;
    }

    void Traversal::exit(const Coordinate &c, Side s) {
        add(c);
        m_exit = s;
    }

    bool Traversal::is_closed_ring() const {
        return m_coords.size() >= 3 && m_coords[0] == m_coords[m_coords.size() - 1];
    }

    bool Traversal::entered() const {
        return m_entry != Side::NONE;
    }

    bool Traversal::exited() const {
        return m_exit != Side::NONE;
    }

    bool Traversal::multiple_unique_coordinates() const {
        for (size_t i = 1; i < m_coords.size(); i++) {
            if (m_coords[0] != m_coords[i]) {
                return true;
            }
        }

        return false;
    }

    bool Traversal::traversed() const {
        return entered() && exited();
    }

    const Coordinate &Traversal::last_coordinate() const {
        return m_coords.at(m_coords.size() - 1);
    }

    const Coordinate &Traversal::exit_coordinate() const {
        if (!exited()) {
            throw std::runtime_error("Can't get exit coordinate from incomplete traversal.");
        }

        return last_coordinate();
    }

}

// Copyright (c) 2019-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 "weighted_quantiles.h"

namespace exactextract {

    void WeightedQuantiles::prepare() const {
        std::sort(m_elems.begin(), m_elems.end(), [](const elem_t &a, const elem_t &b) {
            return a.x < b.x;
        });

        m_sum_w = 0;
        // relies on map being sorted which it is no
        for (size_t i = 0; i < m_elems.size(); i++) {
            m_sum_w += m_elems[i].w;

            if (i == 0) {
                m_elems[i].s = 0;
                m_elems[i].cumsum = m_elems[i].w;
            } else {
                m_elems[i].cumsum = m_elems[i - 1].cumsum + m_elems[i].w;
                m_elems[i].s = i * m_elems[i].w + (static_cast<double>(m_elems.size()) - 1) * m_elems[i - 1].cumsum;
            }
        }

        m_ready_to_query = true;
    }

    double WeightedQuantiles::quantile(double q) const {
        if (!std::isfinite(q) || q < 0 || q > 1) {
            throw std::runtime_error("Quantile must be between 0 and 1.");
        }

        if (!m_ready_to_query) {
            prepare();
        }

        auto sn = m_sum_w * (static_cast<double>(m_elems.size()) - 1);

        elem_t lookup(0, 0); // create a dummy element to use with std::upper_bound
        lookup.s = q*sn;

        // get first element that is greater than the lookup value (q * sn)
        auto right = std::upper_bound(m_elems.cbegin(), m_elems.cend(), lookup, [](const elem_t & a, const elem_t & b) {
            return a.s < b.s;
        });

        // since the minimum value of "lookup" is zero, and the first value of "s" is zero,
        // we are guaranteed to have at least one element to the left of "right"
        auto left = std::prev(right, 1);

        if (right == m_elems.cend()) {
            return left->x;
        }

        return left->x + (q*sn - left->s)*(right->x - left->x)/(right->s - left->s);
    }

}

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		#include <stdexcept>
15
16		#include "area.h"
17		#include "cell.h"
18		#include "crossing.h"
19		#include "traversal_areas.h"
20		#include "geos_utils.h"
21
22		namespace exactextract {
23
24	!	double Cell::height() const {
25	!	return m_box.height();
26		}
27
28	!	double Cell::width() const {
29	!	return m_box.width();
30		}
31
32	1604x	double Cell::area() const {
33	1604x	return m_box.area();
34		}
35
36	2172x	Side Cell::side(const Coordinate &c) const {
37	2172x	if (c.x == m_box.xmin) {
38	556x	return Side::LEFT;
39	1616x	} else if (c.x == m_box.xmax) {
40	653x	return Side::RIGHT;
41	963x	} else if (c.y == m_box.ymin) {
42	417x	return Side::BOTTOM;
43	546x	} else if (c.y == m_box.ymax) {
44	512x	return Side::TOP;
45		}
46
47	34x	return Side::NONE;
48		}
49
50	2055x	void Cell::force_exit() {
51	2055x	if (last_traversal().exited()) {
52	1938x	return;
53		}
54
55	117x	const Coordinate &last = last_traversal().last_coordinate();
56
57	117x	if (location(last) == Location::BOUNDARY) {
58	117x	last_traversal().force_exit(side(last));
59		}
60		}
61
62	15459x	Cell::Location Cell::location(const Coordinate &c) const {
63	15459x	if (m_box.strictly_contains(c)) {
64	13138x	return Cell::Location::INSIDE;
65		}
66
67	2321x	if (m_box.contains(c)) {
68	383x	return Cell::Location::BOUNDARY;
69		}
70
71	1938x	return Cell::Location::OUTSIDE;
72		}
73
74	17397x	Traversal &Cell::traversal_in_progress() {
75	17397x	if (m_traversals.empty() \|\| m_traversals[m_traversals.size() - 1].exited()) {
76	2055x	m_traversals.emplace_back();
77		}
78
79	17397x	return m_traversals[m_traversals.size() - 1];
80		}
81
82	27823x	Traversal &Cell::last_traversal() {
83	27823x	return m_traversals.at(m_traversals.size() - 1);
84		}
85
86	17397x	bool Cell::take(const Coordinate &c) {
87	17397x	Traversal &t = traversal_in_progress();
88
89	17397x	if (t.empty()) {
90		//std::cout << "Entering " << m_box << " from " << side(c) << " at " << c << std::endl;
91
92	2055x	t.enter(c, side(c));
93	2055x	return true;
94		}
95
96	15342x	if (location(c) != Cell::Location::OUTSIDE) {
97		//std::cout << "Still in " << m_box << " with " << c << std::endl;
98
99	13404x	t.add(c);
100	13404x	return true;
101		}
102
103	1938x	Crossing x = m_box.crossing(t.last_coordinate(), c);
104	1938x	t.exit(x.coord(), x.side());
105
106		//std::cout << "Leaving " << m_box << " from " << x.side() << " at " << x.coord();
107		//std::cout << " on the way to " << c << std::endl;
108
109	1938x	return false;
110		}
111
112	1604x	double Cell::covered_fraction() const {
113		// Handle the special case of a ring that is enclosed within a
114		// single pixel of our raster
115		if (m_traversals.size() == 1 && m_traversals[0].is_closed_ring()) {
116	!	return exactextract::area(m_traversals[0].coords()) / area();
117		}
118
119		// TODO consider porting in simplified single-traversal area calculations
120		// from Java code. Do they really make a performance difference?
121		//if (m_traversals.size() == 1) {
122		// double a = area();
123
124		// return (a - area_right_of(m_traversals.at(m_traversals.size() - 1))) / a;
125		//}
126
127	1604x	std::vector<const std::vector<Coordinate> *> coord_lists;
128
129	3301x	for (const auto &t : m_traversals) {
130	1697x	if (!t.traversed() \|\| !t.multiple_unique_coordinates()) {
131	86x	continue;
132		}
133
134	1611x	coord_lists.push_back(&t.coords());
135		}
136
137	1604x	return left_hand_area(m_box, coord_lists) / area();
138		}
139
140		#if 0
141		Crossing Cell::crossing(const Coordinate & c1, const Coordinate & c2) const {
142		Coordinate result(0, 0);
143
144		if (c2.y > c1.y && segment_intersection(c1, c2, m_box.upper_left(), m_box.upper_right(), result)) {
145		return Crossing(Side::TOP, result);
146		}
147
148		if (c2.y < c1.y && segment_intersection(c1, c2, m_box.lower_right(), m_box.lower_left(), result)) {
149		return Crossing(Side::BOTTOM, result);
150		}
151
152		if (c2.x < c1.x && segment_intersection(c1, c2, m_box.lower_left(), m_box.upper_left(), result)) {
153		return Crossing(Side::LEFT, result);
154		}
155
156		if (c2.x > c1.x && segment_intersection(c1, c2, m_box.lower_right(), m_box.upper_right(), result)) {
157		return Crossing(Side::RIGHT, result);
158		}
159
160		throw std::runtime_error("Never get here!");
161		}
162		#endif
163
164
165		}

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		if (!isGeneric("exact_extract")) {
15	136x	setGeneric("exact_extract", function(x, y, ...)
16	136x	standardGeneric("exact_extract"))
17		}
18
19		#' Extract or summarize values from Raster* objects
20		#'
21		#' Extracts the values of cells in a Raster* that are covered by a
22		#' simple feature collection containing polygonal geometries, as well as the
23		#' fraction of each cell that is covered by the polygon. Returns either
24		#' the result of a summary operation or function applied to the values
25		#' and coverage fractions (if \code{fun} is specified), or a data frame
26		#' containing the values and coverage fractions themselves (if \code{fun}
27		#' is \code{NULL}.)
28		#'
29		#' The value of \code{fun} may be set to a string (or vector of strings)
30		#' representing summary operations supported by the exactextract library.
31		#' If the input raster has a single layer and a single summary operation
32		#' is specified, \code{exact_extract} will return a vector with the result
33		#' of the summary operation for each feature in the input. If the input
34		#' raster has multiple layers, or if multiple summary operations are specified,
35		#' \code{exact_extract} will return a data frame with a row for each feature
36		#' and a column for each summary operation / layer combination. (The
37		#' \code{force_df} can be used to always return a data frame instead of a vector.)
38		#' In all of the summary operations, \code{NA} values in the raster are ignored
39		#' (i.e., \code{na.rm = TRUE}.)
40		#'
41		#' The following summary operations are supported:
42		#'
43		#' \itemize{
44		#' \item{\code{min} - the minimum defined value in any raster cell wholly or
45		#' partially covered by the polygon}
46		#' \item{\code{max} - the maximum defined value in any raster cell wholly or
47		#' partially covered by the polygon}
48		#' \item{\code{count} - the sum of fractions of raster cells with defined values
49		#' covered by the polygon}
50		#' \item{\code{sum} - the sum of defined raster cell values, multiplied by
51		#' the fraction of the cell that is covered by the polygon}
52		#' \item{\code{mean} - the mean cell value, weighted by the fraction of each cell
53		#' that is covered by the polygon}
54		#' \item{\code{median} - the median cell value, weighted by the fraction of each
55		#' cell that is covered by the polygon}
56		#' \item{\code{quantile} - arbitrary quantile(s) of cell values, specified in
57		#' \code{quantiles}, weighted by the fraction of each
58		#' cell that is covered by the polygon}
59		#' \item{\code{mode} - the most common cell value, weighted by the fraction of
60		#' each cell that is covered by the polygon. Where multiple
61		#' values occupy the same maximum number of weighted cells,
62		#' the largest value will be returned.}
63		#' \item{\code{majority} - synonym for \code{mode}}
64		#' \item{\code{minority} - the least common cell value, weighted by the fraction
65		#' of each cell that is covered by the polygon. Where
66		#' multiple values occupy the same minimum number of
67		#' weighted cells, the smallest value will be returned.}
68		#' \item{\code{variety} - the number of distinct values in cells that are wholly
69		#' or partially covered by the polygon.}
70		#' \item{\code{variance} - the population variance of cell values, weighted by the
71		#' fraction of each cell that is covered by the polygon.}
72		#' \item{\code{stdev} - the population standard deviation of cell values, weighted
73		#' by the fraction of each cell that is covered by the polygon.}
74		#' \item{\code{coefficient_of_variation} - the population coefficient of variation of
75		#' cell values, weighted by the fraction of each cell that is
76		#' covered by the polygon.}
77		#' \item{\code{weighted_mean} - the mean cell value, weighted by the product of
78		#' the fraction of each cell covered by the polygon
79		#' and the value of a second weighting raster provided
80		#' as \code{weights}}
81		#' \item{\code{weighted_sum} - the sum of defined raster cell values, multiplied by
82		#' the fraction of each cell that is covered by the polygon
83		#' and the value of a second weighting raster provided
84		#' as \code{weights}}
85		#' }
86		#'
87		#' Alternatively, an R function may be provided as \code{fun}. The function will be
88		#' called for each feature with with vectors of cell values and weights as arguments.
89		#' \code{exact_extract} will then return a vector of the return values of \code{fun}.
90		#'
91		#' If \code{fun} is not specified, \code{exact_extract} will return a list with
92		#' one data frame for each feature in the input feature collection. The data
93		#' frame will contain a column with values from each layer in the input `Raster*`,
94		#' and a final column indicating the fraction of the cell that is covered by the
95		#' polygon.
96		#'
97		#' @param x a \code{RasterLayer}, \code{RasterStack}, or \code{RasterBrick}
98		#' @param y a sf object with polygonal geometries
99		#' @param fun an optional function or character vector, as described below
100		#' @param weights a weighting raster to be used with the \code{weighted_mean}
101		#' and \code{weighted_sum} summary operations.
102		#' @param quantiles quantiles to be computed when \code{fun == 'quantile'}
103		#' @param append_cols when \code{fun} is not \code{NULL}, an optional
104		#' character vector of columns from \code{y} to be
105		#' included in returned data frame.
106		#' @param force_df always return a data frame instead of a vector, even if
107		#' \code{x} has only one layer and \code{fun} has length 1
108		#' @param full_colnames include the names of \code{x} in the names of the
109		#' returned data frame, even if \code{x} has only one
110		#' layer. This is useful when the results of multiple
111		#' calls to \code{exact_extract} are combined with
112		#' \code{cbind}.
113		#' @param include_cell if \code{TRUE}, and \code{fun} is \code{NULL}, augment
114		#' the returned data frame for each feature with a column
115		#' for the cell index (\code{cell}). If \code{TRUE} and
116		#' \code{fun} is not \code{NULL}, add \code{cell} to the
117		#' data frame passed to \code{fun} for each feature.
118		#' @param include_cols an optional character vector of column names in
119		#' \code{y} to be added to the data frame for each
120		#' feature that is either returned (when \code{fun} is
121		#' \code{NULL}) or passed to \code{fun}.
122		#' @param include_xy if \code{TRUE}, and \code{fun} is \code{NULL}, augment
123		#' the returned data frame for each feature with columns
124		#' for cell center coordinates (\code{x} and \code{y}). If
125		#' \code{TRUE} and \code{fun} is not \code{NULL}, add
126		#' \code{x} and {y} to the data frame passed to \code{fun}
127		#' for each feature.
128		#' @param stack_apply if \code{TRUE}, apply \code{fun} to each layer of
129		#' \code{x} independently. If \code{FALSE}, apply \code{fun}
130		#' to all layers of \code{x} simultaneously.
131		#' @param max_cells_in_memory the maximum number of raster cells to load at
132		#' a given time when using a named summary operation
133		#' for \code{fun} (as opposed to a function defined using
134		#' R code). If a polygon covers more than \code{max_cells_in_memory}
135		#' raster cells, it will be processed in multiple chunks.
136		#' @param progress if \code{TRUE}, display a progress bar during processing
137		#' @param ... additional arguments to pass to \code{fun}
138		#' @return a vector or list of data frames, depending on the type of \code{x} and the
139		#' value of \code{fun} (see Details)
140		#' @examples
141		#' rast <- raster::raster(matrix(1:100, ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
142		#' poly <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))')
143		#'
144		#' # named summary operation on RasterLayer, returns vector
145		#' exact_extract(rast, poly, 'mean')
146		#'
147		#' # two named summary operations on RasterLayer, returns data frame
148		#' exact_extract(rast, poly, c('min', 'max'))
149		#'
150		#' # named summary operation on RasterStack, returns data frame
151		#' stk <- raster::stack(list(a=rast, b=sqrt(rast)))
152		#' exact_extract(stk, poly, 'mean')
153		#'
154		#' # named weighted summary operation, returns vector
155		#' weights <- raster::raster(matrix(runif(100), ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
156		#' exact_extract(rast, poly, 'weighted_mean', weights=weights)
157		#'
158		#' # custom summary function, returns vector
159		#' exact_extract(rast, poly, function(value, cov_frac) length(value[cov_frac > 0.9]))
160		#'
161		#' @name exact_extract
162		NULL
163
164		#' @import sf
165		#' @import raster
166		#' @useDynLib exactextractr
167		#' @rdname exact_extract
168		#' @export
169		setMethod('exact_extract', signature(x='Raster', y='sf'),
170		function(x, y, fun=NULL, ...,
171		include_xy=FALSE,
172		progress=TRUE,
173		max_cells_in_memory=30000000,
174		include_cell=FALSE,
175		force_df=FALSE,
176		full_colnames=FALSE,
177		stack_apply=FALSE,
178		append_cols=NULL,
179		include_cols=NULL,
180		quantiles=NULL) {
181	9x	.exact_extract(x, y, fun=fun, ...,
182	9x	include_xy=include_xy,
183	9x	progress=progress,
184	9x	max_cells_in_memory=max_cells_in_memory,
185	9x	include_cell=include_cell,
186	9x	force_df=force_df,
187	9x	full_colnames=full_colnames,
188	9x	stack_apply=stack_apply,
189	9x	append_cols=append_cols,
190	9x	include_cols=include_cols,
191	9x	quantiles=quantiles)
192		})
193
194		# Return the number of standard (non-...) arguments in a supplied function that
195		# do not have a default value. This is used to fail if the summary function
196		# provided by the user cannot accept arguments of values and weights.
197		.num_expected_args <- function(fun) {
198	28x	a <- formals(args(fun))
199	28x	a <- a[names(a) != '...']
200	28x	sum(sapply(a, nchar) == 0)
201		}
202
203		emptyVector <- function(rast) {
204	12x	switch(substr(raster::dataType(rast), 1, 3),
205	!	LOG=logical(),
206	3x	INT=integer(),
207	9x	numeric())
208		}
209
210		.exact_extract <- function(x, y, fun=NULL, ...,
211		weights=NULL,
212		include_xy=FALSE,
213		progress=TRUE,
214		max_cells_in_memory=30000000,
215		include_cell=FALSE,
216		force_df=FALSE,
217		full_colnames=FALSE,
218		stack_apply=FALSE,
219		append_cols=NULL,
220		include_cols=NULL,
221		quantiles=NULL) {
222	131x	if(!is.null(append_cols)) {
223	2x	if (!inherits(y, 'sf')) {
224	!	stop(sprintf('append_cols only supported for sf arguments (received %s)',
225	!	paste(class(y), collapse = ' ')))
226		}
227
228	2x	force_df <- TRUE
229		}
230
231	131x	if(sf::st_geometry_type(y, by_geometry = FALSE) == 'GEOMETRY') {
232	2x	if (!all(sf::st_dimension(y) == 2)) {
233	1x	stop("Features in sfc_GEOMETRY must be polygonal")
234		}
235	1x	y <- sf::st_cast(y, 'MULTIPOLYGON')
236		}
237
238	130x	if(!is.null(weights)) {
239	28x	if (!startsWith(class(weights), 'Raster')) {
240	1x	stop("Weights must be a Raster object.")
241		}
242
243	27x	if (!is.character(fun)) {
244	1x	stop("Weighting raster can only be used with named summary operations.")
245		}
246
247	26x	if (!any(startsWith(fun, "weighted"))) {
248	9x	warning("Weights provided but no requested operations use them.")
249		}
250
251	26x	if (!is.na(sf::st_crs(x))) {
252	17x	if (is.na(sf::st_crs(weights))) {
253	1x	warning("No CRS specified for weighting raster; assuming it has the same CRS as the value raster.")
254	16x	} else if (sf::st_crs(x) != sf::st_crs(weights)) {
255	1x	stop("Weighting raster does not have the same CRS as value raster.")
256		}
257		}
258		}
259
260	127x	if(is.na(sf::st_crs(x)) && !is.na(sf::st_crs(y))) {
261	1x	warning("No CRS specified for raster; assuming it has the same CRS as the polygons.")
262	126x	} else if(is.na(sf::st_crs(y)) && !is.na(sf::st_crs(x))) {
263	2x	warning("No CRS specified for polygons; assuming they have the same CRS as the raster.")
264	124x	} else if(sf::st_crs(x) != sf::st_crs(y)) {
265	1x	y <- sf::st_transform(y, sf::st_crs(x))
266	1x	warning("Polygons transformed to raster CRS (EPSG:", sf::st_crs(x)$epsg, ")")
267		}
268
269	127x	if (!is.null(fun) && !is.character(fun) && .num_expected_args(fun) < 2) {
270	5x	stop("exact_extract was called with a function that does not appear to ",
271	5x	"be of the form `function(values, coverage_fractions, ...)`")
272		}
273
274	122x	if (is.character(fun)) {
275	87x	if (length(list(...)) > 0) {
276	1x	stop("exact_extract was called with a named summary operation that ",
277	1x	"does not accept additional arguments ...")
278		}
279
280	86x	if (include_xy) {
281	1x	stop("include_xy must be FALSE for named summary operations")
282		}
283
284	85x	if (include_cell) {
285	1x	stop("include_cell must be FALSE for named summary operations")
286		}
287
288	84x	if (!is.null(include_cols)) {
289	1x	stop("include_cols not supported for named_summary operations (see argument append_cols)")
290		}
291		}
292
293	118x	if (progress && length(y) > 1) {
294	1x	n <- length(y)
295	1x	pb <- utils::txtProgressBar(min = 0, max = n, initial=0, style=3)
296	1x	update_progress <- function() {
297	2x	i <- 1 + utils::getTxtProgressBar(pb)
298	2x	utils::setTxtProgressBar(pb, i)
299	2x	if (i == n) {
300	!	close(pb)
301		}
302		}
303		} else {
304	117x	update_progress <- function() {}
305		}
306
307	118x	tryCatch({
308	118x	x <- raster::readStart(x)
309	118x	if (!is.null(weights)) {
310	25x	weights <- raster::readStart(weights)
311		}
312
313	118x	if (is.character(fun)) {
314		# Compute all stats in C++
315	83x	results <- sapply(sf::st_as_binary(sf::st_geometry(y), EWKB=TRUE), function(wkb) {
316	91x	ret <- CPP_stats(x, weights, wkb, fun, max_cells_in_memory, quantiles)
317	81x	update_progress()
318	81x	return(ret)
319		})
320
321	73x	if (length(fun) == 1 &&
322	73x	raster::nlayers(x) == 1 &&
323	73x	(is.null(quantiles) \|\| length(quantiles) == 1) &&
324	73x	!force_df) {
325		# Just return a vector of stat results
326	60x	return(as.vector(results))
327		} else {
328		# Return a data frame with a column for each stat
329	13x	colnames <- .cppStatColNames(x, fun, full_colnames, quantiles)
330
331	13x	if (is.matrix(results)) {
332	11x	results <- t(results)
333		} else {
334	2x	results <- matrix(results, nrow=length(results))
335		}
336
337	13x	dimnames(results) <- list(NULL, colnames)
338	13x	ret <- as.data.frame(results)
339
340	13x	if (!is.null(append_cols)) {
341	1x	ret <- cbind(sf::st_drop_geometry(y[, append_cols]), ret)
342		}
343
344	13x	return(ret)
345		}
346		} else {
347	35x	geoms <- sf::st_geometry(y)
348
349	35x	ret <- lapply(seq_along(geoms), function(feature_num) {
350	45x	wkb <- sf::st_as_binary(geoms[[feature_num]], EWKB=TRUE)
351
352	45x	ret <- CPP_exact_extract(x, wkb)
353
354	45x	if (length(ret$weights) > 0) {
355	39x	vals <- raster::getValuesBlock(x,
356	39x	row=ret$row,
357	39x	col=ret$col,
358	39x	nrow=nrow(ret$weights),
359	39x	ncol=ncol(ret$weights))
360
361	39x	if(is.matrix(vals)) {
362	7x	vals <- as.data.frame(vals)
363		} else {
364	32x	vals <- data.frame(value=vals)
365		}
366		} else {
367		# Polygon does not intersect raster.
368		# Construct a zero-row data frame with correct column names/types.
369	6x	vals <- do.call(data.frame, lapply(seq_len(raster::nlayers(x)),
370	6x	function(i) emptyVector(x[[i]])))
371	6x	if (raster::nlayers(x) == 1) {
372	3x	names(vals) <- 'value'
373		} else {
374	3x	names(vals) <- names(x)
375		}
376		}
377
378	45x	if (include_xy) {
379	10x	vals <- .appendXY(vals, x, ret$row, nrow(ret$weights), ret$col, ncol(ret$weights))
380		}
381
382	45x	if (include_cell) {
383	4x	vals <- .appendCell(vals, x, ret$row, nrow(ret$weights), ret$col, ncol(ret$weights))
384		}
385
386	45x	if (!is.null(include_cols)) {
387		# use vals as first argument to cbind, then rearrange names so that
388		# include_cols come first
389	2x	vals <- cbind(vals,
390	2x	sf::st_drop_geometry(y[feature_num, include_cols]),
391	2x	row.names = NULL)[, c(include_cols, names(vals))]
392		}
393
394	45x	cov_fracs <- as.vector(t(ret$weights))
395	45x	vals <- vals[cov_fracs > 0, , drop=FALSE]
396	45x	cov_fracs <- cov_fracs[cov_fracs > 0]
397
398	45x	update_progress()
399
400	45x	if (is.null(fun)) {
401	13x	vals$coverage_fraction <- cov_fracs
402	13x	return(vals)
403		} else {
404	32x	if (ncol(vals) == 1) {
405		# Only one layer, nothing appended (cells or XY)
406	24x	return(fun(vals[,1], cov_fracs, ...))
407		} else {
408	8x	if (stack_apply) {
409		# Pass each layer in stack to callback individually
410	4x	nlay <- raster::nlayers(x)
411	4x	appended_cols <- seq_len(ncol(vals))[-seq_len(nlay)]
412
413	4x	if (length(appended_cols) == 0) {
414	2x	result <- lapply(seq_len(nlay), function(z)
415	2x	fun(vals[, z], cov_fracs, ...))
416		} else {
417	2x	result <- lapply(seq_len(nlay), function(z)
418	2x	fun(cbind(data.frame(value=vals[, z]),
419	2x	vals[, c(appended_cols)]), cov_fracs, ...))
420		}
421
422	4x	names(result) <- paste('fun', names(x), sep='.')
423	4x	return(do.call(data.frame, result))
424		} else {
425		# Pass all layers to callback, to be handled together
426	4x	return(fun(vals, cov_fracs, ...))
427		}
428		}
429		}
430		})
431
432	34x	if (!is.null(fun)) {
433	22x	if (class(ret[[1]]) == 'data.frame') {
434	4x	if (requireNamespace('dplyr', quietly = TRUE)) {
435	4x	ret <- dplyr::bind_rows(ret) # handle column name mismatches
436		} else {
437	!	ret <- do.call(rbind, ret)
438		}
439		} else {
440	18x	ret <- simplify2array(ret)
441
442	18x	if (force_df) {
443	3x	ret <- data.frame(result = ret)
444		}
445		}
446		}
447
448	34x	if (!is.null(append_cols)) {
449	1x	ret <- cbind(sf::st_drop_geometry(y[, append_cols]), ret)
450		}
451
452	34x	return(ret)
453		}
454	118x	}, finally={
455	118x	raster::readStop(x)
456	118x	if (!is.null(weights)) {
457	25x	raster::readStop(weights)
458		}
459		})
460		}
461
462		.appendXY <- function(vals_df, rast, first_row, nrow, first_col, ncol) {
463	10x	if (nrow(vals_df) == 0) {
464	2x	vals_df$x <- numeric()
465	2x	vals_df$y <- numeric()
466		} else {
467	8x	x_coords <- raster::xFromCol(rast, col=seq(first_col, first_col + ncol - 1))
468	8x	y_coords <- raster::yFromRow(rast, row=seq(first_row, first_row + nrow - 1))
469
470	8x	vals_df$x <- rep(x_coords, times=nrow)
471	8x	vals_df$y <- rep(y_coords, each=ncol)
472		}
473
474	10x	return(vals_df)
475		}
476
477		.appendCell <- function(vals_df, rast, first_row, nrow, first_col, ncol) {
478	4x	if (nrow(vals_df) == 0) {
479	2x	vals_df$cell <- numeric()
480		} else {
481	2x	rows <- rep(seq(first_row, first_row + nrow - 1), each = ncol)
482	2x	cols <- rep(seq(first_col, first_col + ncol - 1), times = nrow)
483
484	2x	vals_df$cell <- raster::cellFromRowCol(rast, row=rows, col=cols)
485		}
486
487	4x	return(vals_df)
488		}
489
490		.cppStatColNames <- function(rast, stat_names, full_colnames, quantiles) {
491	13x	quantile_index = which(stat_names == 'quantile')
492	13x	if (length(quantile_index) != 0) {
493	1x	stat_names <- c(stat_names[seq_along(stat_names) < quantile_index],
494	1x	sprintf('q%02d', as.integer(100 * quantiles)),
495	1x	stat_names[seq_along(stat_names) > quantile_index])
496		}
497
498	13x	if (raster::nlayers(rast) > 1 \|\| full_colnames) {
499	6x	z <- expand.grid(names(rast), stat_names, stringsAsFactors=TRUE)
500	6x	mapply(paste, z[[2]], z[[1]], MoreArgs=list(sep='.'))
501		} else {
502	7x	stat_names
503		}
504		}
505
506		#' @useDynLib exactextractr
507		#' @rdname exact_extract
508		#' @export
509		setMethod('exact_extract', signature(x='Raster', y='sfc_MULTIPOLYGON'), .exact_extract)
510
511		#' @useDynLib exactextractr
512		#' @rdname exact_extract
513		#' @export
514		setMethod('exact_extract', signature(x='Raster', y='sfc_POLYGON'), .exact_extract)
515
516		#' @useDynLib exactextractr
517		#' @rdname exact_extract
518		#' @export
519		setMethod('exact_extract', signature(x='Raster', y='sfc_GEOMETRY'), .exact_extract)

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		if (!isGeneric("exact_resample")) {
15	6x	setGeneric("exact_resample", function(x, y, ...)
16	6x	standardGeneric("exact_resample"))
17		}
18
19		#' Resample a raster to a new grid
20		#'
21		#' @param x a \code{RasterLayer} to be resampled
22		#' @param y a \code{RasterLayer} with a grid definition to which \code{x}
23		#' should be resampled
24		#' @param fun a named summary operation to be used for the resampling
25		#' @return a resampled version of \code{x}
26		#'
27		#' @name exact_resample
28		NULL
29
30		#' @import sf
31		#' @import raster
32		#' @useDynLib exactextractr
33		#' @rdname exact_resample
34		#' @export
35		setMethod('exact_resample',
36		signature(x='RasterLayer', y='RasterLayer'),
37		function(x, y, fun) {
38	6x	if (sf::st_crs(x) != sf::st_crs(y)) {
39	3x	if (is.na(sf::st_crs(x))) {
40	1x	warning("No CRS specified for source raster; assuming it has the same CRS as destination raster.")
41	2x	} else if (is.na(sf::st_crs(y))) {
42	1x	warning("No CRS specified for destination raster; assuming it has the same CRS as source raster.")
43		} else {
44	1x	stop('Destination raster must have same CRS as source.')
45		}
46		}
47
48	5x	x <- raster::readStart(x)
49	5x	tryCatch({
50	5x	CPP_resample(x, y, fun)
51	5x	}, finally={
52	5x	raster::readStop(x)
53		})
54		})

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		// [[Rcpp::plugins("cpp14")]]
15		#include <memory>
16
17		#include <Rcpp.h>
18		#include <geos_c.h>
19
20		#include "exactextract/src/geos_utils.h"
21		#include "exactextract/src/grid.h"
22		#include "exactextract/src/matrix.h"
23		#include "exactextract/src/raster_cell_intersection.h"
24		#include "exactextract/src/raster_source.h"
25		#include "exactextract/src/raster_stats.h"
26
27		using geom_ptr= std::unique_ptr<GEOSGeometry, std::function<void(GEOSGeometry*)>>;
28		using wkb_reader_ptr = std::unique_ptr<GEOSWKBReader, std::function<void(GEOSWKBReader*)>>;
29
30		using exactextract::Box;
31		using exactextract::Matrix;
32		using exactextract::Grid;
33		using exactextract::subdivide;
34		using exactextract::bounded_extent;
35		using exactextract::Raster;
36		using exactextract::RasterView;
37		using exactextract::raster_cell_intersection;
38		using exactextract::RasterStats;
39		using exactextract::RasterSource;
40
41	183x	static Grid<bounded_extent> make_grid(const Rcpp::S4 & rast) {
42	549x	Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster");
43
44	366x	Rcpp::S4 extent = rast.slot("extent");
45
46	549x	Rcpp::Function resFn = raster["res"];
47
48	183x	Rcpp::NumericVector res = resFn(rast);
49
50		return {{
51	366x	extent.slot("xmin"),
52	366x	extent.slot("ymin"),
53	366x	extent.slot("xmax"),
54	366x	extent.slot("ymax"),
55		},
56	183x	res[0],
57	366x	res[1]
58		};
59		}
60
61		// Construct a Raster using an R vector for storage
62		class NumericVectorRaster : public exactextract::AbstractRaster<double> {
63		public:
64	663x	NumericVectorRaster(const Rcpp::NumericVector & vec, const Grid<bounded_extent> & g) :
65		AbstractRaster<double>(g),
66	663x	m_vec(vec)
67		{}
68
69	538823x	double operator()(size_t row, size_t col) const final {
70	538823x	return m_vec[row*cols() + col];
71		}
72
73		private:
74		const Rcpp::NumericVector m_vec;
75		};
76
77		// Read raster values from an R raster object
78		class S4RasterSource {
79		public:
80	118x	S4RasterSource(Rcpp::S4 rast) :
81	118x	m_grid(Grid<bounded_extent>::make_empty()),
82		m_rast(rast),
83	118x	m_last_box(0, 0, 0, 0)
84		{
85	118x	m_grid = make_grid(rast);
86		}
87
88	109x	const Grid<bounded_extent> &grid() const {
89	109x	return m_grid;
90		}
91
92	663x	std::unique_ptr<exactextract::AbstractRaster<double>> read_box(const exactextract::Box & box, int layer) {
93	663x	auto cropped_grid = m_grid.crop(box);
94
95	663x	if (!(box == m_last_box)) {
96	658x	m_last_box = box;
97
98	1974x	Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster");
99	1974x	Rcpp::Function getValuesBlockFn = raster["getValuesBlock"];
100
101	658x	if (cropped_grid.empty()) {
102	77x	m_rast_values = Rcpp::no_init(0, 0);
103		} else {
104		// Instead of reading only values for the requested band, we read values
105		// for all requested bands and then cache them to return from subsequent
106		// calls to read_box. There are two reasons for this:
107		// 1) Use of the 'lyrs` argument to getValuesBlock does not work for
108		// a single band in the CRAN version of raster
109		// 2) Reading everything once avoids the very significant per-call
110		// overhead of getValuesBlock, which largely comes from operations
111		// like setting names on the result vector.
112	581x	m_rast_values = getValuesBlockFn(m_rast,
113	!	1 + cropped_grid.row_offset(m_grid),
114	581x	cropped_grid.rows(),
115	581x	1 + cropped_grid.col_offset(m_grid),
116	581x	cropped_grid.cols(),
117	1162x	Rcpp::Named("format", "m"));
118		}
119		}
120
121	663x	if (cropped_grid.empty()) {
122	77x	return std::make_unique<NumericVectorRaster>(m_rast_values, cropped_grid);
123		} else {
124	1172x	return std::make_unique<NumericVectorRaster>(m_rast_values(Rcpp::_, layer),
125	586x	cropped_grid);
126		}
127		}
128
129		private:
130		Grid<bounded_extent> m_grid;
131		Rcpp::S4 m_rast;
132		Rcpp::NumericMatrix m_rast_values;
133		exactextract::Box m_last_box;
134		};
135
136		// GEOS warning handler
137	!	static void geos_warn(const char* fmt, ...) {
138	!	char buf[BUFSIZ] = { '\0' };
139
140		va_list msg;
141	!	va_start(msg, fmt);
142	!	vsnprintf(buf, BUFSIZ*sizeof(char), fmt, msg);
143	!	va_end(msg);
144
145	!	Rcpp::Function warning("warning");
146	!	warning(buf);
147		}
148
149		// GEOS error handler
150	!	static void geos_error(const char* fmt, ...) {
151	!	char buf[BUFSIZ] = { '\0' };
152
153		va_list msg;
154	!	va_start(msg, fmt);
155	!	vsnprintf(buf, BUFSIZ*sizeof(char), fmt, msg);
156	!	va_end(msg);
157
158	!	Rcpp::stop(buf);
159		}
160
161		// GEOSContextHandle wrapper to ensure finishGEOS is called.
162		struct GEOSAutoHandle {
163	148x	GEOSAutoHandle() {
164	148x	handle = initGEOS_r(geos_warn, geos_error);
165		}
166
167	296x	~GEOSAutoHandle() {
168	148x	finishGEOS_r(handle);
169		}
170
171		GEOSContextHandle_t handle;
172		};
173
174		// Return a smart pointer to a Geometry, given WKB input
175	142x	static geom_ptr read_wkb(const GEOSContextHandle_t & context, const Rcpp::RawVector & wkb) {
176	426x	wkb_reader_ptr wkb_reader{ GEOSWKBReader_create_r(context), [context](GEOSWKBReader* r) { GEOSWKBReader_destroy_r(context, r); } };
177
178		geom_ptr geom{ GEOSWKBReader_read_r(context,
179		wkb_reader.get(),
180	284x	std::addressof(wkb[0]),
181	142x	wkb.size()),
182	284x	[context](GEOSGeometry* g) { GEOSGeom_destroy_r(context, g); } };
183
184	142x	if (geom.get() == nullptr) {
185	!	Rcpp::stop("Failed to parse WKB geometry");
186		}
187
188	284x	return geom;
189		}
190
191		// [[Rcpp::export]]
192	45x	Rcpp::List CPP_exact_extract(Rcpp::S4 & rast, const Rcpp::RawVector & wkb) {
193	90x	GEOSAutoHandle geos;
194
195	45x	auto grid = make_grid(rast);
196	90x	auto coverage_fractions = raster_cell_intersection(grid, geos.handle, read_wkb(geos.handle, wkb).get());
197
198	45x	size_t nrow = coverage_fractions.rows();
199	45x	size_t ncol = coverage_fractions.cols();
200
201	45x	Rcpp::NumericMatrix weights = Rcpp::no_init(nrow, ncol);
202	159x	for (size_t i = 0; i < nrow; i++) {
203	507x	for (size_t j = 0; j < ncol; j++) {
204	393x	weights(i, j) = coverage_fractions(i ,j);
205		}
206		}
207
208	45x	if (nrow > 0) {
209	39x	size_t row_us = 1 + coverage_fractions.grid().row_offset(grid);
210	39x	if (row_us > static_cast<size_t>(std::numeric_limits<int>::max())) {
211
212		}
213		}
214
215	45x	int row = NA_INTEGER;
216	45x	int col = NA_INTEGER;
217	45x	if (nrow > 0) {
218	39x	size_t row_us = (1 + coverage_fractions.grid().row_offset(grid));
219	39x	if (row_us > static_cast<size_t>(std::numeric_limits<int>::max())) {
220	!	throw std::runtime_error("Cannot represent row offset as an R integer");
221		}
222	39x	row = static_cast<int>(row_us);
223		}
224	45x	if (ncol > 0) {
225	39x	size_t col_us = (1 + coverage_fractions.grid().col_offset(grid));
226	39x	if (col_us > static_cast<size_t>(std::numeric_limits<int>::max())) {
227	!	throw std::runtime_error("Cannot represent column offset as an R integer");
228		}
229	39x	col = static_cast<int>(col_us);
230		}
231
232		return Rcpp::List::create(
233	90x	Rcpp::Named("row") = row,
234	90x	Rcpp::Named("col") = col,
235	90x	Rcpp::Named("weights") = weights
236		);
237		}
238
239		// [[Rcpp::export]]
240	12x	Rcpp::S4 CPP_coverage_fraction(Rcpp::S4 & rast, const Rcpp::RawVector & wkb, bool crop)
241		{
242	24x	GEOSAutoHandle geos;
243	36x	Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster");
244	36x	Rcpp::Function rasterFn = raster["raster"];
245	36x	Rcpp::Function crsFn = raster["crs"];
246
247	12x	auto grid = make_grid(rast);
248	24x	auto coverage_fraction = raster_cell_intersection(grid, geos.handle, read_wkb(geos.handle, wkb).get());
249
250	12x	if (crop) {
251	2x	grid = coverage_fraction.grid();
252		}
253	24x	RasterView<float> coverage_view(coverage_fraction, grid);
254
255	12x	Rcpp::NumericMatrix weights{static_cast<int>(grid.rows()),
256	12x	static_cast<int>(grid.cols())};
257
258	629x	for (size_t i = 0; i < grid.rows(); i++) {
259	159018x	for (size_t j = 0; j < grid.cols(); j++) {
260	158401x	weights(i, j) = coverage_view(i, j);
261	158401x	if (!crop && std::isnan(weights(i, j))) {
262	139692x	weights(i, j) = 0;
263		}
264		}
265		}
266
267		return rasterFn(weights,
268	24x	Rcpp::Named("xmn")=grid.xmin(),
269	24x	Rcpp::Named("xmx")=grid.xmax(),
270	24x	Rcpp::Named("ymn")=grid.ymin(),
271	24x	Rcpp::Named("ymx")=grid.ymax(),
272	48x	Rcpp::Named("crs")=crsFn(rast));
273		}
274
275	115x	static int get_nlayers(Rcpp::S4 & rast) {
276	345x	Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster");
277	345x	Rcpp::Function nlayersFn = raster["nlayers"];
278
279	115x	Rcpp::NumericVector nlayersVec = nlayersFn(rast);
280
281	230x	return static_cast<int>(nlayersVec[0]);
282		}
283
284	500250x	static double get_stat_value(const RasterStats<double> & stats, const std::string & stat_name) {
285		if (stat_name == "mean") return stats.mean();
286
287	500250x	else if (stat_name == "sum") return stats.sum();
288	!	else if (stat_name == "count") return stats.count();
289
290	!	else if (stat_name == "min") return stats.min().value_or(NA_REAL);
291	!	else if (stat_name == "max") return stats.max().value_or(NA_REAL);
292
293	!	else if (stat_name == "mode") return stats.mode().value_or(NA_REAL);
294	!	else if (stat_name == "majority") return stats.mode().value_or(NA_REAL);
295	!	else if (stat_name == "minority") return stats.minority().value_or(NA_REAL);
296
297	!	else if (stat_name == "variety") return stats.variety();
298	!	else if (stat_name == "weighted_mean") return stats.weighted_mean();
299	!	else if (stat_name == "weighted_sum") return stats.weighted_sum();
300
301	!	else if (stat_name == "variance") return stats.variance();
302	!	else if (stat_name == "stdev") return stats.stdev();
303	!	else if (stat_name == "coefficient_of_variation") return stats.coefficient_of_variation();
304
305	!	else Rcpp::stop("Unknown stat: " + stat_name);
306		}
307
308
309		// [[Rcpp::export]]
310	91x	Rcpp::NumericMatrix CPP_stats(Rcpp::S4 & rast,
311		Rcpp::Nullable<Rcpp::S4> weights,
312		const Rcpp::RawVector & wkb,
313		const Rcpp::StringVector & stats,
314		int max_cells_in_memory,
315		const Rcpp::Nullable<Rcpp::NumericVector> & quantiles) {
316		try {
317	182x	GEOSAutoHandle geos;
318
319	91x	if (max_cells_in_memory < 1) {
320	1x	Rcpp::stop("Invalid value for max_cells_in_memory: %d", max_cells_in_memory);
321		}
322
323	90x	int nlayers = get_nlayers(rast);
324
325	180x	S4RasterSource rsrc(rast);
326
327	90x	std::unique_ptr<S4RasterSource> rweights;
328	90x	bool weighted = false;
329	90x	if (weights.isNotNull()) {
330	26x	Rcpp::S4 weights_s4 = weights.get();
331
332	25x	if (get_nlayers(weights_s4) != 1) {
333	1x	Rcpp::stop("Weighting raster must have only a single layer.");
334		}
335
336	24x	rweights = std::make_unique<S4RasterSource>(weights_s4);
337	24x	weighted = true;
338		}
339
340	89x	bool store_values = false;
341	89x	int stat_result_size = 0;
342	187x	for (const auto & stat : stats) {
343		// explicit construction of std::string seems necessary to avoid ambiguous overload error
344	204x	if (stat == std::string("mode") \|\|
345	199x	stat == std::string("majority") \|\|
346	196x	stat == std::string("minority") \|\|
347	192x	stat == std::string("variety") \|\|
348	484x	stat == std::string("median") \|\|
349	183x	stat == std::string("quantile")) {
350	29x	store_values = true;
351		}
352
353	179x	if (!weighted &&
354	179x	(stat == std::string("weighted_mean") \|\|
355	178x	stat == std::string("weighted_sum"))) {
356	2x	Rcpp::stop("Weighted stat requested but no weights provided.");
357		}
358
359	100x	if (stat == std::string("quantile")) {
360	8x	int num_quantiles = 0;
361
362	8x	if (quantiles.isNotNull()) {
363	7x	Rcpp::NumericVector qvec = quantiles.get();
364	7x	num_quantiles = qvec.size();
365		}
366
367	8x	if (num_quantiles == 0) {
368	2x	Rcpp::stop("Quantiles not specified.");
369		}
370
371	6x	stat_result_size += num_quantiles;
372		} else {
373	92x	stat_result_size += 1;
374		}
375		}
376
377	170x	std::vector<RasterStats<double>> raster_stats;
378	85x	raster_stats.reserve(nlayers);
379	175x	for (int i = 0; i < nlayers; i++) {
380	90x	raster_stats.emplace_back(store_values);
381		}
382
383	170x	Rcpp::NumericMatrix stat_results = Rcpp::no_init(nlayers, stat_result_size);
384
385	170x	auto geom = read_wkb(geos.handle, wkb);
386	85x	auto bbox = exactextract::geos_get_box(geos.handle, geom.get());
387
388	85x	auto grid = weighted ? rsrc.grid().common_grid(rweights->grid()) : rsrc.grid();
389
390	84x	if (bbox.intersects(grid.extent())) {
391	75x	auto cropped_grid = grid.crop(bbox);
392
393	185x	for (const auto &subgrid : subdivide(cropped_grid, max_cells_in_memory)) {
394	220x	auto coverage_fraction = raster_cell_intersection(subgrid, geos.handle, geom.get());
395	110x	auto& cov_grid = coverage_fraction.grid();
396
397	110x	if (!cov_grid.empty()) {
398	110x	if (weighted) {
399	46x	auto weights = rweights->read_box(cov_grid.extent(), 0);
400
401	47x	for (int i = 0; i < nlayers; i++) {
402	48x	auto values = rsrc.read_box(cov_grid.extent(), i);
403	24x	raster_stats[i].process(coverage_fraction, values, weights);
404		}
405		} else {
406	178x	for (int i = 0; i < nlayers; i++) {
407	182x	auto values = rsrc.read_box(cov_grid.extent(), i);
408	91x	raster_stats[i].process(coverage_fraction, *values);
409		}
410		}
411		}
412		}
413		}
414
415	170x	for (int j = 0; j < nlayers; j++) {
416	89x	const auto& rs = raster_stats[j];
417
418	89x	int i = 0;
419	190x	for(const auto& stat : stats) {
420	104x	if (stat == std::string("mean")) stat_results(j, i++) = rs.mean();
421
422	83x	else if (stat == std::string("sum")) stat_results(j, i++) = rs.sum();
423	68x	else if (stat == std::string("count")) stat_results(j, i++) = rs.count();
424
425	62x	else if (stat == std::string("min")) stat_results(j, i++) = rs.min().value_or(NA_REAL);
426	56x	else if (stat == std::string("max")) stat_results(j, i++) = rs.max().value_or(NA_REAL);
427
428	50x	else if (stat == std::string("median")) stat_results(j, i++) = rs.quantile(0.5).value_or(NA_REAL);
429
430	49x	else if (stat == std::string("mode")) stat_results(j, i++) = rs.mode().value_or(NA_REAL);
431	44x	else if (stat == std::string("majority")) stat_results(j, i++) = rs.mode().value_or(NA_REAL);
432	41x	else if (stat == std::string("minority")) stat_results(j, i++) = rs.minority().value_or(NA_REAL);
433
434	37x	else if (stat == std::string("variety")) stat_results(j, i++) = rs.variety();
435	25x	else if (stat == std::string("weighted_mean")) stat_results(j, i++) = rs.weighted_mean();
436	13x	else if (stat == std::string("weighted_sum")) stat_results(j, i++) = rs.weighted_sum();
437
438	10x	else if (stat == std::string("variance")) stat_results(j, i++) = rs.variance();
439	9x	else if (stat == std::string("stdev")) stat_results(j, i++) = rs.stdev();
440	8x	else if (stat == std::string("coefficient_of_variation")) stat_results(j, i++) = rs.coefficient_of_variation();
441
442	7x	else if (stat == std::string("quantile")) {
443	12x	Rcpp::NumericVector qvec = quantiles.get();
444	15x	for (double q : qvec) {
445	11x	stat_results(j, i++) = rs.quantile(q).value_or(NA_REAL);
446		}
447		}
448
449	4x	else Rcpp::stop("Unknown stat: " + stat);
450		}
451		}
452
453	162x	return stat_results;
454	20x	} catch (std::exception & e) {
455		// throw predictable exception class
456	10x	Rcpp::stop(e.what());
457		}
458		}
459
460		// [[Rcpp::export]]
461	5x	Rcpp::S4 CPP_resample(Rcpp::S4 & rast_in,
462		Rcpp::S4 & rast_out,
463		const Rcpp::StringVector & stat) {
464	15x	Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster");
465	15x	Rcpp::Function rasterFn = raster["raster"];
466	15x	Rcpp::Function valuesFn = raster["values<-"];
467
468	5x	if (stat.size() != 1) {
469	1x	Rcpp::stop("Only a single operation may be used for resampling.");
470		}
471
472	8x	S4RasterSource rsrc(rast_in);
473
474	4x	Rcpp::S4 out = rasterFn(rast_out);
475
476	4x	auto grid_in = make_grid(rast_in);
477	4x	auto grid_out = make_grid(rast_out);
478
479	8x	std::string stat_name = Rcpp::as<std::string>(stat[0]);
480
481	8x	Rcpp::NumericMatrix values_out = Rcpp::no_init(grid_out.rows(), grid_out.cols());
482
483	529x	for (size_t row = 0; row < grid_out.rows(); row++) {
484		// Read enough source raster data to process an entire destination row at
485		// a time, since getValuesBlock calls have a lot of overhead.
486	525x	auto y = grid_out.y_for_row(row);
487	525x	auto ymin = y - grid_out.dy();
488	525x	auto ymax = y + grid_out.dy();
489
490	525x	Box row_box{ grid_out.xmin(), ymin, grid_out.xmax(), ymax };
491	1050x	auto values = rsrc.read_box(row_box, 0);
492
493	500775x	for (size_t col = 0; col < grid_out.cols(); col++) {
494	1000500x	RasterStats<double> stats;
495
496	500250x	Box cell = grid_cell(grid_out, row, col);
497	500250x	auto coverage_fraction = raster_cell_intersection(grid_in, cell);
498
499	500250x	auto& cov_grid = coverage_fraction.grid();
500
501	500250x	if (!cov_grid.empty()) {
502	382698x	stats.process(coverage_fraction, *values);
503		}
504
505	500250x	values_out(row, col) = get_stat_value(stats, stat_name);
506		}
507		}
508
509	4x	out = valuesFn(out, values_out);
510	8x	return out;
511		}

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	8379714x	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	8379714x	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	167x	static Box make_empty() {
47	167x	return {0, 0, 0, 0};
48		}
49
50	6112695x	double width() const {
51	6112695x	return xmax - xmin;
52		}
53
54	6112695x	double height() const {
55	6112695x	return ymax - ymin;
56		}
57
58	6109487x	double area() const {
59	6109487x	return width() * height();
60		}
61
62	1604x	double perimeter() const {
63	1604x	return 2 * width() + 2 * height();
64		}
65
66	383901x	bool intersects(const Box & other) const {
67	383901x	if (other.ymin > ymax)
68	48x	return false;
69	383853x	if (other.ymax < ymin)
70	44x	return false;
71	383809x	if (other.xmin > xmax)
72	!	return false;
73	383809x	if (other.xmax < xmin)
74	!	return false;
75
76	383809x	return true;
77		}
78
79	3937853x	Box intersection(const Box & other) const {
80		return {
81	3937853x	std::max(xmin, other.xmin),
82	3937853x	std::max(ymin, other.ymin),
83	3937853x	std::min(xmax, other.xmax),
84	3937853x	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	500587x	bool empty() const {
118	500587x	return xmin >= xmax \|\| ymin >= ymax;
119		}
120
121	3x	Box expand_to_include(const Box & other) const {
122	3x	if (empty()) {
123	!	return other;
124		}
125
126	3x	if (other.empty()) {
127	!	return *this;
128		}
129
130	3x	return { std::min(xmin, other.xmin),
131	3x	std::min(ymin, other.ymin),
132	3x	std::max(xmax, other.xmax),
133	3x	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	833x	bool operator==(const Box& other) const {
143	833x	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) 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;
31		double t = 0;
32
33		public:
34		/** \brief Update variance estimate with another value
35		*
36		* @param x value to add
37		* @param w weight of `x`
38		*/
39	537255x	void process(double x, double w) {
40	537255x	double mean_old = mean;
41
42	537255x	sum_w += w;
43	537255x	mean += (w / sum_w) * (x - mean_old);
44	537255x	t += w * (x - mean_old) * (x - mean);
45		}
46
47		/** \brief Return the population variance.
48		*/
49	3x	constexpr double variance() const {
50	3x	return t / sum_w;
51		}
52
53		/** \brief Return the population standard deviation
54		*/
55	2x	double stdev() const {
56	2x	return std::sqrt(variance());
57		}
58
59		/** \brief Return the population coefficient of variation
60		*/
61	1x	double coefficent_of_variation() const {
62	1x	return stdev() / mean;
63		}
64
65		};
66
67		}
68
69		#endif //EXACTEXTRACT_VARIANCE_H

1		// Copyright (c) 2019-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_WEIGHTED_QUANTILES_H
15		#define EXACTEXTRACT_WEIGHTED_QUANTILES_H
16
17		#include <algorithm>
18		#include <cmath>
19		#include <stdexcept>
20		#include <vector>
21
22		namespace exactextract {
23
24		class WeightedQuantiles {
25		// Compute weighted quantiles based on https://stats.stackexchange.com/a/13223
26
27		public:
28
29	55x	void process(double x, double w) {
30	55x	if (w < 0) {
31	!	throw std::runtime_error("Weighted quantile calculation does not support negative weights.");
32		}
33
34	55x	if (!std::isfinite(w)) {
35	!	throw std::runtime_error("Weighted quantile does not support non-finite weights.");
36		}
37
38	55x	m_ready_to_query = false;
39
40	55x	m_elems.emplace_back(x, w);
41		}
42
43		double quantile(double q) const;
44
45		private:
46		struct elem_t {
47	65x	elem_t(double _x, double _w) : x(_x), w(_w), cumsum(0), s(0) {}
48
49		double x;
50		double w;
51		double cumsum;
52		double s;
53		};
54
55		void prepare() const;
56
57		mutable std::vector<elem_t> m_elems;
58		mutable double m_sum_w;
59		mutable bool m_ready_to_query;
60		};
61
62		}
63
64		#endif //EXACTEXTRACT_WEIGHTED_QUANTILES_H