exactextractr coverage - 87.53%

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('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-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 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{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     include_xy if \code{TRUE}, augment the returned data frame with
#'                        columns for cell center coordinates (\code{x} and
#'                        \code{y}) or pass them to \code{fun}
#' @param     include_cell if \code{TRUE}, augment the returned data frame with
#'                        column for cell index
#' @param     fun an optional function or character vector, as described below
#' @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     weights  a weighting raster to be used with the \code{weighted_mean}
#'                     and \code{weighted_sum} summary operations.
#' @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) {
  exact_extract(x, sf::st_geometry(y), fun=fun, ..., include_xy=include_xy, progress=progress, max_cells_in_memory=max_cells_in_memory, include_cell=include_cell)
})

# 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) {
  if(inherits(y, 'sfc_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) && length(list(...)) > 0) {
    stop("exact_extract was called with a named summary operation that ",
         "does not accept additional arguments ...")
  }

  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 <- readStart(x)
    if (!is.null(weights)) {
      weights <- readStart(weights)
    }

    if (is.character(fun)) {
      results <- sapply(sf::st_as_binary(y, EWKB=TRUE), function(wkb) {
        update_progress()
        CPP_stats(x, weights, wkb, fun, max_cells_in_memory)
      })

      if (length(fun) == 1 && raster::nlayers(x) == 1) {
        # Just return a vector of stat results
        return(as.vector(results))
      } else {
        # Return a data frame with a column for each stat
        if (raster::nlayers(x) > 1) {
          z <- expand.grid(names(x), fun, stringsAsFactors=TRUE)
          colnames <- mapply(paste, z[[2]], z[[1]], MoreArgs=list(sep='.'))
        } else {
          colnames <- fun
        }

        results <- t(results)
        dimnames(results) <- list(NULL, colnames)
        return(as.data.frame(results))
      }
    } else {
      if (is.null(fun)) {
        appfn <- lapply # return list of data frames
      } else {
        appfn <- sapply
      }

      appfn(sf::st_as_binary(y, EWKB=TRUE), function(wkb) {
        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) {
          if (nrow(vals) == 0) {
            vals$x <- numeric()
            vals$y <- numeric()
          } else {
            x_coords <- raster::xFromCol(x, col=ret$col:(ret$col+ncol(ret$weights) - 1))
            y_coords <- raster::yFromRow(x, row=ret$row:(ret$row+nrow(ret$weights) - 1))

            vals$x <- rep.int(x_coords, times=nrow(ret$weights))
            vals$y <- rep(y_coords, each=ncol(ret$weights))
          }
        }

        if (include_cell) {
          if (nrow(vals) == 0) {
            vals$cell <- numeric()
          } else {
            rows <- rep(ret$row:(ret$row+nrow(ret$weights) - 1), each=ncol(ret$weights))
            cols <- rep.int(ret$col:(ret$col+ncol(ret$weights) - 1),
                                     times = nrow(ret$weights))
            vals$cell <- raster::cellFromRowCol(x,
                                           row=rows,
                                           col=cols)
          }
        }
        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) {
            return(fun(vals[,1], cov_fracs, ...))
          } else {
            return(fun(vals, cov_fracs, ...))
          }
        }
      })
    }
  }, finally={
    readStop(x)
    if (!is.null(weights)) {
      readStop(weights)
    }
  })
}

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

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

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

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('coverage_fraction')) {
15	12x	setGeneric('coverage_fraction', function(x, y, crop=FALSE, ...)
16	12x	standardGeneric('coverage_fraction'))
17		}
18
19		.coverage_fraction <- function(x, y, crop) {
20	12x	if(is.na(sf::st_crs(x)) && !is.na(sf::st_crs(y))) {
21	1x	warning("No CRS specified for raster; assuming it has the same CRS as the polygons.")
22	11x	} else if(is.na(sf::st_crs(y)) && !is.na(sf::st_crs(x))) {
23	1x	warning("No CRS specified for polygons; assuming they have the same CRS as the raster.")
24	10x	} else if(sf::st_crs(x) != sf::st_crs(y)) {
25	1x	y <- sf::st_transform(y, sf::st_crs(x))
26	1x	warning("Polygons transformed to raster CRS (EPSG:", sf::st_crs(x)$epsg, ")")
27		}
28
29	12x	lapply(sf::st_as_binary(y, EWKB=TRUE), function(wkb) {
30	12x	CPP_coverage_fraction(x, wkb, crop)
31		})
32		}
33
34		#' Compute the fraction of raster cells covered by a polygon
35		#'
36		#' @param x a (possibly empty) \code{RasterLayer} whose resolution and
37		#' extent will be used for the generated \code{RasterLayer}.
38		#' @param y a \code{sf} object with polygonal geometries
39		#' @param crop if \code{TRUE}, each generated \code{RasterLayer} will be
40		#' cropped to the extent of its associated feature.
41		#' @return a list with a \code{RasterLayer} for each feature in \code{y}.
42		#' Values of the raster represent the fraction of each
43		#' cell in \code{x} that is covered by \code{y}.
44		#' @examples
45		#' rast <- raster::raster(matrix(1:100, ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
46		#' poly <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))')
47		#'
48		#' cov_frac <- coverage_fraction(rast, poly)[[1]]
49		#' @name coverage_fraction
50		NULL
51
52		#' @import sf
53		#' @import raster
54		#' @useDynLib exactextractr
55		#' @rdname coverage_fraction
56		#' @export
57		setMethod('coverage_fraction', signature(x='RasterLayer', y='sf'), function(x, y, crop=FALSE) {
58	!	coverage_fraction(x, sf::st_geometry(y), crop)
59		})
60
61		#' @rdname coverage_fraction
62		#' @export
63		setMethod('coverage_fraction', signature(x='RasterLayer', y='sfc_MULTIPOLYGON'), .coverage_fraction)
64
65		#' @rdname coverage_fraction
66		#' @export
67		setMethod('coverage_fraction', signature(x='RasterLayer', y='sfc_POLYGON'), .coverage_fraction)
68

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::Matrix;
31		using exactextract::Grid;
32		using exactextract::subdivide;
33		using exactextract::bounded_extent;
34		using exactextract::Raster;
35		using exactextract::RasterView;
36		using exactextract::raster_cell_intersection;
37		using exactextract::RasterStats;
38		using exactextract::RasterSource;
39
40	130x	static Grid<bounded_extent> make_grid(const Rcpp::S4 & rast) {
41	390x	Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster");
42
43	260x	Rcpp::S4 extent = rast.slot("extent");
44
45	390x	Rcpp::Function resFn = raster["res"];
46
47	130x	Rcpp::NumericVector res = resFn(rast);
48
49		return {{
50	260x	extent.slot("xmin"),
51	260x	extent.slot("ymin"),
52	260x	extent.slot("xmax"),
53	260x	extent.slot("ymax"),
54		},
55	130x	res[0],
56	260x	res[1]
57		};
58		}
59
60		// Construct a Raster using an R vector for storage
61		class NumericVectorRaster : public exactextract::AbstractRaster<double> {
62		public:
63	121x	NumericVectorRaster(const Rcpp::NumericVector & vec, const Grid<bounded_extent> & g) :
64		AbstractRaster<double>(g),
65	121x	m_vec(vec)
66		{}
67
68	3529x	double operator()(size_t row, size_t col) const final {
69	3529x	return m_vec[row*cols() + col];
70		}
71
72		private:
73		const Rcpp::NumericVector m_vec;
74		};
75
76		// Read raster values from an R raster object
77		class S4RasterSource {
78		public:
79	95x	S4RasterSource(Rcpp::S4 rast) :
80	95x	m_grid(Grid<bounded_extent>::make_empty()),
81		m_rast(rast),
82	95x	m_last_box(0, 0, 0, 0)
83		{
84	95x	m_grid = make_grid(rast);
85		}
86
87	92x	const Grid<bounded_extent> &grid() const {
88	92x	return m_grid;
89		}
90
91	121x	std::unique_ptr<exactextract::AbstractRaster<double>> read_box(const exactextract::Box & box, int layer) {
92	121x	auto cropped_grid = m_grid.crop(box);
93
94	121x	if (!(box == m_last_box)) {
95	116x	m_last_box = box;
96
97	348x	Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster");
98	348x	Rcpp::Function getValuesBlockFn = raster["getValuesBlock"];
99
100	116x	if (cropped_grid.empty()) {
101	1x	m_rast_values = Rcpp::no_init(0, 0);
102		} else {
103		// Instead of reading only values for the requested band, we read values
104		// for all requested bands and then cache them to return from subsequent
105		// calls to read_box. There are two reasons for this:
106		// 1) Use of the 'lyrs` argument to getValuesBlock does not work for
107		// a single band in the CRAN version of raster
108		// 2) Reading everything once avoids the very significant per-call
109		// overhead of getValuesBlock, which largely comes from operations
110		// like setting names on the result vector.
111	115x	m_rast_values = getValuesBlockFn(m_rast,
112	!	1 + cropped_grid.row_offset(m_grid),
113	115x	cropped_grid.rows(),
114	115x	1 + cropped_grid.col_offset(m_grid),
115	115x	cropped_grid.cols(),
116	230x	Rcpp::Named("format", "m"));
117		}
118		}
119
120	121x	if (cropped_grid.empty()) {
121	1x	return std::make_unique<NumericVectorRaster>(m_rast_values, cropped_grid);
122		} else {
123	240x	return std::make_unique<NumericVectorRaster>(m_rast_values(Rcpp::_, layer),
124	120x	cropped_grid);
125		}
126		}
127
128		private:
129		Grid<bounded_extent> m_grid;
130		Rcpp::S4 m_rast;
131		Rcpp::NumericMatrix m_rast_values;
132		exactextract::Box m_last_box;
133		};
134
135		// GEOS warning handler
136	!	static void geos_warn(const char* fmt, ...) {
137	!	char buf[BUFSIZ] = { '\0' };
138
139		va_list msg;
140	!	va_start(msg, fmt);
141	!	vsnprintf(buf, BUFSIZ*sizeof(char), fmt, msg);
142	!	va_end(msg);
143
144	!	Rcpp::Function warning("warning");
145	!	warning(buf);
146		}
147
148		// GEOS error handler
149	!	static void geos_error(const char* fmt, ...) {
150	!	char buf[BUFSIZ] = { '\0' };
151
152		va_list msg;
153	!	va_start(msg, fmt);
154	!	vsnprintf(buf, BUFSIZ*sizeof(char), fmt, msg);
155	!	va_end(msg);
156
157	!	Rcpp::stop(buf);
158		}
159
160		// GEOSContextHandle wrapper to ensure finishGEOS is called.
161		struct GEOSAutoHandle {
162	107x	GEOSAutoHandle() {
163	107x	handle = initGEOS_r(geos_warn, geos_error);
164		}
165
166	214x	~GEOSAutoHandle() {
167	107x	finishGEOS_r(handle);
168		}
169
170		GEOSContextHandle_t handle;
171		};
172
173		// Return a smart pointer to a Geometry, given WKB input
174	103x	static geom_ptr read_wkb(const GEOSContextHandle_t & context, const Rcpp::RawVector & wkb) {
175	309x	wkb_reader_ptr wkb_reader{ GEOSWKBReader_create_r(context), [context](GEOSWKBReader* r) { GEOSWKBReader_destroy_r(context, r); } };
176
177		geom_ptr geom{ GEOSWKBReader_read_r(context,
178		wkb_reader.get(),
179	206x	std::addressof(wkb[0]),
180	103x	wkb.size()),
181	206x	[context](GEOSGeometry* g) { GEOSGeom_destroy_r(context, g); } };
182
183	103x	if (geom.get() == nullptr) {
184	!	Rcpp::stop("Failed to parse WKB geometry");
185		}
186
187	206x	return geom;
188		}
189
190		// [[Rcpp::export]]
191	23x	Rcpp::List CPP_exact_extract(Rcpp::S4 & rast, const Rcpp::RawVector & wkb) {
192	46x	GEOSAutoHandle geos;
193
194	23x	auto grid = make_grid(rast);
195	46x	auto coverage_fractions = raster_cell_intersection(grid, geos.handle, read_wkb(geos.handle, wkb).get());
196
197	23x	size_t nrow = coverage_fractions.rows();
198	23x	size_t ncol = coverage_fractions.cols();
199
200	23x	Rcpp::NumericMatrix weights = Rcpp::no_init(nrow, ncol);
201	69x	for (size_t i = 0; i < nrow; i++) {
202	199x	for (size_t j = 0; j < ncol; j++) {
203	153x	weights(i, j) = coverage_fractions(i ,j);
204		}
205		}
206
207	23x	if (nrow > 0) {
208	17x	size_t row_us = 1 + coverage_fractions.grid().row_offset(grid);
209	17x	if (row_us > static_cast<size_t>(std::numeric_limits<int>::max())) {
210
211		}
212		}
213
214	23x	int row = NA_INTEGER;
215	23x	int col = NA_INTEGER;
216	23x	if (nrow > 0) {
217	17x	size_t row_us = (1 + coverage_fractions.grid().row_offset(grid));
218	17x	if (row_us > static_cast<size_t>(std::numeric_limits<int>::max())) {
219	!	throw std::runtime_error("Cannot represent row offset as an R integer");
220		}
221	17x	row = static_cast<int>(row_us);
222		}
223	23x	if (ncol > 0) {
224	17x	size_t col_us = (1 + coverage_fractions.grid().col_offset(grid));
225	17x	if (col_us > static_cast<size_t>(std::numeric_limits<int>::max())) {
226	!	throw std::runtime_error("Cannot represent column offset as an R integer");
227		}
228	17x	col = static_cast<int>(col_us);
229		}
230
231		return Rcpp::List::create(
232	46x	Rcpp::Named("row") = row,
233	46x	Rcpp::Named("col") = col,
234	46x	Rcpp::Named("weights") = weights
235		);
236		}
237
238		// [[Rcpp::export]]
239	12x	Rcpp::S4 CPP_coverage_fraction(Rcpp::S4 & rast, const Rcpp::RawVector & wkb, bool crop)
240		{
241	24x	GEOSAutoHandle geos;
242	36x	Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster");
243	36x	Rcpp::Function rasterFn = raster["raster"];
244	36x	Rcpp::Function crsFn = raster["crs"];
245
246	12x	auto grid = make_grid(rast);
247	24x	auto coverage_fraction = raster_cell_intersection(grid, geos.handle, read_wkb(geos.handle, wkb).get());
248
249	12x	if (crop) {
250	2x	grid = coverage_fraction.grid();
251		}
252	24x	RasterView<float> coverage_view(coverage_fraction, grid);
253
254	12x	Rcpp::NumericMatrix weights{static_cast<int>(grid.rows()),
255	12x	static_cast<int>(grid.cols())};
256
257	629x	for (size_t i = 0; i < grid.rows(); i++) {
258	159018x	for (size_t j = 0; j < grid.cols(); j++) {
259	158401x	weights(i, j) = coverage_view(i, j);
260	158401x	if (!crop && std::isnan(weights(i, j))) {
261	139692x	weights(i, j) = 0;
262		}
263		}
264		}
265
266		return rasterFn(weights,
267	24x	Rcpp::Named("xmn")=grid.xmin(),
268	24x	Rcpp::Named("xmx")=grid.xmax(),
269	24x	Rcpp::Named("ymn")=grid.ymin(),
270	24x	Rcpp::Named("ymx")=grid.ymax(),
271	48x	Rcpp::Named("crs")=crsFn(rast));
272		}
273
274	96x	static int get_nlayers(Rcpp::S4 & rast) {
275	288x	Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster");
276	288x	Rcpp::Function nlayersFn = raster["nlayers"];
277
278	96x	Rcpp::NumericVector nlayersVec = nlayersFn(rast);
279
280	192x	return static_cast<int>(nlayersVec[0]);
281		}
282
283		// [[Rcpp::export]]
284	72x	Rcpp::NumericMatrix CPP_stats(Rcpp::S4 & rast,
285		Rcpp::Nullable<Rcpp::S4> weights,
286		const Rcpp::RawVector & wkb,
287		const Rcpp::StringVector & stats,
288		int max_cells_in_memory) {
289		try {
290	144x	GEOSAutoHandle geos;
291
292	72x	if (max_cells_in_memory < 1) {
293	1x	Rcpp::stop("Invalid value for max_cells_in_memory: %d", max_cells_in_memory);
294		}
295
296	71x	int nlayers = get_nlayers(rast);
297
298	142x	S4RasterSource rsrc(rast);
299
300	71x	std::unique_ptr<S4RasterSource> rweights;
301	71x	bool weighted = false;
302	71x	if (weights.isNotNull()) {
303	26x	Rcpp::S4 weights_s4 = weights.get();
304
305	25x	if (get_nlayers(weights_s4) != 1) {
306	1x	Rcpp::stop("Weighting raster must have only a single layer.");
307		}
308
309	24x	rweights = std::make_unique<S4RasterSource>(weights_s4);
310	24x	weighted = true;
311		}
312
313	70x	bool store_values = false;
314
315	147x	for (const auto & stat : stats) {
316		// explicit construction of std::string seems necessary to avoid ambiguous overload error
317	158x	if (stat == std::string("mode") \|\|
318	153x	stat == std::string("majority") \|\|
319	378x	stat == std::string("minority") \|\|
320	146x	stat == std::string("variety")) {
321	20x	store_values = true;
322		}
323
324	133x	if (!weighted &&
325	133x	(stat == std::string("weighted_mean") \|\|
326	132x	stat == std::string("weighted_sum"))) {
327	2x	Rcpp::stop("Weighted stat requested but no weights provided.");
328		}
329		}
330
331	136x	std::vector<RasterStats<double>> raster_stats;
332	68x	raster_stats.reserve(nlayers);
333	141x	for (int i = 0; i < nlayers; i++) {
334	73x	raster_stats.emplace_back(store_values);
335		}
336
337	136x	Rcpp::NumericMatrix stat_results = Rcpp::no_init(nlayers, stats.size());
338
339	136x	auto geom = read_wkb(geos.handle, wkb);
340	68x	auto bbox = exactextract::geos_get_box(geos.handle, geom.get());
341
342	68x	auto grid = weighted ? rsrc.grid().common_grid(rweights->grid()) : rsrc.grid();
343
344	67x	if (bbox.intersects(grid.extent())) {
345	58x	auto cropped_grid = grid.crop(bbox);
346
347	151x	for (const auto &subgrid : subdivide(cropped_grid, max_cells_in_memory)) {
348	186x	auto coverage_fraction = raster_cell_intersection(subgrid, geos.handle, geom.get());
349	93x	auto& cov_grid = coverage_fraction.grid();
350
351	93x	if (!cov_grid.empty()) {
352	93x	if (weighted) {
353	46x	auto weights = rweights->read_box(cov_grid.extent(), 0);
354
355	47x	for (int i = 0; i < nlayers; i++) {
356	48x	auto values = rsrc.read_box(cov_grid.extent(), i);
357	24x	raster_stats[i].process(coverage_fraction, values, weights);
358		}
359		} else {
360	144x	for (int i = 0; i < nlayers; i++) {
361	148x	auto values = rsrc.read_box(cov_grid.extent(), i);
362	74x	raster_stats[i].process(coverage_fraction, *values);
363		}
364		}
365		}
366		}
367		}
368
369	138x	for (int j = 0; j < nlayers; j++) {
370	154x	for(int i = 0; i < stats.size(); i++) {
371	83x	if (stats[i] == std::string("mean")) stat_results(j, i) = raster_stats[j].mean();
372
373	72x	else if (stats[i] == std::string("sum")) stat_results(j, i) = raster_stats[j].sum();
374	61x	else if (stats[i] == std::string("count")) stat_results(j, i) = raster_stats[j].count();
375
376	55x	else if (stats[i] == std::string("min")) stat_results(j, i) = raster_stats[j].min().value_or(NA_REAL);
377	49x	else if (stats[i] == std::string("max")) stat_results(j, i) = raster_stats[j].max().value_or(NA_REAL);
378
379	43x	else if (stats[i] == std::string("mode")) stat_results(j, i) = raster_stats[j].mode().value_or(NA_REAL);
380	38x	else if (stats[i] == std::string("majority")) stat_results(j, i) = raster_stats[j].mode().value_or(NA_REAL);
381	35x	else if (stats[i] == std::string("minority")) stat_results(j, i) = raster_stats[j].minority().value_or(NA_REAL);
382
383	31x	else if (stats[i] == std::string("variety")) stat_results(j, i) = raster_stats[j].variety();
384	19x	else if (stats[i] == std::string("weighted_mean")) stat_results(j, i) = raster_stats[j].weighted_mean();
385	7x	else if (stats[i] == std::string("weighted_sum")) stat_results(j, i) = raster_stats[j].weighted_sum();
386
387	4x	else if (stats[i] == std::string("variance")) stat_results(j, i) = raster_stats[j].variance();
388	3x	else if (stats[i] == std::string("stdev")) stat_results(j, i) = raster_stats[j].stdev();
389	2x	else if (stats[i] == std::string("coefficient_of_variation")) stat_results(j, i) = raster_stats[j].coefficient_of_variation();
390
391	4x	else Rcpp::stop("Unknown stat: " + stats[i]);
392		}
393		}
394
395	132x	return stat_results;
396	12x	} catch (std::exception & e) {
397		// throw predictible exception class
398	6x	Rcpp::stop(e.what());
399		}
400		}

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

1		// Copyright (c) 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	1961x	void process(double x, double w) {
40	1961x	double mean_old = mean;
41
42	1961x	sum_w += w;
43	1961x	mean += (w / sum_w) * (x - mean_old);
44	1961x	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