Line data Source code
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 : {
43 : return min_coverage_fraction == other.min_coverage_fraction &&
44 : calc_median == other.calc_median &&
45 : calc_variance == other.calc_variance &&
46 : store_histogram == other.store_histogram &&
47 : store_values == other.store_values &&
48 : store_weights == other.store_weights &&
49 : store_coverage_fraction == other.store_coverage_fraction &&
50 : store_xy == other.store_xy &&
51 : include_nodata == other.include_nodata &&
52 : (default_weight == other.default_weight ||
53 : (std::isnan(default_weight) &&
54 : std::isnan(other.default_weight)));
55 : }
56 :
57 : bool operator!=(const RasterStatsOptions &other) const
58 : {
59 : return !(*this == other);
60 : }
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 248 : void process(double x, double w)
83 : {
84 248 : if (w == 0)
85 : {
86 1 : return;
87 : }
88 :
89 247 : double mean_old = mean;
90 :
91 247 : sum_w += w;
92 247 : mean += (w / sum_w) * (x - mean_old);
93 247 : t += w * (x - mean_old) * (x - mean);
94 : }
95 :
96 : /** \brief Return the population variance.
97 : */
98 18 : constexpr double variance() const
99 : {
100 18 : return t / sum_w;
101 : }
102 :
103 : /** \brief Return the population standard deviation
104 : */
105 9 : double stdev() const
106 : {
107 9 : return std::sqrt(variance());
108 : }
109 :
110 : /** \brief Return the population coefficient of variation
111 : */
112 : double coefficent_of_variation() const
113 : {
114 : return stdev() / mean;
115 : }
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 13811 : explicit RasterStats(const RasterStatsOptions &options)
127 13811 : : m_min{std::numeric_limits<ValueType>::max()},
128 13811 : m_max{std::numeric_limits<ValueType>::lowest()}, m_sum_ciwi{0},
129 13811 : m_sum_ci{0}, m_sum_xici{0}, m_sum_xiciwi{0}, m_options{options}
130 : {
131 13811 : }
132 :
133 : // All pixels covered 100%
134 13600 : 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 : {
138 27200 : for (size_t i = 0; i < nX * nY; i++)
139 : {
140 13600 : if (pabyMask[i] == 255)
141 : {
142 13464 : if (padfX && padfY)
143 : {
144 10000 : process_location(padfX[i % nX], padfY[i / nX]);
145 : }
146 24264 : const double dfWeight =
147 : pabyWeightsMask
148 10800 : ? (pabyWeightsMask[i] == 255
149 : ? padfWeights[i]
150 0 : : std::numeric_limits<double>::quiet_NaN())
151 : : 1.0;
152 13464 : process_value(pValues[i], 1.0, dfWeight);
153 : }
154 : }
155 13600 : }
156 :
157 : // Pixels covered 0% or 100%
158 762 : 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 : {
163 14006 : for (size_t i = 0; i < nX * nY; i++)
164 : {
165 13244 : if (pabyMask[i] == 255 && pabyCov[i])
166 : {
167 7125 : if (pdfX && pdfY)
168 : {
169 1298 : process_location(pdfX[i % nX], pdfY[i / nX]);
170 : }
171 7274 : const double dfWeight =
172 : pabyWeightsMask
173 149 : ? (pabyWeightsMask[i] == 255
174 : ? padfWeights[i]
175 1 : : std::numeric_limits<double>::quiet_NaN())
176 : : 1.0;
177 7125 : process_value(pValues[i], 1.0, dfWeight);
178 : }
179 : }
180 762 : }
181 :
182 : // Pixels fractionally covered
183 312 : 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 : {
188 4864 : for (size_t i = 0; i < nX * nY; i++)
189 : {
190 4552 : if (pabyMask[i] == 255 &&
191 4504 : pfCov[i] >= m_options.min_coverage_fraction)
192 : {
193 3080 : if (pdfX && pdfY)
194 : {
195 764 : process_location(pdfX[i % nX], pdfY[i / nX]);
196 : }
197 3488 : const double dfWeight =
198 : pabyWeightsMask
199 408 : ? (pabyWeightsMask[i] == 255
200 : ? padfWeights[i]
201 0 : : std::numeric_limits<double>::quiet_NaN())
202 : : 1.0;
203 3080 : process_value(pValues[i], pfCov[i], dfWeight);
204 : }
205 : }
206 312 : }
207 :
208 12062 : void process_location(double x, double y)
209 : {
210 12062 : if (m_options.store_xy)
211 : {
212 12062 : m_cell_x.push_back(x);
213 12062 : m_cell_y.push_back(y);
214 : }
215 12062 : }
216 :
217 23669 : void process_value(const ValueType &val, float coverage, double weight)
218 : {
219 23669 : if (m_options.store_coverage_fraction)
220 : {
221 16 : m_cell_cov.push_back(coverage);
222 : }
223 :
224 23669 : m_sum_ci += static_cast<double>(coverage);
225 23669 : m_sum_xici += static_cast<double>(val) * static_cast<double>(coverage);
226 :
227 23669 : double ciwi = static_cast<double>(coverage) * weight;
228 23669 : m_sum_ciwi += ciwi;
229 23669 : m_sum_xiciwi += static_cast<double>(val) * ciwi;
230 :
231 23669 : if (m_options.calc_variance)
232 : {
233 124 : m_variance.process(static_cast<double>(val),
234 : static_cast<double>(coverage));
235 124 : m_weighted_variance.process(static_cast<double>(val), ciwi);
236 : }
237 :
238 23669 : if (val < m_min)
239 : {
240 657 : m_min = val;
241 657 : if (m_options.store_xy)
242 : {
243 166 : m_min_xy = {m_cell_x.back(), m_cell_y.back()};
244 : }
245 : }
246 :
247 23669 : if (val > m_max)
248 : {
249 2366 : m_max = val;
250 2366 : if (m_options.store_xy)
251 : {
252 353 : m_max_xy = {m_cell_x.back(), m_cell_y.back()};
253 : }
254 : }
255 :
256 23669 : if (m_options.store_histogram)
257 : {
258 164 : auto &entry = m_freq[val];
259 164 : entry.m_sum_ci += static_cast<double>(coverage);
260 164 : entry.m_sum_ciwi += ciwi;
261 : }
262 :
263 23669 : if (m_options.store_values || m_options.calc_median)
264 : {
265 1616 : m_cell_values.push_back(val);
266 : }
267 :
268 23669 : if (m_options.store_weights)
269 : {
270 36 : m_cell_weights.push_back(weight);
271 : }
272 23669 : }
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 94 : double mean() const
279 : {
280 94 : if (count() > 0)
281 : {
282 92 : return sum() / count();
283 : }
284 : else
285 : {
286 2 : return std::numeric_limits<double>::quiet_NaN();
287 : }
288 : }
289 :
290 : /** The median value of cells touched by this polygon. The cell
291 : * coverage fraction is not taken into account. */
292 20 : double median() const
293 : {
294 20 : auto n = m_cell_values.size();
295 20 : if (n == 0)
296 : {
297 0 : return std::numeric_limits<double>::quiet_NaN();
298 : }
299 20 : if (n == 1)
300 : {
301 4 : return static_cast<double>(m_cell_values[0]);
302 : }
303 16 : if (n == 2)
304 : {
305 4 : return 0.5 * (static_cast<double>(m_cell_values[0]) +
306 4 : static_cast<double>(m_cell_values[1]));
307 : }
308 :
309 12 : std::vector<ValueType> *values = nullptr;
310 12 : std::unique_ptr<std::vector<ValueType>> values_copy;
311 :
312 12 : if (m_options.store_values)
313 : {
314 6 : values_copy =
315 6 : std::make_unique<std::vector<ValueType>>(m_cell_values);
316 6 : values = values_copy.get();
317 : }
318 : else
319 : {
320 6 : values = const_cast<std::vector<ValueType> *>(&m_cell_values);
321 : }
322 :
323 12 : auto mid = std::next(values->begin(), n / 2);
324 12 : std::nth_element(values->begin(), mid, values->end());
325 :
326 12 : if (n % 2 != 0)
327 : {
328 4 : return static_cast<double>(*mid);
329 : }
330 :
331 8 : auto lhs_max = std::max_element(values->begin(), mid);
332 :
333 : return 0.5 *
334 8 : (static_cast<double>(*lhs_max) + static_cast<double>(*mid));
335 : }
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 43 : double weighted_mean() const
347 : {
348 43 : if (weighted_count() > 0)
349 : {
350 28 : return weighted_sum() / weighted_count();
351 : }
352 : else
353 : {
354 15 : return std::numeric_limits<double>::quiet_NaN();
355 : }
356 : }
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 46 : std::optional<ValueType> mode() const
374 : {
375 46 : auto it = std::max_element(
376 : m_freq.cbegin(), m_freq.cend(),
377 14 : [](const auto &a, const auto &b)
378 : {
379 24 : return a.second.m_sum_ci < b.second.m_sum_ci ||
380 10 : (a.second.m_sum_ci == b.second.m_sum_ci &&
381 18 : a.first < b.first);
382 : });
383 46 : if (it == m_freq.end())
384 : {
385 39 : return std::nullopt;
386 : }
387 7 : return it->first;
388 : }
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 2 : std::optional<ValueType> min() const
395 : {
396 2 : if (m_sum_ci == 0)
397 : {
398 0 : return std::nullopt;
399 : }
400 2 : return m_min;
401 : }
402 :
403 : /// XY values corresponding to the center of the cell whose value
404 : /// is returned by min()
405 4 : std::optional<std::pair<double, double>> min_xy() const
406 : {
407 4 : if (m_sum_ci == 0)
408 : {
409 0 : return std::nullopt;
410 : }
411 4 : return m_min_xy;
412 : }
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 34 : std::optional<ValueType> max() const
419 : {
420 34 : if (m_sum_ci == 0)
421 : {
422 0 : return std::nullopt;
423 : }
424 34 : return m_max;
425 : }
426 :
427 : /// XY values corresponding to the center of the cell whose value
428 : /// is returned by max()
429 68 : std::optional<std::pair<double, double>> max_xy() const
430 : {
431 68 : if (m_sum_ci == 0)
432 : {
433 0 : return std::nullopt;
434 : }
435 68 : return m_max_xy;
436 : }
437 :
438 : /**
439 : * The sum of raster cells covered by the polygon, with each raster
440 : * value weighted by its coverage fraction.
441 : */
442 334 : double sum() const
443 : {
444 334 : return m_sum_xici;
445 : }
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 52 : double weighted_sum() const
456 : {
457 52 : return m_sum_xiciwi;
458 : }
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 245 : double count() const
466 : {
467 245 : return m_sum_ci;
468 : }
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 2 : double variance() const
526 : {
527 2 : return m_variance.variance();
528 : }
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 7 : double weighted_variance() const
536 : {
537 7 : return m_weighted_variance.variance();
538 : }
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 2 : double stdev() const
547 : {
548 2 : return m_variance.stdev();
549 : }
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 7 : double weighted_stdev() const
557 : {
558 7 : return m_weighted_variance.stdev();
559 : }
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 71 : double weighted_count() const
571 : {
572 71 : return m_sum_ciwi;
573 : }
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 2 : std::optional<ValueType> minority() const
605 : {
606 2 : auto it = std::min_element(
607 : m_freq.cbegin(), m_freq.cend(),
608 14 : [](const auto &a, const auto &b)
609 : {
610 28 : return a.second.m_sum_ci < b.second.m_sum_ci ||
611 14 : (a.second.m_sum_ci == b.second.m_sum_ci &&
612 20 : a.first < b.first);
613 : });
614 2 : if (it == m_freq.end())
615 : {
616 0 : return std::nullopt;
617 : }
618 2 : return it->first;
619 : }
620 :
621 : /**
622 : * The number of distinct defined raster values in cells wholly
623 : * or partially covered by the polygon.
624 : */
625 2 : std::uint64_t variety() const
626 : {
627 2 : return m_freq.size();
628 : }
629 :
630 12 : const std::vector<ValueType> &values() const
631 : {
632 12 : return m_cell_values;
633 : }
634 :
635 : const std::vector<bool> &values_defined() const
636 : {
637 : return m_cell_values_defined;
638 : }
639 :
640 2 : const std::vector<float> &coverage_fractions() const
641 : {
642 2 : return m_cell_cov;
643 : }
644 :
645 5 : const std::vector<double> &weights() const
646 : {
647 5 : return m_cell_weights;
648 : }
649 :
650 : const std::vector<bool> &weights_defined() const
651 : {
652 : return m_cell_weights_defined;
653 : }
654 :
655 10 : const std::vector<double> ¢er_x() const
656 : {
657 10 : return m_cell_x;
658 : }
659 :
660 10 : const std::vector<double> ¢er_y() const
661 : {
662 10 : return m_cell_y;
663 : }
664 :
665 4 : const auto &freq() const
666 : {
667 4 : return m_freq;
668 : }
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
|
