/src/gdal/alg/raster_stats.h

Source
/******************************************************************************
*
 * Project:  GDAL
 * Purpose:  GDALZonalStats implementation
 * Author:   Dan Baston
 *
 ******************************************************************************
 * Copyright (c) 2018-2025, ISciences LLC
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#pragma once

#include <algorithm>
#include <cmath>
#include <limits>
#include <optional>
#include <unordered_map>

//! @cond Doxygen_Suppress

namespace gdal
{
struct RasterStatsOptions
{
    static constexpr float min_coverage_fraction_default =
        std::numeric_limits<float>::min();  // ~1e-38

    float min_coverage_fraction = min_coverage_fraction_default;
    bool calc_median = false;
    bool calc_variance = false;
    bool store_histogram = false;
    bool store_values = false;
    bool store_weights = false;
    bool store_coverage_fraction = false;
    bool store_xy = false;
    bool include_nodata = false;
    double default_weight = std::numeric_limits<double>::quiet_NaN();

    bool operator==(const RasterStatsOptions &other) const
    {
        return min_coverage_fraction == other.min_coverage_fraction &&
               calc_median == other.calc_median &&
               calc_variance == other.calc_variance &&
               store_histogram == other.store_histogram &&
               store_values == other.store_values &&
               store_weights == other.store_weights &&
               store_coverage_fraction == other.store_coverage_fraction &&
               store_xy == other.store_xy &&
               include_nodata == other.include_nodata &&
               (default_weight == other.default_weight ||
                (std::isnan(default_weight) &&
                 std::isnan(other.default_weight)));
    }

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

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)
    {
        if (w == 0)
        {
            return;
        }

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

template <typename ValueType> 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(const RasterStatsOptions &options)
        : m_min{std::numeric_limits<ValueType>::max()},
          m_max{std::numeric_limits<ValueType>::lowest()}, m_sum_ciwi{0},
          m_sum_ci{0}, m_sum_xici{0}, m_sum_xiciwi{0}, m_options{options}
    {
    }

    // All pixels covered 100%
    void process(const ValueType *pValues, const GByte *pabyMask,
                 const double *padfWeights, const GByte *pabyWeightsMask,
                 const double *padfX, const double *padfY, size_t nX, size_t nY)
    {
        for (size_t i = 0; i < nX * nY; i++)
        {
            if (pabyMask[i] == 255)
            {
                if (padfX && padfY)
                {
                    process_location(padfX[i % nX], padfY[i / nX]);
                }
                const double dfWeight =
                    pabyWeightsMask
                        ? (pabyWeightsMask[i] == 255
                               ? padfWeights[i]
                               : std::numeric_limits<double>::quiet_NaN())
                        : 1.0;
                process_value(pValues[i], 1.0, dfWeight);
            }
        }
    }

    // Pixels covered 0% or 100%
    void process(const ValueType *pValues, const GByte *pabyMask,
                 const double *padfWeights, const GByte *pabyWeightsMask,
                 const GByte *pabyCov, const double *pdfX, const double *pdfY,
                 size_t nX, size_t nY)
    {
        for (size_t i = 0; i < nX * nY; i++)
        {
            if (pabyMask[i] == 255 && pabyCov[i])
            {
                if (pdfX && pdfY)
                {
                    process_location(pdfX[i % nX], pdfY[i / nX]);
                }
                const double dfWeight =
                    pabyWeightsMask
                        ? (pabyWeightsMask[i] == 255
                               ? padfWeights[i]
                               : std::numeric_limits<double>::quiet_NaN())
                        : 1.0;
                process_value(pValues[i], 1.0, dfWeight);
            }
        }
    }

    // Pixels fractionally covered
    void process(const ValueType *pValues, const GByte *pabyMask,
                 const double *padfWeights, const GByte *pabyWeightsMask,
                 const float *pfCov, const double *pdfX, const double *pdfY,
                 size_t nX, size_t nY)
    {
        for (size_t i = 0; i < nX * nY; i++)
        {
            if (pabyMask[i] == 255 &&
                pfCov[i] >= m_options.min_coverage_fraction)
            {
                if (pdfX && pdfY)
                {
                    process_location(pdfX[i % nX], pdfY[i / nX]);
                }
                const double dfWeight =
                    pabyWeightsMask
                        ? (pabyWeightsMask[i] == 255
                               ? padfWeights[i]
                               : std::numeric_limits<double>::quiet_NaN())
                        : 1.0;
                process_value(pValues[i], pfCov[i], dfWeight);
            }
        }
    }

    void process_location(double x, double y)
    {
        if (m_options.store_xy)
        {
            m_cell_x.push_back(x);
            m_cell_y.push_back(y);
        }
    }

    void process_value(const ValueType &val, float coverage, double weight)
    {
        if (m_options.store_coverage_fraction)
        {
            m_cell_cov.push_back(coverage);
        }

        m_sum_ci += static_cast<double>(coverage);
        m_sum_xici += static_cast<double>(val) * static_cast<double>(coverage);

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

        if (m_options.calc_variance)
        {
            m_variance.process(static_cast<double>(val),
                               static_cast<double>(coverage));
            m_weighted_variance.process(static_cast<double>(val), ciwi);
        }

        if (val < m_min)
        {
            m_min = val;
            if (m_options.store_xy)
            {
                m_min_xy = {m_cell_x.back(), m_cell_y.back()};
            }
        }

        if (val > m_max)
        {
            m_max = val;
            if (m_options.store_xy)
            {
                m_max_xy = {m_cell_x.back(), m_cell_y.back()};
            }
        }

        if (m_options.store_histogram)
        {
            auto &entry = m_freq[val];
            entry.m_sum_ci += static_cast<double>(coverage);
            entry.m_sum_ciwi += ciwi;
        }

        if (m_options.store_values || m_options.calc_median)
        {
            m_cell_values.push_back(val);
        }

        if (m_options.store_weights)
        {
            m_cell_weights.push_back(weight);
        }
    }

    /**
     * The mean value of cells covered by this polygon, weighted
     * by the percent of the cell that is covered.
     */
    double mean() const
    {
        if (count() > 0)
        {
            return sum() / count();
        }
        else
        {
            return std::numeric_limits<double>::quiet_NaN();
        }
    }

    /** The median value of cells touched by this polygon. The cell
     *  coverage fraction is not taken into account. */
    double median() const
    {
        auto n = m_cell_values.size();
        if (n == 0)
        {
            return std::numeric_limits<double>::quiet_NaN();
        }
        if (n == 1)
        {
            return static_cast<double>(m_cell_values[0]);
        }
        if (n == 2)
        {
            return 0.5 * (static_cast<double>(m_cell_values[0]) +
                          static_cast<double>(m_cell_values[1]));
        }

        std::vector<ValueType> *values = nullptr;
        std::unique_ptr<std::vector<ValueType>> values_copy;

        if (m_options.store_values)
        {
            values_copy =
                std::make_unique<std::vector<ValueType>>(m_cell_values);
            values = values_copy.get();
        }
        else
        {
            values = const_cast<std::vector<ValueType> *>(&m_cell_values);
        }

        auto mid = std::next(values->begin(), n / 2);
        std::nth_element(values->begin(), mid, values->end());

        if (n % 2 != 0)
        {
            return static_cast<double>(*mid);
        }

        auto lhs_max = std::max_element(values->begin(), mid);

        return 0.5 *
               (static_cast<double>(*lhs_max) + static_cast<double>(*mid));
    }

    /**
     * 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.
     */
    double weighted_mean() const
    {
        if (weighted_count() > 0)
        {
            return weighted_sum() / weighted_count();
        }
        else
        {
            return std::numeric_limits<double>::quiet_NaN();
        }
    }

    /** The fraction of weighted cells to unweighted cells.
     *  Meaningful only when the values of the weighting
     *  raster are between 0 and 1.
     */
    double 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.
     */
    std::optional<ValueType> mode() const
    {
        auto it = std::max_element(
            m_freq.cbegin(), m_freq.cend(),
            [](const auto &a, const auto &b)
            {
                return a.second.m_sum_ci < b.second.m_sum_ci ||
                       (a.second.m_sum_ci == b.second.m_sum_ci &&
                        a.first < b.first);
            });
        if (it == m_freq.end())
        {
            return std::nullopt;
        }
        return it->first;
    }

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

    /// XY values corresponding to the center of the cell whose value
    /// is returned by min()
    std::optional<std::pair<double, double>> min_xy() const
    {
        if (m_sum_ci == 0)
        {
            return std::nullopt;
        }
        return m_min_xy;
    }

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

    /// XY values corresponding to the center of the cell whose value
    /// is returned by max()
    std::optional<std::pair<double, double>> max_xy() const
    {
        if (m_sum_ci == 0)
        {
            return std::nullopt;
        }
        return m_max_xy;
    }

    /**
     * The sum of raster cells covered by the polygon, with each raster
     * value weighted by its coverage fraction.
     */
    double sum() const
    {
        return 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.
     */
    double weighted_sum() const
    {
        return m_sum_xiciwi;
    }

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

    /**
     * The number of raster cells with a specific value
     * covered by the polygon. Weights are not taken
     * into account.
     */
    std::optional<double> count(const ValueType &value) const
    {
        const auto &entry = m_freq.find(value);

        if (entry == m_freq.end())
        {
            return std::nullopt;
        }

        return entry->second.m_sum_ci;
    }

    /**
     * The fraction of defined raster cells covered by the polygon with
     * a value that equals the specified value.
     * Weights are not taken into account.
     */
    std::optional<double> frac(const ValueType &value) const
    {
        auto count_for_value = count(value);

        if (!count_for_value.has_value())
        {
            return count_for_value;
        }

        return count_for_value.value() / count();
    }

    /**
     * The weighted fraction of defined raster cells covered by the polygon with
     * a value that equals the specified value.
     */
    std::optional<double> weighted_frac(const ValueType &value) const
    {
        auto count_for_value = weighted_count(value);

        if (!count_for_value.has_value())
        {
            return count_for_value;
        }

        return count_for_value.value() / weighted_count();
    }

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

    /**
     * The population variance of raster cells touched
     * by the polygon, taking into account cell coverage
     * fractions and values of a weighting raster.
     */
    double weighted_variance() const
    {
        return m_weighted_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.
     */
    double stdev() const
    {
        return m_variance.stdev();
    }

    /**
     * The population standard deviation of raster cells
     * touched by the polygon, taking into account cell
     * coverage fractions and values of a weighting raster.
     */
    double weighted_stdev() const
    {
        return m_weighted_variance.stdev();
    }

    /**
     * The sum of weights for each cell covered by the
     * polygon, with each weight multiplied by the 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.
     */
    double weighted_count() const
    {
        return m_sum_ciwi;
    }

    /**
     * The sum of weights for each cell of a specific value covered by the
     * polygon, with each weight multiplied by the 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.
     */
    std::optional<double> weighted_count(const ValueType &value) const
    {
        const auto &entry = m_freq.find(value);

        if (entry == m_freq.end())
        {
            return std::nullopt;
        }

        return entry->second.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.
     */
    std::optional<ValueType> minority() const
    {
        auto it = std::min_element(
            m_freq.cbegin(), m_freq.cend(),
            [](const auto &a, const auto &b)
            {
                return a.second.m_sum_ci < b.second.m_sum_ci ||
                       (a.second.m_sum_ci == b.second.m_sum_ci &&
                        a.first < b.first);
            });
        if (it == m_freq.end())
        {
            return std::nullopt;
        }
        return it->first;
    }

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

    const std::vector<ValueType> &values() const
    {
        return m_cell_values;
    }

    const std::vector<bool> &values_defined() const
    {
        return m_cell_values_defined;
    }

    const std::vector<float> &coverage_fractions() const
    {
        return m_cell_cov;
    }

    const std::vector<double> &weights() const
    {
        return m_cell_weights;
    }

    const std::vector<bool> &weights_defined() const
    {
        return m_cell_weights_defined;
    }

    const std::vector<double> &center_x() const
    {
        return m_cell_x;
    }

    const std::vector<double> &center_y() const
    {
        return m_cell_y;
    }

    const auto &freq() const
    {
        return m_freq;
    }

  private:
    ValueType m_min{};
    ValueType m_max{};
    std::pair<double, double> m_min_xy{
        std::numeric_limits<double>::quiet_NaN(),
        std::numeric_limits<double>::quiet_NaN()};
    std::pair<double, double> m_max_xy{
        std::numeric_limits<double>::quiet_NaN(),
        std::numeric_limits<double>::quiet_NaN()};

    // ci: coverage fraction of pixel i
    // wi: weight of pixel i
    // xi: value of pixel i
    double m_sum_ciwi{0};
    double m_sum_ci{0};
    double m_sum_xici{0};
    double m_sum_xiciwi{0};

    WestVariance m_variance{};
    WestVariance m_weighted_variance{};

    struct ValueFreqEntry
    {
        double m_sum_ci = 0;
        double m_sum_ciwi = 0;
    };

    std::unordered_map<ValueType, ValueFreqEntry> m_freq{};

    std::vector<float> m_cell_cov{};
    std::vector<ValueType> m_cell_values{};
    std::vector<double> m_cell_weights{};
    std::vector<double> m_cell_x{};
    std::vector<double> m_cell_y{};
    std::vector<bool> m_cell_values_defined{};
    std::vector<bool> m_cell_weights_defined{};

    RasterStatsOptions m_options;
};

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

}  // namespace gdal

//! @endcond

Coverage Report

Created: 2026-10-28 07:09

Line	Count	Source
1		/******************************************************************************
2		*
3		* Project: GDAL
4		* Purpose: GDALZonalStats implementation
5		* Author: Dan Baston
6		*
7		******************************************************************************
8		* Copyright (c) 2018-2025, ISciences LLC
9		*
10		* SPDX-License-Identifier: MIT
11		****************************************************************************/
12
13		#pragma once
14
15		#include <algorithm>
16		#include <cmath>
17		#include <limits>
18		#include <optional>
19		#include <unordered_map>
20
21		//! @cond Doxygen_Suppress
22
23		namespace gdal
24		{
25		struct RasterStatsOptions
26		{
27		static constexpr float min_coverage_fraction_default =
28		std::numeric_limits<float>::min(); // ~1e-38
29
30		float min_coverage_fraction = min_coverage_fraction_default;
31		bool calc_median = false;
32		bool calc_variance = false;
33		bool store_histogram = false;
34		bool store_values = false;
35		bool store_weights = false;
36		bool store_coverage_fraction = false;
37		bool store_xy = false;
38		bool include_nodata = false;
39		double default_weight = std::numeric_limits<double>::quiet_NaN();
40
41		bool operator==(const RasterStatsOptions &other) const
42	0	{
43	0	return min_coverage_fraction == other.min_coverage_fraction &&
44	0	calc_median == other.calc_median &&
45	0	calc_variance == other.calc_variance &&
46	0	store_histogram == other.store_histogram &&
47	0	store_values == other.store_values &&
48	0	store_weights == other.store_weights &&
49	0	store_coverage_fraction == other.store_coverage_fraction &&
50	0	store_xy == other.store_xy &&
51	0	include_nodata == other.include_nodata &&
52	0	(default_weight == other.default_weight \|\|
53	0	(std::isnan(default_weight) &&
54	0	std::isnan(other.default_weight)));
55	0	}
56
57		bool operator!=(const RasterStatsOptions &other) const
58	0	{
59	0	return !(*this == other);
60	0	}
61		};
62
63		class WestVariance
64		{
65		/** \brief Implements an incremental algorithm for weighted standard
66		* deviation, variance, and coefficient of variation, as described in
67		* formula WV2 of West, D.H.D. (1979) "Updating Mean and Variance
68		* Estimates: An Improved Method". Communications of the ACM 22(9).
69		*/
70
71		private:
72		double sum_w = 0;
73		double mean = 0;
74		double t = 0;
75
76		public:
77		/** \brief Update variance estimate with another value
78		*
79		* @param x value to add
80		* @param w weight of `x`
81		*/
82		void process(double x, double w)
83	0	{
84	0	if (w == 0)
85	0	{
86	0	return;
87	0	}
88
89	0	double mean_old = mean;
90
91	0	sum_w += w;
92	0	mean += (w / sum_w) * (x - mean_old);
93	0	t += w * (x - mean_old) * (x - mean);
94	0	}
95
96		/** \brief Return the population variance.
97		*/
98		constexpr double variance() const
99	0	{
100	0	return t / sum_w;
101	0	}
102
103		/** \brief Return the population standard deviation
104		*/
105		double stdev() const
106	0	{
107	0	return std::sqrt(variance());
108	0	}
109
110		/** \brief Return the population coefficient of variation
111		*/
112		double coefficent_of_variation() const
113	0	{
114	0	return stdev() / mean;
115	0	}
116		};
117
118		template <typename ValueType> class RasterStats
119		{
120		public:
121		/**
122		* Compute raster statistics from a Raster representing intersection percentages,
123		* a Raster representing data values, and (optionally) a Raster representing weights.
124		* and a set of raster values.
125		*/
126		explicit RasterStats(const RasterStatsOptions &options)
127	0	: m_min{std::numeric_limits<ValueType>::max()},
128	0	m_max{std::numeric_limits<ValueType>::lowest()}, m_sum_ciwi{0},
129	0	m_sum_ci{0}, m_sum_xici{0}, m_sum_xiciwi{0}, m_options{options}
130	0	{
131	0	}
132
133		// All pixels covered 100%
134		void process(const ValueType pValues, const GByte pabyMask,
135		const double padfWeights, const GByte pabyWeightsMask,
136		const double padfX, const double padfY, size_t nX, size_t nY)
137	0	{
138	0	for (size_t i = 0; i < nX * nY; i++)
139	0	{
140	0	if (pabyMask[i] == 255)
141	0	{
142	0	if (padfX && padfY)
143	0	{
144	0	process_location(padfX[i % nX], padfY[i / nX]);
145	0	}
146	0	const double dfWeight =
147	0	pabyWeightsMask
148	0	? (pabyWeightsMask[i] == 255
149	0	? padfWeights[i]
150	0	: std::numeric_limits<double>::quiet_NaN())
151	0	: 1.0;
152	0	process_value(pValues[i], 1.0, dfWeight);
153	0	}
154	0	}
155	0	}
156
157		// Pixels covered 0% or 100%
158		void process(const ValueType pValues, const GByte pabyMask,
159		const double padfWeights, const GByte pabyWeightsMask,
160		const GByte pabyCov, const double pdfX, const double *pdfY,
161		size_t nX, size_t nY)
162	0	{
163	0	for (size_t i = 0; i < nX * nY; i++)
164	0	{
165	0	if (pabyMask[i] == 255 && pabyCov[i])
166	0	{
167	0	if (pdfX && pdfY)
168	0	{
169	0	process_location(pdfX[i % nX], pdfY[i / nX]);
170	0	}
171	0	const double dfWeight =
172	0	pabyWeightsMask
173	0	? (pabyWeightsMask[i] == 255
174	0	? padfWeights[i]
175	0	: std::numeric_limits<double>::quiet_NaN())
176	0	: 1.0;
177	0	process_value(pValues[i], 1.0, dfWeight);
178	0	}
179	0	}
180	0	}
181
182		// Pixels fractionally covered
183		void process(const ValueType pValues, const GByte pabyMask,
184		const double padfWeights, const GByte pabyWeightsMask,
185		const float pfCov, const double pdfX, const double *pdfY,
186		size_t nX, size_t nY)
187	0	{
188	0	for (size_t i = 0; i < nX * nY; i++)
189	0	{
190	0	if (pabyMask[i] == 255 &&
191	0	pfCov[i] >= m_options.min_coverage_fraction)
192	0	{
193	0	if (pdfX && pdfY)
194	0	{
195	0	process_location(pdfX[i % nX], pdfY[i / nX]);
196	0	}
197	0	const double dfWeight =
198	0	pabyWeightsMask
199	0	? (pabyWeightsMask[i] == 255
200	0	? padfWeights[i]
201	0	: std::numeric_limits<double>::quiet_NaN())
202	0	: 1.0;
203	0	process_value(pValues[i], pfCov[i], dfWeight);
204	0	}
205	0	}
206	0	}
207
208		void process_location(double x, double y)
209	0	{
210	0	if (m_options.store_xy)
211	0	{
212	0	m_cell_x.push_back(x);
213	0	m_cell_y.push_back(y);
214	0	}
215	0	}
216
217		void process_value(const ValueType &val, float coverage, double weight)
218	0	{
219	0	if (m_options.store_coverage_fraction)
220	0	{
221	0	m_cell_cov.push_back(coverage);
222	0	}
223
224	0	m_sum_ci += static_cast<double>(coverage);
225	0	m_sum_xici += static_cast<double>(val) * static_cast<double>(coverage);
226
227	0	double ciwi = static_cast<double>(coverage) * weight;
228	0	m_sum_ciwi += ciwi;
229	0	m_sum_xiciwi += static_cast<double>(val) * ciwi;
230
231	0	if (m_options.calc_variance)
232	0	{
233	0	m_variance.process(static_cast<double>(val),
234	0	static_cast<double>(coverage));
235	0	m_weighted_variance.process(static_cast<double>(val), ciwi);
236	0	}
237
238	0	if (val < m_min)
239	0	{
240	0	m_min = val;
241	0	if (m_options.store_xy)
242	0	{
243	0	m_min_xy = {m_cell_x.back(), m_cell_y.back()};
244	0	}
245	0	}
246
247	0	if (val > m_max)
248	0	{
249	0	m_max = val;
250	0	if (m_options.store_xy)
251	0	{
252	0	m_max_xy = {m_cell_x.back(), m_cell_y.back()};
253	0	}
254	0	}
255
256	0	if (m_options.store_histogram)
257	0	{
258	0	auto &entry = m_freq[val];
259	0	entry.m_sum_ci += static_cast<double>(coverage);
260	0	entry.m_sum_ciwi += ciwi;
261	0	}
262
263	0	if (m_options.store_values \|\| m_options.calc_median)
264	0	{
265	0	m_cell_values.push_back(val);
266	0	}
267
268	0	if (m_options.store_weights)
269	0	{
270	0	m_cell_weights.push_back(weight);
271	0	}
272	0	}
273
274		/**
275		* The mean value of cells covered by this polygon, weighted
276		* by the percent of the cell that is covered.
277		*/
278		double mean() const
279	0	{
280	0	if (count() > 0)
281	0	{
282	0	return sum() / count();
283	0	}
284	0	else
285	0	{
286	0	return std::numeric_limits<double>::quiet_NaN();
287	0	}
288	0	}
289
290		/** The median value of cells touched by this polygon. The cell
291		* coverage fraction is not taken into account. */
292		double median() const
293	0	{
294	0	auto n = m_cell_values.size();
295	0	if (n == 0)
296	0	{
297	0	return std::numeric_limits<double>::quiet_NaN();
298	0	}
299	0	if (n == 1)
300	0	{
301	0	return static_cast<double>(m_cell_values[0]);
302	0	}
303	0	if (n == 2)
304	0	{
305	0	return 0.5 * (static_cast<double>(m_cell_values[0]) +
306	0	static_cast<double>(m_cell_values[1]));
307	0	}
308
309	0	std::vector<ValueType> *values = nullptr;
310	0	std::unique_ptr<std::vector<ValueType>> values_copy;
311
312	0	if (m_options.store_values)
313	0	{
314	0	values_copy =
315	0	std::make_unique<std::vector<ValueType>>(m_cell_values);
316	0	values = values_copy.get();
317	0	}
318	0	else
319	0	{
320	0	values = const_cast<std::vector<ValueType> *>(&m_cell_values);
321	0	}
322
323	0	auto mid = std::next(values->begin(), n / 2);
324	0	std::nth_element(values->begin(), mid, values->end());
325
326	0	if (n % 2 != 0)
327	0	{
328	0	return static_cast<double>(*mid);
329	0	}
330
331	0	auto lhs_max = std::max_element(values->begin(), mid);
332
333	0	return 0.5 *
334	0	(static_cast<double>(lhs_max) + static_cast<double>(mid));
335	0	}
336
337		/**
338		* The mean value of cells covered by this polygon, weighted
339		* by the percent of the cell that is covered and a secondary
340		* weighting raster.
341		*
342		* If any weights are undefined, will return NAN. If this is undesirable,
343		* caller should replace undefined weights with a suitable default
344		* before computing statistics.
345		*/
346		double weighted_mean() const
347	0	{
348	0	if (weighted_count() > 0)
349	0	{
350	0	return weighted_sum() / weighted_count();
351	0	}
352	0	else
353	0	{
354	0	return std::numeric_limits<double>::quiet_NaN();
355	0	}
356	0	}
357
358		/** The fraction of weighted cells to unweighted cells.
359		* Meaningful only when the values of the weighting
360		* raster are between 0 and 1.
361		*/
362		double weighted_fraction() const
363		{
364		return weighted_sum() / sum();
365		}
366
367		/**
368		* The raster value occupying the greatest number of cells
369		* or partial cells within the polygon. When multiple values
370		* cover the same number of cells, the greatest value will
371		* be returned. Weights are not taken into account.
372		*/
373		std::optional<ValueType> mode() const
374	0	{
375	0	auto it = std::max_element(
376	0	m_freq.cbegin(), m_freq.cend(),
377	0	[](const auto &a, const auto &b)
378	0	{
379	0	return a.second.m_sum_ci < b.second.m_sum_ci \|\|
380	0	(a.second.m_sum_ci == b.second.m_sum_ci &&
381	0	a.first < b.first);
382	0	});
383	0	if (it == m_freq.end())
384	0	{
385	0	return std::nullopt;
386	0	}
387	0	return it->first;
388	0	}
389
390		/**
391		* The minimum value in any raster cell wholly or partially covered
392		* by the polygon. Weights are not taken into account.
393		*/
394		std::optional<ValueType> min() const
395	0	{
396	0	if (m_sum_ci == 0)
397	0	{
398	0	return std::nullopt;
399	0	}
400	0	return m_min;
401	0	}
402
403		/// XY values corresponding to the center of the cell whose value
404		/// is returned by min()
405		std::optional<std::pair<double, double>> min_xy() const
406	0	{
407	0	if (m_sum_ci == 0)
408	0	{
409	0	return std::nullopt;
410	0	}
411	0	return m_min_xy;
412	0	}
413
414		/**
415		* The maximum value in any raster cell wholly or partially covered
416		* by the polygon. Weights are not taken into account.
417		*/
418		std::optional<ValueType> max() const
419	0	{
420	0	if (m_sum_ci == 0)
421	0	{
422	0	return std::nullopt;
423	0	}
424	0	return m_max;
425	0	}
426
427		/// XY values corresponding to the center of the cell whose value
428		/// is returned by max()
429		std::optional<std::pair<double, double>> max_xy() const
430	0	{
431	0	if (m_sum_ci == 0)
432	0	{
433	0	return std::nullopt;
434	0	}
435	0	return m_max_xy;
436	0	}
437
438		/**
439		* The sum of raster cells covered by the polygon, with each raster
440		* value weighted by its coverage fraction.
441		*/
442		double sum() const
443	0	{
444	0	return m_sum_xici;
445	0	}
446
447		/**
448		* The sum of raster cells covered by the polygon, with each raster
449		* value weighted by its coverage fraction and weighting raster value.
450		*
451		* If any weights are undefined, will return NAN. If this is undesirable,
452		* caller should replace undefined weights with a suitable default
453		* before computing statistics.
454		*/
455		double weighted_sum() const
456	0	{
457	0	return m_sum_xiciwi;
458	0	}
459
460		/**
461		* The number of raster cells with any defined value
462		* covered by the polygon. Weights are not taken
463		* into account.
464		*/
465		double count() const
466	0	{
467	0	return m_sum_ci;
468	0	}
469
470		/**
471		* The number of raster cells with a specific value
472		* covered by the polygon. Weights are not taken
473		* into account.
474		*/
475		std::optional<double> count(const ValueType &value) const
476		{
477		const auto &entry = m_freq.find(value);
478
479		if (entry == m_freq.end())
480		{
481		return std::nullopt;
482		}
483
484		return entry->second.m_sum_ci;
485		}
486
487		/**
488		* The fraction of defined raster cells covered by the polygon with
489		* a value that equals the specified value.
490		* Weights are not taken into account.
491		*/
492		std::optional<double> frac(const ValueType &value) const
493		{
494		auto count_for_value = count(value);
495
496		if (!count_for_value.has_value())
497		{
498		return count_for_value;
499		}
500
501		return count_for_value.value() / count();
502		}
503
504		/**
505		* The weighted fraction of defined raster cells covered by the polygon with
506		* a value that equals the specified value.
507		*/
508		std::optional<double> weighted_frac(const ValueType &value) const
509		{
510		auto count_for_value = weighted_count(value);
511
512		if (!count_for_value.has_value())
513		{
514		return count_for_value;
515		}
516
517		return count_for_value.value() / weighted_count();
518		}
519
520		/**
521		* The population variance of raster cells touched
522		* by the polygon. Cell coverage fractions are taken
523		* into account; values of a weighting raster are not.
524		*/
525		double variance() const
526	0	{
527	0	return m_variance.variance();
528	0	}
529
530		/**
531		* The population variance of raster cells touched
532		* by the polygon, taking into account cell coverage
533		* fractions and values of a weighting raster.
534		*/
535		double weighted_variance() const
536	0	{
537	0	return m_weighted_variance.variance();
538	0	}
539
540		/**
541		* The population standard deviation of raster cells
542		* touched by the polygon. Cell coverage fractions
543		* are taken into account; values of a weighting
544		* raster are not.
545		*/
546		double stdev() const
547	0	{
548	0	return m_variance.stdev();
549	0	}
550
551		/**
552		* The population standard deviation of raster cells
553		* touched by the polygon, taking into account cell
554		* coverage fractions and values of a weighting raster.
555		*/
556		double weighted_stdev() const
557	0	{
558	0	return m_weighted_variance.stdev();
559	0	}
560
561		/**
562		* The sum of weights for each cell covered by the
563		* polygon, with each weight multiplied by the coverage
564		* fraction of each cell.
565		*
566		* If any weights are undefined, will return NAN. If this is undesirable,
567		* caller should replace undefined weights with a suitable default
568		* before computing statistics.
569		*/
570		double weighted_count() const
571	0	{
572	0	return m_sum_ciwi;
573	0	}
574
575		/**
576		* The sum of weights for each cell of a specific value covered by the
577		* polygon, with each weight multiplied by the coverage fraction
578		* of each cell.
579		*
580		* If any weights are undefined, will return NAN. If this is undesirable,
581		* caller should replace undefined weights with a suitable default
582		* before computing statistics.
583		*/
584		std::optional<double> weighted_count(const ValueType &value) const
585		{
586		const auto &entry = m_freq.find(value);
587
588		if (entry == m_freq.end())
589		{
590		return std::nullopt;
591		}
592
593		return entry->second.m_sum_ciwi;
594		}
595
596		/**
597		* The raster value occupying the least number of cells
598		* or partial cells within the polygon. When multiple values
599		* cover the same number of cells, the lowest value will
600		* be returned.
601		*
602		* Cell weights are not taken into account.
603		*/
604		std::optional<ValueType> minority() const
605	0	{
606	0	auto it = std::min_element(
607	0	m_freq.cbegin(), m_freq.cend(),
608	0	[](const auto &a, const auto &b)
609	0	{
610	0	return a.second.m_sum_ci < b.second.m_sum_ci \|\|
611	0	(a.second.m_sum_ci == b.second.m_sum_ci &&
612	0	a.first < b.first);
613	0	});
614	0	if (it == m_freq.end())
615	0	{
616	0	return std::nullopt;
617	0	}
618	0	return it->first;
619	0	}
620
621		/**
622		* The number of distinct defined raster values in cells wholly
623		* or partially covered by the polygon.
624		*/
625		std::uint64_t variety() const
626	0	{
627	0	return m_freq.size();
628	0	}
629
630		const std::vector<ValueType> &values() const
631	0	{
632	0	return m_cell_values;
633	0	}
634
635		const std::vector<bool> &values_defined() const
636		{
637		return m_cell_values_defined;
638		}
639
640		const std::vector<float> &coverage_fractions() const
641	0	{
642	0	return m_cell_cov;
643	0	}
644
645		const std::vector<double> &weights() const
646	0	{
647	0	return m_cell_weights;
648	0	}
649
650		const std::vector<bool> &weights_defined() const
651		{
652		return m_cell_weights_defined;
653		}
654
655		const std::vector<double> &center_x() const
656	0	{
657	0	return m_cell_x;
658	0	}
659
660		const std::vector<double> &center_y() const
661	0	{
662	0	return m_cell_y;
663	0	}
664
665		const auto &freq() const
666	0	{
667	0	return m_freq;
668	0	}
669
670		private:
671		ValueType m_min{};
672		ValueType m_max{};
673		std::pair<double, double> m_min_xy{
674		std::numeric_limits<double>::quiet_NaN(),
675		std::numeric_limits<double>::quiet_NaN()};
676		std::pair<double, double> m_max_xy{
677		std::numeric_limits<double>::quiet_NaN(),
678		std::numeric_limits<double>::quiet_NaN()};
679
680		// ci: coverage fraction of pixel i
681		// wi: weight of pixel i
682		// xi: value of pixel i
683		double m_sum_ciwi{0};
684		double m_sum_ci{0};
685		double m_sum_xici{0};
686		double m_sum_xiciwi{0};
687
688		WestVariance m_variance{};
689		WestVariance m_weighted_variance{};
690
691		struct ValueFreqEntry
692		{
693		double m_sum_ci = 0;
694		double m_sum_ciwi = 0;
695		};
696
697		std::unordered_map<ValueType, ValueFreqEntry> m_freq{};
698
699		std::vector<float> m_cell_cov{};
700		std::vector<ValueType> m_cell_values{};
701		std::vector<double> m_cell_weights{};
702		std::vector<double> m_cell_x{};
703		std::vector<double> m_cell_y{};
704		std::vector<bool> m_cell_values_defined{};
705		std::vector<bool> m_cell_weights_defined{};
706
707		RasterStatsOptions m_options;
708		};
709
710		template <typename T>
711		std::ostream &operator<<(std::ostream &os, const RasterStats<T> &stats)
712		{
713		os << "{" << std::endl;
714		os << " \"count\" : " << stats.count() << "," << std::endl;
715
716		os << " \"min\" : ";
717		if (stats.min().has_value())
718		{
719		os << stats.min().value();
720		}
721		else
722		{
723		os << "null";
724		}
725		os << "," << std::endl;
726
727		os << " \"max\" : ";
728		if (stats.max().has_value())
729		{
730		os << stats.max().value();
731		}
732		else
733		{
734		os << "null";
735		}
736		os << "," << std::endl;
737
738		os << " \"mean\" : " << stats.mean() << "," << std::endl;
739		os << " \"sum\" : " << stats.sum() << "," << std::endl;
740		os << " \"weighted_mean\" : " << stats.weighted_mean() << "," << std::endl;
741		os << " \"weighted_sum\" : " << stats.weighted_sum();
742		if (stats.stores_values())
743		{
744		os << "," << std::endl;
745		os << " \"mode\" : ";
746		if (stats.mode().has_value())
747		{
748		os << stats.mode().value();
749		}
750		else
751		{
752		os << "null";
753		}
754		os << "," << std::endl;
755
756		os << " \"minority\" : ";
757		if (stats.minority().has_value())
758		{
759		os << stats.minority().value();
760		}
761		else
762		{
763		os << "null";
764		}
765		os << "," << std::endl;
766
767		os << " \"variety\" : " << stats.variety() << std::endl;
768		}
769		else
770		{
771		os << std::endl;
772		}
773		os << "}" << std::endl;
774		return os;
775		}
776
777		} // namespace gdal
778
779		//! @endcond