Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: GDALZonalStats implementation
5 : * Author: Dan Baston
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2025, ISciences LLC
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "cpl_string.h"
14 : #include "gdal_priv.h"
15 : #include "gdal_alg.h"
16 : #include "gdal_utils.h"
17 : #include "ogrsf_frmts.h"
18 : #include "raster_stats.h"
19 :
20 : #include "../frmts/mem/memdataset.h"
21 : #include "../frmts/vrt/vrtdataset.h"
22 :
23 : #include "ogr_geos.h"
24 :
25 : #include <algorithm>
26 : #include <array>
27 : #include <cstring>
28 : #include <limits>
29 : #include <variant>
30 : #include <vector>
31 :
32 : #if GEOS_VERSION_MAJOR > 3 || \
33 : (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 14)
34 : #define GEOS_GRID_INTERSECTION_AVAILABLE 1
35 : #endif
36 :
37 : struct GDALZonalStatsOptions
38 : {
39 147 : CPLErr Init(CSLConstList papszOptions)
40 : {
41 928 : for (const auto &[key, value] : cpl::IterateNameValue(papszOptions))
42 : {
43 781 : if (EQUAL(key, "BANDS"))
44 : {
45 : const CPLStringList aosBands(CSLTokenizeString2(
46 21 : value, ",", CSLT_STRIPLEADSPACES | CSLT_STRIPENDSPACES));
47 49 : for (const char *pszBand : aosBands)
48 : {
49 28 : int nBand = std::atoi(pszBand);
50 28 : if (nBand <= 0)
51 : {
52 0 : CPLError(CE_Failure, CPLE_IllegalArg,
53 : "Invalid band: %s", pszBand);
54 0 : return CE_Failure;
55 : }
56 28 : bands.push_back(nBand);
57 : }
58 : }
59 760 : else if (EQUAL(key, "INCLUDE_FIELDS"))
60 : {
61 : CPLStringList aosFields(CSLTokenizeString2(
62 : value, ",",
63 : CSLT_HONOURSTRINGS | CSLT_STRIPLEADSPACES |
64 10 : CSLT_STRIPENDSPACES));
65 12 : for (const char *pszField : aosFields)
66 : {
67 7 : include_fields.push_back(pszField);
68 : }
69 : }
70 755 : else if (EQUAL(key, "PIXEL_INTERSECTION"))
71 : {
72 147 : if (EQUAL(value, "DEFAULT"))
73 : {
74 83 : pixels = DEFAULT;
75 : }
76 64 : else if (EQUAL(value, "ALL-TOUCHED") ||
77 56 : EQUAL(value, "ALL_TOUCHED"))
78 : {
79 8 : pixels = ALL_TOUCHED;
80 : }
81 56 : else if (EQUAL(value, "FRACTIONAL"))
82 : {
83 56 : pixels = FRACTIONAL;
84 : }
85 : else
86 : {
87 0 : CPLError(CE_Failure, CPLE_IllegalArg,
88 : "Unexpected value of PIXEL_INTERSECTION: %s",
89 : value);
90 0 : return CE_Failure;
91 : }
92 : }
93 608 : else if (EQUAL(key, "RASTER_CHUNK_SIZE_BYTES"))
94 : {
95 147 : char *endptr = nullptr;
96 147 : errno = 0;
97 147 : const auto memory64 = std::strtoull(value, &endptr, 10);
98 294 : bool ok = errno != ERANGE && memory64 != ULLONG_MAX &&
99 147 : endptr == value + strlen(value);
100 : if constexpr (sizeof(memory64) > sizeof(size_t))
101 : {
102 : ok = ok &&
103 : memory64 <= std::numeric_limits<size_t>::max() - 1;
104 : }
105 147 : if (!ok)
106 : {
107 0 : CPLError(CE_Failure, CPLE_IllegalArg,
108 : "Invalid memory size: %s", value);
109 0 : return CE_Failure;
110 : }
111 147 : memory = static_cast<size_t>(memory64);
112 : }
113 461 : else if (EQUAL(key, "STATS"))
114 : {
115 294 : stats = CPLStringList(CSLTokenizeString2(
116 147 : value, ",", CSLT_STRIPLEADSPACES | CSLT_STRIPENDSPACES));
117 : }
118 314 : else if (EQUAL(key, "STRATEGY"))
119 : {
120 147 : if (EQUAL(value, "FEATURE_SEQUENTIAL"))
121 : {
122 98 : strategy = FEATURE_SEQUENTIAL;
123 : }
124 49 : else if (EQUAL(value, "RASTER_SEQUENTIAL"))
125 : {
126 49 : strategy = RASTER_SEQUENTIAL;
127 : }
128 : else
129 : {
130 0 : CPLError(CE_Failure, CPLE_IllegalArg,
131 : "Unexpected value of STRATEGY: %s", value);
132 0 : return CE_Failure;
133 : }
134 : }
135 167 : else if (EQUAL(key, "WEIGHTS_BAND"))
136 : {
137 147 : weights_band = std::atoi(value);
138 147 : if (weights_band <= 0)
139 : {
140 0 : CPLError(CE_Failure, CPLE_IllegalArg,
141 : "Invalid weights band: %s", value);
142 0 : return CE_Failure;
143 : }
144 : }
145 20 : else if (EQUAL(key, "ZONES_BAND"))
146 : {
147 1 : zones_band = std::atoi(value);
148 1 : if (zones_band <= 0)
149 : {
150 0 : CPLError(CE_Failure, CPLE_IllegalArg,
151 : "Invalid zones band: %s", value);
152 0 : return CE_Failure;
153 : }
154 : }
155 19 : else if (EQUAL(key, "ZONES_LAYER"))
156 : {
157 1 : zones_layer = value;
158 : }
159 18 : else if (STARTS_WITH(key, "LCO_"))
160 : {
161 18 : layer_creation_options.SetNameValue(key + strlen("LCO_"),
162 18 : value);
163 : }
164 : else
165 : {
166 0 : CPLError(CE_Failure, CPLE_IllegalArg,
167 : "Unexpected zonal stats option: %s", key);
168 : }
169 : }
170 :
171 147 : return CE_None;
172 : }
173 :
174 : enum PixelIntersection
175 : {
176 : DEFAULT,
177 : ALL_TOUCHED,
178 : FRACTIONAL,
179 : };
180 :
181 : enum Strategy
182 : {
183 : FEATURE_SEQUENTIAL,
184 : RASTER_SEQUENTIAL,
185 : };
186 :
187 : PixelIntersection pixels{DEFAULT};
188 : Strategy strategy{FEATURE_SEQUENTIAL};
189 : std::vector<std::string> stats{};
190 : std::vector<std::string> include_fields{};
191 : std::vector<int> bands{};
192 : std::string zones_layer{};
193 : std::size_t memory{0};
194 : int zones_band{};
195 : int weights_band{};
196 : CPLStringList layer_creation_options{};
197 : };
198 :
199 16 : template <typename T = GByte> auto CreateBuffer()
200 : {
201 16 : return std::unique_ptr<T, VSIFreeReleaser>(nullptr);
202 : }
203 :
204 : template <typename T>
205 533 : void Realloc(T &buf, size_t size1, size_t size2, bool &success)
206 : {
207 533 : if (!success)
208 : {
209 0 : return;
210 : }
211 : if constexpr (sizeof(size_t) < sizeof(uint64_t))
212 : {
213 : if (size1 > std::numeric_limits<size_t>::max() / size2)
214 : {
215 : success = false;
216 : CPLError(CE_Failure, CPLE_OutOfMemory,
217 : "Too big memory allocation attempt");
218 : return;
219 : }
220 : }
221 533 : const auto size = size1 * size2;
222 533 : auto oldBuf = buf.release();
223 : auto newBuf = static_cast<typename T::element_type *>(
224 533 : VSI_REALLOC_VERBOSE(oldBuf, size));
225 533 : if (newBuf == nullptr)
226 : {
227 0 : VSIFree(oldBuf);
228 0 : success = false;
229 : }
230 533 : buf.reset(newBuf);
231 : }
232 :
233 57 : static void CalculateCellCenters(const GDALRasterWindow &window,
234 : const GDALGeoTransform >, double *padfX,
235 : double *padfY)
236 : {
237 : double dfJunk;
238 57 : double x0 = window.nXOff;
239 57 : double y0 = window.nYOff;
240 :
241 765 : for (int i = 0; i < window.nXSize; i++)
242 : {
243 708 : gt.Apply(x0 + i + 0.5, window.nYOff, padfX + i, &dfJunk);
244 : }
245 1176 : for (int i = 0; i < window.nYSize; i++)
246 : {
247 1119 : gt.Apply(x0, y0 + i + 0.5, &dfJunk, padfY + i);
248 : }
249 57 : }
250 :
251 : class GDALZonalStatsImpl
252 : {
253 : public:
254 : enum Stat
255 : {
256 : CENTER_X, // must be first value
257 : CENTER_Y,
258 : COUNT,
259 : COVERAGE,
260 : FRAC,
261 : MAX,
262 : MAX_CENTER_X,
263 : MAX_CENTER_Y,
264 : MEAN,
265 : MIN,
266 : MIN_CENTER_X,
267 : MIN_CENTER_Y,
268 : MINORITY,
269 : MODE,
270 : STDEV,
271 : SUM,
272 : UNIQUE,
273 : VALUES,
274 : VARIANCE,
275 : VARIETY,
276 : WEIGHTED_FRAC,
277 : WEIGHTED_MEAN,
278 : WEIGHTED_SUM,
279 : WEIGHTED_STDEV,
280 : WEIGHTED_VARIANCE,
281 : WEIGHTS,
282 : INVALID, // must be last value
283 : };
284 :
285 134 : static constexpr bool IsWeighted(Stat eStat)
286 : {
287 133 : return eStat == WEIGHTS || eStat == WEIGHTED_FRAC ||
288 132 : eStat == WEIGHTED_MEAN || eStat == WEIGHTED_SUM ||
289 267 : eStat == WEIGHTED_VARIANCE || eStat == WEIGHTED_STDEV;
290 : }
291 :
292 : using BandOrLayer = std::variant<GDALRasterBand *, OGRLayer *>;
293 :
294 144 : GDALZonalStatsImpl(GDALDataset &src, GDALDataset &dst, GDALDataset *weights,
295 : BandOrLayer zones, const GDALZonalStatsOptions &options)
296 144 : : m_src(src), m_weights(weights), m_dst(dst), m_zones(zones),
297 144 : m_coverageDataType(options.pixels == GDALZonalStatsOptions::FRACTIONAL
298 144 : ? GDT_Float32
299 : : GDT_UInt8),
300 : m_options(options),
301 288 : m_maxCells(options.memory /
302 288 : std::max(1, GDALGetDataTypeSizeBytes(m_workingDataType)))
303 : {
304 : #ifdef HAVE_GEOS
305 144 : m_geosContext = OGRGeometry::createGEOSContext();
306 : #endif
307 144 : }
308 :
309 144 : ~GDALZonalStatsImpl()
310 144 : {
311 : #ifdef HAVE_GEOS
312 144 : if (m_geosContext)
313 : {
314 144 : finishGEOS_r(m_geosContext);
315 : }
316 : #endif
317 144 : }
318 :
319 : private:
320 144 : bool Init()
321 : {
322 : #if !(GEOS_GRID_INTERSECTION_AVAILABLE)
323 : if (m_options.pixels == GDALZonalStatsOptions::FRACTIONAL)
324 : {
325 : CPLError(CE_Failure, CPLE_AppDefined,
326 : "Fractional pixel coverage calculation requires a GDAL "
327 : "build against GEOS >= 3.14");
328 : return false;
329 : }
330 : #endif
331 :
332 144 : if (m_options.bands.empty())
333 : {
334 123 : const int nBands = m_src.GetRasterCount();
335 123 : if (nBands == 0)
336 : {
337 0 : CPLError(CE_Failure, CPLE_AppDefined,
338 : "GDALRasterZonalStats: input dataset has no bands");
339 0 : return false;
340 : }
341 123 : m_options.bands.resize(nBands);
342 246 : for (int i = 0; i < nBands; i++)
343 : {
344 123 : m_options.bands[i] = i + 1;
345 : }
346 : }
347 : else
348 : {
349 49 : for (int nBand : m_options.bands)
350 : {
351 28 : if (nBand <= 0 || nBand > m_src.GetRasterCount())
352 : {
353 0 : CPLError(CE_Failure, CPLE_AppDefined,
354 : "GDALRasterZonalStats: Invalid band number: %d",
355 : nBand);
356 0 : return false;
357 : }
358 : }
359 : }
360 :
361 : {
362 144 : const auto eSrcType = m_src.GetRasterBand(m_options.bands.front())
363 144 : ->GetRasterDataType();
364 144 : if (GDALDataTypeIsConversionLossy(eSrcType, m_workingDataType))
365 : {
366 5 : CPLError(CE_Failure, CPLE_AppDefined,
367 : "GDALRasterZonalStats: Source data type %s is not "
368 : "supported",
369 : GDALGetDataTypeName(eSrcType));
370 5 : return false;
371 : }
372 : }
373 :
374 139 : if (m_weights)
375 : {
376 79 : if (m_options.weights_band > m_weights->GetRasterCount())
377 : {
378 1 : CPLError(CE_Failure, CPLE_AppDefined,
379 : "GDALRasterZonalStats: invalid weights band");
380 1 : return false;
381 : }
382 : const auto eWeightsType =
383 78 : m_weights->GetRasterBand(m_options.weights_band)
384 78 : ->GetRasterDataType();
385 78 : if (GDALDataTypeIsConversionLossy(eWeightsType, GDT_Float64))
386 : {
387 5 : CPLError(CE_Failure, CPLE_AppDefined,
388 : "GDALRasterZonalStats: Weights data type %s is not "
389 : "supported",
390 : GDALGetDataTypeName(eWeightsType));
391 5 : return false;
392 : }
393 : }
394 :
395 374 : for (const auto &stat : m_options.stats)
396 : {
397 246 : const auto eStat = GetStat(stat);
398 246 : switch (eStat)
399 : {
400 0 : case INVALID:
401 : {
402 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid stat: %s",
403 : stat.c_str());
404 5 : return false;
405 : }
406 :
407 2 : case COVERAGE:
408 2 : m_stats_options.store_coverage_fraction = true;
409 2 : break;
410 :
411 30 : case VARIETY:
412 : case MODE:
413 : case MINORITY:
414 : case UNIQUE:
415 : case FRAC:
416 : case WEIGHTED_FRAC:
417 30 : m_stats_options.store_histogram = true;
418 30 : break;
419 :
420 20 : case VARIANCE:
421 : case STDEV:
422 : case WEIGHTED_VARIANCE:
423 : case WEIGHTED_STDEV:
424 20 : m_stats_options.calc_variance = true;
425 20 : break;
426 :
427 30 : case CENTER_X:
428 : case CENTER_Y:
429 : case MIN_CENTER_X:
430 : case MIN_CENTER_Y:
431 : case MAX_CENTER_X:
432 : case MAX_CENTER_Y:
433 30 : m_stats_options.store_xy = true;
434 30 : break;
435 :
436 2 : case VALUES:
437 2 : m_stats_options.store_values = true;
438 2 : break;
439 :
440 6 : case WEIGHTS:
441 6 : m_stats_options.store_weights = true;
442 6 : break;
443 :
444 156 : case COUNT:
445 : case MIN:
446 : case MAX:
447 : case SUM:
448 : case MEAN:
449 : case WEIGHTED_SUM:
450 : case WEIGHTED_MEAN:
451 156 : break;
452 : }
453 246 : if (m_weights == nullptr && IsWeighted(eStat))
454 : {
455 5 : CPLError(CE_Failure, CPLE_AppDefined,
456 : "Stat %s requires weights but none were provided",
457 : stat.c_str());
458 5 : return false;
459 : }
460 : }
461 :
462 128 : if (m_src.GetGeoTransform(m_srcGT) != CE_None)
463 : {
464 1 : CPLError(CE_Failure, CPLE_AppDefined,
465 : "Dataset has no geotransform");
466 1 : return false;
467 : }
468 127 : if (!m_srcGT.GetInverse(m_srcInvGT))
469 : {
470 1 : CPLError(CE_Failure, CPLE_AppDefined,
471 : "Dataset geotransform cannot be inverted");
472 1 : return false;
473 : }
474 :
475 126 : const OGRSpatialReference *poRastSRS = m_src.GetSpatialRefRasterOnly();
476 : const OGRSpatialReference *poWeightsSRS =
477 126 : m_weights ? m_weights->GetSpatialRefRasterOnly() : nullptr;
478 126 : const OGRSpatialReference *poZonesSRS = nullptr;
479 :
480 126 : if (ZonesAreFeature())
481 : {
482 109 : const OGRLayer *poSrcLayer = std::get<OGRLayer *>(m_zones);
483 109 : const OGRFeatureDefn *poSrcDefn = poSrcLayer->GetLayerDefn();
484 109 : poZonesSRS = poSrcLayer->GetSpatialRef();
485 :
486 113 : for (const auto &field : m_options.include_fields)
487 : {
488 6 : if (poSrcDefn->GetFieldIndex(field.c_str()) == -1)
489 : {
490 2 : CPLError(CE_Failure, CPLE_AppDefined, "Field %s not found.",
491 : field.c_str());
492 2 : return false;
493 : }
494 : }
495 : }
496 : else
497 : {
498 17 : poZonesSRS = std::get<GDALRasterBand *>(m_zones)
499 17 : ->GetDataset()
500 17 : ->GetSpatialRefRasterOnly();
501 :
502 17 : if (!m_options.include_fields.empty())
503 : {
504 1 : CPLError(CE_Failure, CPLE_AppDefined,
505 : "Cannot include fields from raster zones");
506 1 : return false;
507 : }
508 : }
509 :
510 123 : CPLStringList aosOptions;
511 123 : aosOptions.AddNameValue("IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING", "1");
512 :
513 135 : if (poRastSRS && poZonesSRS &&
514 12 : !poRastSRS->IsSame(poZonesSRS, aosOptions.List()))
515 : {
516 2 : CPLError(CE_Warning, CPLE_AppDefined,
517 : "Inputs and zones do not have the same SRS");
518 : }
519 :
520 125 : if (poWeightsSRS && poZonesSRS &&
521 2 : !poWeightsSRS->IsSame(poZonesSRS, aosOptions.List()))
522 : {
523 2 : CPLError(CE_Warning, CPLE_AppDefined,
524 : "Weights and zones do not have the same SRS");
525 : }
526 :
527 127 : if (poWeightsSRS && poRastSRS &&
528 4 : !poWeightsSRS->IsSame(poRastSRS, aosOptions.List()))
529 : {
530 4 : CPLError(CE_Warning, CPLE_AppDefined,
531 : "Inputs and weights do not have the same SRS");
532 : }
533 :
534 123 : return true;
535 : }
536 :
537 11761 : gdal::RasterStats<double> CreateStats() const
538 : {
539 11761 : return gdal::RasterStats<double>{m_stats_options};
540 : }
541 :
542 122 : OGRLayer *GetOutputLayer(bool createValueField)
543 : {
544 244 : std::string osLayerName = "stats";
545 :
546 : OGRLayer *poLayer =
547 122 : m_dst.CreateLayer(osLayerName.c_str(), nullptr,
548 122 : m_options.layer_creation_options.List());
549 122 : if (!poLayer)
550 0 : return nullptr;
551 :
552 122 : if (createValueField)
553 : {
554 16 : OGRFieldDefn oFieldDefn("value", OFTReal);
555 16 : if (poLayer->CreateField(&oFieldDefn) != OGRERR_NONE)
556 0 : return nullptr;
557 : }
558 :
559 122 : if (!m_options.include_fields.empty())
560 : {
561 : const OGRFeatureDefn *poSrcDefn =
562 2 : std::get<OGRLayer *>(m_zones)->GetLayerDefn();
563 :
564 6 : for (const auto &field : m_options.include_fields)
565 : {
566 4 : const int iField = poSrcDefn->GetFieldIndex(field.c_str());
567 : // Already checked field names during Init()
568 4 : if (poLayer->CreateField(poSrcDefn->GetFieldDefn(iField)) !=
569 : OGRERR_NONE)
570 0 : return nullptr;
571 : }
572 : }
573 :
574 251 : for (int iBand : m_options.bands)
575 : {
576 129 : auto &aiStatFields = m_statFields[iBand];
577 129 : aiStatFields.fill(-1);
578 :
579 382 : for (const auto &stat : m_options.stats)
580 : {
581 253 : const Stat eStat = GetStat(stat);
582 :
583 253 : std::string osFieldName;
584 253 : if (m_options.bands.size() > 1)
585 : {
586 36 : osFieldName = CPLSPrintf("%s_band_%d", stat.c_str(), iBand);
587 : }
588 : else
589 : {
590 217 : osFieldName = stat;
591 : }
592 :
593 : OGRFieldDefn oFieldDefn(osFieldName.c_str(),
594 253 : GetFieldType(eStat));
595 253 : if (poLayer->CreateField(&oFieldDefn) != OGRERR_NONE)
596 0 : return nullptr;
597 : const int iNewField =
598 253 : poLayer->GetLayerDefn()->GetFieldIndex(osFieldName.c_str());
599 253 : aiStatFields[eStat] = iNewField;
600 : }
601 : }
602 :
603 122 : return poLayer;
604 : }
605 :
606 6812 : static const char *GetString(Stat s)
607 : {
608 6812 : switch (s)
609 : {
610 499 : case CENTER_X:
611 499 : return "center_x";
612 495 : case CENTER_Y:
613 495 : return "center_y";
614 491 : case COUNT:
615 491 : return "count";
616 447 : case COVERAGE:
617 447 : return "coverage";
618 443 : case FRAC:
619 443 : return "frac";
620 439 : case MAX:
621 439 : return "max";
622 416 : case MAX_CENTER_X:
623 416 : return "max_center_x";
624 393 : case MAX_CENTER_Y:
625 393 : return "max_center_y";
626 370 : case MEAN:
627 370 : return "mean";
628 305 : case MIN:
629 305 : return "min";
630 301 : case MIN_CENTER_X:
631 301 : return "min_center_x";
632 297 : case MIN_CENTER_Y:
633 297 : return "min_center_y";
634 293 : case MINORITY:
635 293 : return "minority";
636 289 : case MODE:
637 289 : return "mode";
638 245 : case STDEV:
639 245 : return "stdev";
640 241 : case SUM:
641 241 : return "sum";
642 128 : case UNIQUE:
643 128 : return "unique";
644 124 : case VALUES:
645 124 : return "values";
646 120 : case VARIANCE:
647 120 : return "variance";
648 116 : case VARIETY:
649 116 : return "variety";
650 112 : case WEIGHTED_FRAC:
651 112 : return "weighted_frac";
652 112 : case WEIGHTED_MEAN:
653 112 : return "weighted_mean";
654 58 : case WEIGHTED_SUM:
655 58 : return "weighted_sum";
656 41 : case WEIGHTED_STDEV:
657 41 : return "weighted_stdev";
658 26 : case WEIGHTED_VARIANCE:
659 26 : return "weighted_variance";
660 11 : case WEIGHTS:
661 11 : return "weights";
662 0 : case INVALID:
663 0 : break;
664 : }
665 0 : return "invalid";
666 : }
667 :
668 499 : static Stat GetStat(const std::string &stat)
669 : {
670 6812 : for (Stat s = CENTER_X; s < INVALID; s = static_cast<Stat>(s + 1))
671 : {
672 6812 : if (stat == GetString(s))
673 499 : return s;
674 : }
675 0 : return INVALID;
676 : }
677 :
678 253 : static OGRFieldType GetFieldType(Stat stat)
679 : {
680 253 : switch (stat)
681 : {
682 17 : case CENTER_X:
683 : case CENTER_Y:
684 : case COVERAGE:
685 : case FRAC:
686 : case UNIQUE:
687 : case VALUES:
688 : case WEIGHTS:
689 17 : return OFTRealList;
690 2 : case VARIETY:
691 2 : return OFTInteger;
692 234 : case COUNT:
693 : case MAX:
694 : case MAX_CENTER_X:
695 : case MAX_CENTER_Y:
696 : case MEAN:
697 : case MIN:
698 : case MIN_CENTER_X:
699 : case MIN_CENTER_Y:
700 : case MINORITY:
701 : case MODE:
702 : case STDEV:
703 : case SUM:
704 : case VARIANCE:
705 : case WEIGHTED_FRAC:
706 : case WEIGHTED_MEAN:
707 : case WEIGHTED_SUM:
708 : case WEIGHTED_STDEV:
709 : case WEIGHTED_VARIANCE:
710 : case INVALID:
711 234 : break;
712 : }
713 234 : return OFTReal;
714 : }
715 :
716 7182 : int GetFieldIndex(int iBand, Stat eStat) const
717 : {
718 7182 : auto it = m_statFields.find(iBand);
719 7182 : if (it == m_statFields.end())
720 : {
721 0 : return -1;
722 : }
723 :
724 7182 : return it->second[eStat];
725 : }
726 :
727 382 : OGREnvelope ToEnvelope(const GDALRasterWindow &window) const
728 : {
729 382 : OGREnvelope oSnappedGeomExtent;
730 382 : m_srcGT.Apply(window, oSnappedGeomExtent);
731 382 : return oSnappedGeomExtent;
732 : }
733 :
734 266 : void SetStatFields(OGRFeature &feature, int iBand,
735 : const gdal::RasterStats<double> &stats) const
736 : {
737 266 : if (auto iField = GetFieldIndex(iBand, CENTER_X); iField != -1)
738 : {
739 2 : const auto ¢er_x = stats.center_x();
740 2 : feature.SetField(iField, static_cast<int>(center_x.size()),
741 : center_x.data());
742 : }
743 266 : if (auto iField = GetFieldIndex(iBand, CENTER_Y); iField != -1)
744 : {
745 2 : const auto ¢er_y = stats.center_y();
746 2 : feature.SetField(iField, static_cast<int>(center_y.size()),
747 : center_y.data());
748 : }
749 266 : if (auto iField = GetFieldIndex(iBand, COUNT); iField != -1)
750 : {
751 46 : feature.SetField(iField, stats.count());
752 : }
753 266 : if (auto iField = GetFieldIndex(iBand, COVERAGE); iField != -1)
754 : {
755 2 : const auto &cov = stats.coverage_fractions();
756 4 : std::vector<double> doubleCov(cov.begin(), cov.end());
757 : // TODO: Add float* overload to Feature::SetField to avoid this copy
758 2 : feature.SetField(iField, static_cast<int>(doubleCov.size()),
759 2 : doubleCov.data());
760 : }
761 266 : if (auto iField = GetFieldIndex(iBand, FRAC); iField != -1)
762 : {
763 2 : const auto count = stats.count();
764 2 : const auto &freq = stats.freq();
765 4 : std::vector<double> values;
766 2 : values.reserve(freq.size());
767 18 : for (const auto &[_, valueCount] : freq)
768 : {
769 16 : values.push_back(valueCount.m_sum_ci / count);
770 : }
771 2 : feature.SetField(iField, static_cast<int>(values.size()),
772 2 : values.data());
773 : }
774 266 : if (auto iField = GetFieldIndex(iBand, MAX); iField != -1)
775 : {
776 34 : const auto &max = stats.max();
777 34 : if (max.has_value())
778 34 : feature.SetField(iField, max.value());
779 : }
780 266 : if (auto iField = GetFieldIndex(iBand, MAX_CENTER_X); iField != -1)
781 : {
782 34 : const auto &loc = stats.max_xy();
783 34 : if (loc.has_value())
784 34 : feature.SetField(iField, loc.value().first);
785 : }
786 266 : if (auto iField = GetFieldIndex(iBand, MAX_CENTER_Y); iField != -1)
787 : {
788 34 : const auto &loc = stats.max_xy();
789 34 : if (loc.has_value())
790 34 : feature.SetField(iField, loc.value().second);
791 : }
792 266 : if (auto iField = GetFieldIndex(iBand, MEAN); iField != -1)
793 : {
794 94 : feature.SetField(iField, stats.mean());
795 : }
796 266 : if (auto iField = GetFieldIndex(iBand, MIN); iField != -1)
797 : {
798 2 : const auto &min = stats.min();
799 2 : if (min.has_value())
800 2 : feature.SetField(iField, min.value());
801 : }
802 266 : if (auto iField = GetFieldIndex(iBand, MINORITY); iField != -1)
803 : {
804 2 : const auto &minority = stats.minority();
805 2 : if (minority.has_value())
806 2 : feature.SetField(iField, minority.value());
807 : }
808 266 : if (auto iField = GetFieldIndex(iBand, MIN_CENTER_X); iField != -1)
809 : {
810 2 : const auto &loc = stats.min_xy();
811 2 : if (loc.has_value())
812 2 : feature.SetField(iField, loc.value().first);
813 : }
814 266 : if (auto iField = GetFieldIndex(iBand, MIN_CENTER_Y); iField != -1)
815 : {
816 2 : const auto &loc = stats.min_xy();
817 2 : if (loc.has_value())
818 2 : feature.SetField(iField, loc.value().second);
819 : }
820 266 : if (auto iField = GetFieldIndex(iBand, MODE); iField != -1)
821 : {
822 46 : const auto &mode = stats.mode();
823 46 : if (mode.has_value())
824 7 : feature.SetField(iField, mode.value());
825 : }
826 266 : if (auto iField = GetFieldIndex(iBand, STDEV); iField != -1)
827 : {
828 2 : feature.SetField(iField, stats.stdev());
829 : }
830 266 : if (auto iField = GetFieldIndex(iBand, SUM); iField != -1)
831 : {
832 172 : feature.SetField(iField, stats.sum());
833 : }
834 266 : if (auto iField = GetFieldIndex(iBand, UNIQUE); iField != -1)
835 : {
836 2 : const auto &freq = stats.freq();
837 4 : std::vector<double> values;
838 2 : values.reserve(freq.size());
839 18 : for (const auto &[value, _] : freq)
840 : {
841 16 : values.push_back(value);
842 : }
843 :
844 2 : feature.SetField(iField, static_cast<int>(values.size()),
845 2 : values.data());
846 : }
847 266 : if (auto iField = GetFieldIndex(iBand, VALUES); iField != -1)
848 : {
849 2 : const auto &values = stats.values();
850 2 : feature.SetField(iField, static_cast<int>(values.size()),
851 : values.data());
852 : }
853 266 : if (auto iField = GetFieldIndex(iBand, VARIANCE); iField != -1)
854 : {
855 2 : feature.SetField(iField, stats.variance());
856 : }
857 266 : if (auto iField = GetFieldIndex(iBand, VARIETY); iField != -1)
858 : {
859 2 : feature.SetField(iField, static_cast<GIntBig>(stats.variety()));
860 : }
861 266 : if (auto iField = GetFieldIndex(iBand, WEIGHTED_FRAC); iField != -1)
862 : {
863 0 : const auto count = stats.count();
864 0 : const auto &freq = stats.freq();
865 0 : std::vector<double> values;
866 0 : values.reserve(freq.size());
867 0 : for (const auto &[_, valueCount] : freq)
868 : {
869 : // Add std::numeric_limits<double>::min() to please Coverity Scan
870 0 : values.push_back(valueCount.m_sum_ciwi /
871 0 : (count + std::numeric_limits<double>::min()));
872 : }
873 0 : feature.SetField(iField, static_cast<int>(values.size()),
874 0 : values.data());
875 : }
876 266 : if (auto iField = GetFieldIndex(iBand, WEIGHTED_MEAN); iField != -1)
877 : {
878 43 : feature.SetField(iField, stats.weighted_mean());
879 : }
880 266 : if (auto iField = GetFieldIndex(iBand, WEIGHTED_STDEV); iField != -1)
881 : {
882 7 : feature.SetField(iField, stats.weighted_stdev());
883 : }
884 266 : if (auto iField = GetFieldIndex(iBand, WEIGHTED_SUM); iField != -1)
885 : {
886 12 : feature.SetField(iField, stats.weighted_sum());
887 : }
888 266 : if (auto iField = GetFieldIndex(iBand, WEIGHTED_VARIANCE); iField != -1)
889 : {
890 7 : feature.SetField(iField, stats.weighted_variance());
891 : }
892 266 : if (auto iField = GetFieldIndex(iBand, WEIGHTED_SUM); iField != -1)
893 : {
894 12 : feature.SetField(iField, stats.weighted_sum());
895 : }
896 266 : if (auto iField = GetFieldIndex(iBand, WEIGHTS); iField != -1)
897 : {
898 5 : const auto &weights = stats.weights();
899 5 : feature.SetField(iField, static_cast<int>(weights.size()),
900 : weights.data());
901 : }
902 266 : }
903 :
904 : public:
905 270 : bool ZonesAreFeature() const
906 : {
907 270 : return std::holds_alternative<OGRLayer *>(m_zones);
908 : }
909 :
910 144 : bool Process(GDALProgressFunc pfnProgress, void *pProgressData)
911 : {
912 144 : if (ZonesAreFeature())
913 : {
914 125 : if (m_options.strategy == GDALZonalStatsOptions::RASTER_SEQUENTIAL)
915 : {
916 49 : return ProcessVectorZonesByChunk(pfnProgress, pProgressData);
917 : }
918 :
919 76 : return ProcessVectorZonesByFeature(pfnProgress, pProgressData);
920 : }
921 :
922 19 : return ProcessRasterZones(pfnProgress, pProgressData);
923 : }
924 :
925 : private:
926 : static std::unique_ptr<GDALDataset>
927 89 : GetVRT(GDALDataset &src, const GDALDataset &dst, bool &resampled)
928 : {
929 89 : resampled = false;
930 :
931 89 : GDALGeoTransform srcGT, dstGT;
932 89 : if (src.GetGeoTransform(srcGT) != CE_None)
933 : {
934 0 : return nullptr;
935 : }
936 89 : if (dst.GetGeoTransform(dstGT) != CE_None)
937 : {
938 0 : return nullptr;
939 : }
940 :
941 178 : CPLStringList aosOptions;
942 89 : aosOptions.AddString("-of");
943 89 : aosOptions.AddString("VRT");
944 :
945 89 : aosOptions.AddString("-ot");
946 89 : aosOptions.AddString("Float64");
947 :
948 : // Prevent warning message about Computed -srcwin outside source raster extent.
949 : // We've already tested for this an issued a more understandable message.
950 89 : aosOptions.AddString("--no-warn-about-outside-window");
951 :
952 123 : if (srcGT != dstGT || src.GetRasterXSize() != dst.GetRasterXSize() ||
953 34 : src.GetRasterYSize() != dst.GetRasterYSize())
954 : {
955 : const double dfColOffset =
956 55 : std::fmod(std::abs(srcGT.xorig - dstGT.xorig), dstGT.xscale);
957 : const double dfRowOffset =
958 55 : std::fmod(std::abs(srcGT.yorig - dstGT.yorig), dstGT.yscale);
959 :
960 55 : OGREnvelope oDstEnv;
961 55 : dst.GetExtent(&oDstEnv);
962 :
963 55 : aosOptions.AddString("-projwin");
964 55 : aosOptions.AddString(CPLSPrintf("%.17g", oDstEnv.MinX));
965 55 : aosOptions.AddString(CPLSPrintf("%.17g", oDstEnv.MaxY));
966 55 : aosOptions.AddString(CPLSPrintf("%.17g", oDstEnv.MaxX));
967 55 : aosOptions.AddString(CPLSPrintf("%.17g", oDstEnv.MinY));
968 :
969 106 : if (srcGT.xscale != dstGT.xscale || srcGT.yscale != dstGT.yscale ||
970 161 : std::abs(dfColOffset) > 1e-4 || std::abs(dfRowOffset) > 1e-4)
971 : {
972 2 : resampled = true;
973 2 : aosOptions.AddString("-r");
974 2 : aosOptions.AddString("average");
975 : }
976 :
977 55 : aosOptions.AddString("-tr");
978 55 : aosOptions.AddString(CPLSPrintf("%.17g", dstGT.xscale));
979 55 : aosOptions.AddString(CPLSPrintf("%.17g", std::abs(dstGT.yscale)));
980 : }
981 :
982 89 : std::unique_ptr<GDALDataset> ret;
983 :
984 : GDALTranslateOptions *psOptions =
985 89 : GDALTranslateOptionsNew(aosOptions.List(), nullptr);
986 89 : ret.reset(GDALDataset::FromHandle(GDALTranslate(
987 : "", GDALDataset::ToHandle(&src), psOptions, nullptr)));
988 89 : GDALTranslateOptionsFree(psOptions);
989 :
990 89 : return ret;
991 : }
992 :
993 16 : void WarnIfZonesNotCovered(const GDALRasterBand *poZonesBand) const
994 : {
995 16 : OGREnvelope oZonesEnv;
996 16 : poZonesBand->GetDataset()->GetExtent(&oZonesEnv);
997 :
998 : {
999 16 : OGREnvelope oSrcEnv;
1000 16 : m_src.GetExtent(&oSrcEnv);
1001 :
1002 16 : if (!oZonesEnv.Intersects(oSrcEnv))
1003 : {
1004 : // TODO: Make this an error? Or keep it as a warning but short-circuit to avoid reading pixels?
1005 2 : CPLError(CE_Warning, CPLE_AppDefined,
1006 : "Source raster does not intersect zones raster");
1007 : }
1008 14 : else if (!oSrcEnv.Contains(oZonesEnv))
1009 : {
1010 : int bHasNoData;
1011 2 : m_src.GetRasterBand(m_options.bands.front())
1012 2 : ->GetNoDataValue(&bHasNoData);
1013 2 : if (bHasNoData)
1014 : {
1015 1 : CPLError(CE_Warning, CPLE_AppDefined,
1016 : "Source raster does not fully cover zones raster."
1017 : "Pixels that do not intersect the values raster "
1018 : "will be treated as having a NoData value.");
1019 : }
1020 : else
1021 : {
1022 1 : CPLError(CE_Warning, CPLE_AppDefined,
1023 : "Source raster does not fully cover zones raster. "
1024 : "Pixels that do not intersect the value raster "
1025 : "will be treated as having value of zero.");
1026 : }
1027 : }
1028 : }
1029 :
1030 16 : if (!m_weights)
1031 : {
1032 5 : return;
1033 : }
1034 :
1035 11 : OGREnvelope oWeightsEnv;
1036 11 : m_weights->GetExtent(&oWeightsEnv);
1037 :
1038 11 : if (!oZonesEnv.Intersects(oWeightsEnv))
1039 : {
1040 : // TODO: Make this an error? Or keep it as a warning but short-circuit to avoid reading pixels?
1041 0 : CPLError(CE_Warning, CPLE_AppDefined,
1042 : "Weighting raster does not intersect zones raster");
1043 : }
1044 11 : else if (!oWeightsEnv.Contains(oZonesEnv))
1045 : {
1046 : int bHasNoData;
1047 1 : m_src.GetRasterBand(m_options.bands.front())
1048 1 : ->GetNoDataValue(&bHasNoData);
1049 1 : if (bHasNoData)
1050 : {
1051 0 : CPLError(CE_Warning, CPLE_AppDefined,
1052 : "Weighting raster does not fully cover zones raster."
1053 : "Pixels that do not intersect the weighting raster "
1054 : "will be treated as having a NoData weight.");
1055 : }
1056 : else
1057 : {
1058 1 : CPLError(CE_Warning, CPLE_AppDefined,
1059 : "Weighting raster does not fully cover zones raster. "
1060 : "Pixels that do not intersect the weighting raster "
1061 : "will be treated as having a weight of zero.");
1062 : }
1063 : }
1064 : }
1065 :
1066 19 : bool ProcessRasterZones(GDALProgressFunc pfnProgress, void *pProgressData)
1067 : {
1068 19 : if (!Init())
1069 : {
1070 3 : return false;
1071 : }
1072 :
1073 16 : GDALRasterBand *poZonesBand = std::get<GDALRasterBand *>(m_zones);
1074 16 : WarnIfZonesNotCovered(poZonesBand);
1075 :
1076 16 : OGRLayer *poDstLayer = GetOutputLayer(true);
1077 16 : if (!poDstLayer)
1078 0 : return false;
1079 :
1080 : // Align the src dataset to the zones.
1081 : bool resampled;
1082 : std::unique_ptr<GDALDataset> poAlignedValuesDS =
1083 32 : GetVRT(m_src, *poZonesBand->GetDataset(), resampled);
1084 16 : if (resampled)
1085 : {
1086 0 : CPLError(CE_Warning, CPLE_AppDefined,
1087 : "Resampled source raster to match zones using average "
1088 : "resampling.");
1089 : }
1090 :
1091 : // Align the weighting dataset to the zones.
1092 16 : std::unique_ptr<GDALDataset> poAlignedWeightsDS;
1093 16 : GDALRasterBand *poWeightsBand = nullptr;
1094 16 : if (m_weights)
1095 : {
1096 : poAlignedWeightsDS =
1097 11 : GetVRT(*m_weights, *poZonesBand->GetDataset(), resampled);
1098 11 : if (!poAlignedWeightsDS)
1099 : {
1100 0 : return false;
1101 : }
1102 11 : if (resampled)
1103 : {
1104 0 : CPLError(CE_Warning, CPLE_AppDefined,
1105 : "Resampled weighting raster to match zones using "
1106 : "average resampling.");
1107 : }
1108 :
1109 : poWeightsBand =
1110 11 : poAlignedWeightsDS->GetRasterBand(m_options.weights_band);
1111 : }
1112 :
1113 32 : std::map<double, std::vector<gdal::RasterStats<double>>> stats;
1114 :
1115 32 : auto pabyZonesBuf = CreateBuffer();
1116 16 : size_t nBufSize = 0;
1117 :
1118 : const auto windowIteratorWrapper =
1119 16 : poAlignedValuesDS->GetRasterBand(1)->IterateWindows(m_maxCells);
1120 16 : const auto nIterCount = windowIteratorWrapper.count();
1121 16 : uint64_t iWindow = 0;
1122 53 : for (const auto &oWindow : windowIteratorWrapper)
1123 : {
1124 37 : const auto nWindowSize = static_cast<size_t>(oWindow.nXSize) *
1125 37 : static_cast<size_t>(oWindow.nYSize);
1126 :
1127 37 : if (nBufSize < nWindowSize)
1128 : {
1129 16 : bool bAllocSuccess = true;
1130 16 : Realloc(m_pabyValuesBuf, nWindowSize,
1131 16 : GDALGetDataTypeSizeBytes(m_workingDataType),
1132 : bAllocSuccess);
1133 16 : Realloc(pabyZonesBuf, nWindowSize,
1134 16 : GDALGetDataTypeSizeBytes(m_zonesDataType),
1135 : bAllocSuccess);
1136 16 : Realloc(m_pabyMaskBuf, nWindowSize,
1137 16 : GDALGetDataTypeSizeBytes(m_maskDataType),
1138 : bAllocSuccess);
1139 :
1140 16 : if (m_stats_options.store_xy)
1141 : {
1142 3 : Realloc(m_padfX, oWindow.nXSize,
1143 3 : GDALGetDataTypeSizeBytes(GDT_Float64),
1144 : bAllocSuccess);
1145 3 : Realloc(m_padfY, oWindow.nYSize,
1146 3 : GDALGetDataTypeSizeBytes(GDT_Float64),
1147 : bAllocSuccess);
1148 : }
1149 :
1150 16 : if (poWeightsBand)
1151 : {
1152 11 : Realloc(m_padfWeightsBuf, nWindowSize,
1153 11 : GDALGetDataTypeSizeBytes(GDT_Float64),
1154 : bAllocSuccess);
1155 11 : Realloc(m_pabyWeightsMaskBuf, nWindowSize,
1156 11 : GDALGetDataTypeSizeBytes(m_maskDataType),
1157 : bAllocSuccess);
1158 : }
1159 16 : if (!bAllocSuccess)
1160 : {
1161 0 : return false;
1162 : }
1163 :
1164 16 : nBufSize = nWindowSize;
1165 : }
1166 :
1167 37 : if (m_padfX && m_padfY)
1168 : {
1169 24 : CalculateCellCenters(oWindow, m_srcGT, m_padfX.get(),
1170 : m_padfY.get());
1171 : }
1172 :
1173 37 : if (!ReadWindow(*poZonesBand, oWindow, pabyZonesBuf.get(),
1174 : m_zonesDataType))
1175 : {
1176 0 : return false;
1177 : }
1178 :
1179 37 : if (poWeightsBand)
1180 : {
1181 32 : if (!ReadWindow(
1182 : *poWeightsBand, oWindow,
1183 32 : reinterpret_cast<GByte *>(m_padfWeightsBuf.get()),
1184 : GDT_Float64))
1185 : {
1186 0 : return false;
1187 : }
1188 32 : if (!ReadWindow(*poWeightsBand->GetMaskBand(), oWindow,
1189 : m_pabyWeightsMaskBuf.get(), GDT_UInt8))
1190 : {
1191 0 : return false;
1192 : }
1193 : }
1194 :
1195 82 : for (size_t i = 0; i < m_options.bands.size(); i++)
1196 : {
1197 45 : const int iBand = m_options.bands[i];
1198 :
1199 : GDALRasterBand *poBand =
1200 45 : poAlignedValuesDS->GetRasterBand(iBand);
1201 :
1202 45 : if (!ReadWindow(*poBand, oWindow, m_pabyValuesBuf.get(),
1203 45 : m_workingDataType))
1204 : {
1205 0 : return false;
1206 : }
1207 :
1208 45 : if (!ReadWindow(*poBand->GetMaskBand(), oWindow,
1209 45 : m_pabyMaskBuf.get(), m_maskDataType))
1210 : {
1211 0 : return false;
1212 : }
1213 :
1214 45 : size_t ipx = 0;
1215 985 : for (int k = 0; k < oWindow.nYSize; k++)
1216 : {
1217 12540 : for (int j = 0; j < oWindow.nXSize; j++)
1218 : {
1219 : // TODO use inner loop to search for a block of constant pixel values.
1220 : double zone =
1221 11600 : reinterpret_cast<double *>(pabyZonesBuf.get())[ipx];
1222 :
1223 11600 : auto &aoStats = stats[zone];
1224 11600 : aoStats.resize(m_options.bands.size(), CreateStats());
1225 :
1226 68800 : aoStats[i].process(
1227 11600 : reinterpret_cast<double *>(m_pabyValuesBuf.get()) +
1228 : ipx,
1229 11600 : m_pabyMaskBuf.get() + ipx,
1230 11600 : m_padfWeightsBuf.get()
1231 10800 : ? m_padfWeightsBuf.get() + ipx
1232 : : nullptr,
1233 11600 : m_pabyWeightsMaskBuf.get()
1234 10800 : ? m_pabyWeightsMaskBuf.get() + ipx
1235 : : nullptr,
1236 21600 : m_padfX ? m_padfX.get() + j : nullptr,
1237 21600 : m_padfY ? m_padfY.get() + k : nullptr, 1, 1);
1238 :
1239 11600 : ipx++;
1240 : }
1241 : }
1242 : }
1243 :
1244 37 : if (pfnProgress != nullptr)
1245 : {
1246 0 : ++iWindow;
1247 0 : pfnProgress(static_cast<double>(iWindow) /
1248 0 : static_cast<double>(nIterCount),
1249 : "", pProgressData);
1250 : }
1251 : }
1252 :
1253 81 : for (const auto &[dfValue, zoneStats] : stats)
1254 : {
1255 65 : OGRFeature oFeature(poDstLayer->GetLayerDefn());
1256 65 : oFeature.SetField("value", dfValue);
1257 135 : for (size_t i = 0; i < m_options.bands.size(); i++)
1258 : {
1259 70 : const auto iBand = m_options.bands[i];
1260 70 : SetStatFields(oFeature, iBand, zoneStats[i]);
1261 : }
1262 65 : if (poDstLayer->CreateFeature(&oFeature) != OGRERR_NONE)
1263 : {
1264 0 : return false;
1265 : }
1266 : }
1267 :
1268 16 : return true;
1269 : }
1270 :
1271 717 : static bool ReadWindow(GDALRasterBand &band,
1272 : const GDALRasterWindow &oWindow, GByte *pabyBuf,
1273 : GDALDataType dataType)
1274 : {
1275 1434 : return band.RasterIO(GF_Read, oWindow.nXOff, oWindow.nYOff,
1276 717 : oWindow.nXSize, oWindow.nYSize, pabyBuf,
1277 717 : oWindow.nXSize, oWindow.nYSize, dataType, 0, 0,
1278 717 : nullptr) == CE_None;
1279 : }
1280 :
1281 : #ifndef HAVE_GEOS
1282 : bool ProcessVectorZonesByChunk(GDALProgressFunc, void *)
1283 : {
1284 : CPLError(CE_Failure, CPLE_AppDefined,
1285 : "The GEOS library is required to iterate over blocks of the "
1286 : "input rasters. Processing can be performed by iterating over "
1287 : "the input features instead.");
1288 : return false;
1289 : #else
1290 49 : bool ProcessVectorZonesByChunk(GDALProgressFunc pfnProgress,
1291 : void *pProgressData)
1292 : {
1293 49 : if (!Init())
1294 : {
1295 1 : return false;
1296 : }
1297 :
1298 48 : std::unique_ptr<GDALDataset> poAlignedWeightsDS;
1299 : // Align the weighting dataset to the values.
1300 48 : if (m_weights)
1301 : {
1302 27 : bool resampled = false;
1303 27 : poAlignedWeightsDS = GetVRT(*m_weights, m_src, resampled);
1304 27 : if (!poAlignedWeightsDS)
1305 : {
1306 0 : return false;
1307 : }
1308 27 : if (resampled)
1309 : {
1310 1 : CPLError(CE_Warning, CPLE_AppDefined,
1311 : "Resampled weights to match source raster using "
1312 : "average resampling.");
1313 : }
1314 : }
1315 :
1316 48 : auto TreeDeleter = [this](GEOSSTRtree *tree)
1317 48 : { GEOSSTRtree_destroy_r(m_geosContext, tree); };
1318 :
1319 : std::unique_ptr<GEOSSTRtree, decltype(TreeDeleter)> tree(
1320 96 : GEOSSTRtree_create_r(m_geosContext, 10), TreeDeleter);
1321 :
1322 96 : std::vector<std::unique_ptr<OGRFeature>> features;
1323 96 : std::map<int, std::vector<gdal::RasterStats<double>>> statsMap;
1324 :
1325 : // Construct spatial index of all input features, storing the index
1326 : // of the feature.
1327 : {
1328 48 : OGREnvelope oGeomExtent;
1329 121 : for (auto &poFeatureIn : *std::get<OGRLayer *>(m_zones))
1330 : {
1331 74 : features.emplace_back(poFeatureIn.release());
1332 :
1333 74 : const OGRGeometry *poGeom = features.back()->GetGeometryRef();
1334 :
1335 74 : if (poGeom == nullptr || poGeom->IsEmpty())
1336 : {
1337 9 : continue;
1338 : }
1339 :
1340 65 : if (poGeom->getDimension() != 2)
1341 : {
1342 1 : CPLError(CE_Failure, CPLE_AppDefined,
1343 : "Non-polygonal geometry encountered.");
1344 1 : return false;
1345 : }
1346 :
1347 64 : poGeom->getEnvelope(&oGeomExtent);
1348 64 : GEOSGeometry *poEnv = CreateGEOSEnvelope(oGeomExtent);
1349 64 : if (poEnv == nullptr)
1350 : {
1351 0 : return false;
1352 : }
1353 :
1354 64 : GEOSSTRtree_insert_r(
1355 : m_geosContext, tree.get(), poEnv,
1356 64 : reinterpret_cast<void *>(features.size() - 1));
1357 64 : GEOSGeom_destroy_r(m_geosContext, poEnv);
1358 : }
1359 : }
1360 :
1361 97 : for (int iBand : m_options.bands)
1362 : {
1363 50 : statsMap[iBand].resize(features.size(), CreateStats());
1364 : }
1365 :
1366 94 : std::vector<void *> aiHits;
1367 92 : auto addHit = [](void *hit, void *hits)
1368 92 : { static_cast<std::vector<void *> *>(hits)->push_back(hit); };
1369 47 : size_t nBufSize = 0;
1370 :
1371 : const auto windowIteratorWrapper =
1372 47 : m_src.GetRasterBand(m_options.bands.front())
1373 47 : ->IterateWindows(m_maxCells);
1374 47 : const auto nIterCount = windowIteratorWrapper.count();
1375 47 : uint64_t iWindow = 0;
1376 157 : for (const auto &oChunkWindow : windowIteratorWrapper)
1377 : {
1378 110 : const size_t nWindowSize =
1379 110 : static_cast<size_t>(oChunkWindow.nXSize) *
1380 110 : static_cast<size_t>(oChunkWindow.nYSize);
1381 110 : const OGREnvelope oChunkExtent = ToEnvelope(oChunkWindow);
1382 :
1383 110 : aiHits.clear();
1384 :
1385 : {
1386 110 : GEOSGeometry *poEnv = CreateGEOSEnvelope(oChunkExtent);
1387 110 : if (poEnv == nullptr)
1388 : {
1389 0 : return false;
1390 : }
1391 :
1392 110 : GEOSSTRtree_query_r(m_geosContext, tree.get(), poEnv, addHit,
1393 : &aiHits);
1394 110 : GEOSGeom_destroy_r(m_geosContext, poEnv);
1395 : }
1396 :
1397 110 : if (!aiHits.empty())
1398 : {
1399 74 : if (nBufSize < nWindowSize)
1400 : {
1401 38 : bool bAllocSuccess = true;
1402 38 : Realloc(m_pabyValuesBuf, nWindowSize,
1403 38 : GDALGetDataTypeSizeBytes(m_workingDataType),
1404 : bAllocSuccess);
1405 38 : Realloc(m_pabyCoverageBuf, nWindowSize,
1406 38 : GDALGetDataTypeSizeBytes(m_coverageDataType),
1407 : bAllocSuccess);
1408 38 : Realloc(m_pabyMaskBuf, nWindowSize,
1409 38 : GDALGetDataTypeSizeBytes(m_maskDataType),
1410 : bAllocSuccess);
1411 38 : if (m_stats_options.store_xy)
1412 : {
1413 9 : Realloc(m_padfX, oChunkWindow.nXSize,
1414 9 : GDALGetDataTypeSizeBytes(GDT_Float64),
1415 : bAllocSuccess);
1416 9 : Realloc(m_padfY, oChunkWindow.nYSize,
1417 9 : GDALGetDataTypeSizeBytes(GDT_Float64),
1418 : bAllocSuccess);
1419 : }
1420 38 : if (m_weights != nullptr)
1421 : {
1422 27 : Realloc(m_padfWeightsBuf, nWindowSize,
1423 27 : GDALGetDataTypeSizeBytes(GDT_Float64),
1424 : bAllocSuccess);
1425 27 : Realloc(m_pabyWeightsMaskBuf, nWindowSize,
1426 27 : GDALGetDataTypeSizeBytes(m_maskDataType),
1427 : bAllocSuccess);
1428 : }
1429 38 : if (!bAllocSuccess)
1430 : {
1431 0 : return false;
1432 : }
1433 38 : nBufSize = nWindowSize;
1434 : }
1435 :
1436 74 : if (m_padfX && m_padfY)
1437 : {
1438 21 : CalculateCellCenters(oChunkWindow, m_srcGT, m_padfX.get(),
1439 : m_padfY.get());
1440 : }
1441 :
1442 74 : if (m_weights != nullptr)
1443 : {
1444 : GDALRasterBand *poWeightsBand =
1445 27 : poAlignedWeightsDS->GetRasterBand(
1446 : m_options.weights_band);
1447 :
1448 27 : if (!ReadWindow(
1449 : *poWeightsBand, oChunkWindow,
1450 27 : reinterpret_cast<GByte *>(m_padfWeightsBuf.get()),
1451 : GDT_Float64))
1452 : {
1453 0 : return false;
1454 : }
1455 27 : if (!ReadWindow(*poWeightsBand->GetMaskBand(), oChunkWindow,
1456 : m_pabyWeightsMaskBuf.get(), GDT_UInt8))
1457 : {
1458 0 : return false;
1459 : }
1460 : }
1461 :
1462 163 : for (int iBand : m_options.bands)
1463 : {
1464 :
1465 89 : GDALRasterBand *poBand = m_src.GetRasterBand(iBand);
1466 :
1467 178 : if (!(ReadWindow(*poBand, oChunkWindow,
1468 : m_pabyValuesBuf.get(),
1469 89 : m_workingDataType) &&
1470 89 : ReadWindow(*poBand->GetMaskBand(), oChunkWindow,
1471 89 : m_pabyMaskBuf.get(), m_maskDataType)))
1472 : {
1473 0 : return false;
1474 : }
1475 :
1476 : GDALRasterWindow oGeomWindow;
1477 89 : OGREnvelope oGeomExtent;
1478 199 : for (const void *hit : aiHits)
1479 : {
1480 110 : const size_t iHit = reinterpret_cast<size_t>(hit);
1481 110 : const auto poGeom = features[iHit]->GetGeometryRef();
1482 :
1483 : // Trim the chunk window to the portion that intersects
1484 : // the geometry being processed.
1485 110 : poGeom->getEnvelope(&oGeomExtent);
1486 110 : oGeomExtent.Intersect(oChunkExtent);
1487 110 : if (!m_srcInvGT.Apply(oGeomExtent, oGeomWindow))
1488 : {
1489 0 : return false;
1490 : }
1491 110 : oGeomWindow.nXOff =
1492 110 : std::max(oGeomWindow.nXOff, oChunkWindow.nXOff);
1493 110 : oGeomWindow.nYOff =
1494 110 : std::max(oGeomWindow.nYOff, oChunkWindow.nYOff);
1495 110 : oGeomWindow.nXSize =
1496 110 : std::min(oGeomWindow.nXSize,
1497 220 : oChunkWindow.nXOff + oChunkWindow.nXSize -
1498 110 : oGeomWindow.nXOff);
1499 110 : oGeomWindow.nYSize =
1500 110 : std::min(oGeomWindow.nYSize,
1501 220 : oChunkWindow.nYOff + oChunkWindow.nYSize -
1502 110 : oGeomWindow.nYOff);
1503 110 : if (oGeomWindow.nXSize <= 0 || oGeomWindow.nYSize <= 0)
1504 0 : continue;
1505 : const OGREnvelope oTrimmedEnvelope =
1506 110 : ToEnvelope(oGeomWindow);
1507 :
1508 110 : if (!CalculateCoverage(
1509 : poGeom, oTrimmedEnvelope, oGeomWindow.nXSize,
1510 : oGeomWindow.nYSize, m_pabyCoverageBuf.get()))
1511 : {
1512 0 : return false;
1513 : }
1514 :
1515 : // Because the window used for polygon coverage is not the
1516 : // same as the window used for raster values, iterate
1517 : // over partial scanlines on the raster window.
1518 110 : const auto nCoverageXOff =
1519 110 : oGeomWindow.nXOff - oChunkWindow.nXOff;
1520 110 : const auto nCoverageYOff =
1521 110 : oGeomWindow.nYOff - oChunkWindow.nYOff;
1522 859 : for (int iRow = 0; iRow < oGeomWindow.nYSize; iRow++)
1523 : {
1524 749 : const auto nFirstPx =
1525 749 : (nCoverageYOff + iRow) * oChunkWindow.nXSize +
1526 : nCoverageXOff;
1527 3745 : UpdateStats(
1528 749 : statsMap[iBand][iHit],
1529 1498 : m_pabyValuesBuf.get() +
1530 749 : nFirstPx * GDALGetDataTypeSizeBytes(
1531 749 : m_workingDataType),
1532 749 : m_pabyMaskBuf.get() +
1533 749 : nFirstPx * GDALGetDataTypeSizeBytes(
1534 749 : m_maskDataType),
1535 : m_padfWeightsBuf
1536 80 : ? m_padfWeightsBuf.get() + nFirstPx
1537 749 : : nullptr,
1538 : m_pabyWeightsMaskBuf
1539 80 : ? m_pabyWeightsMaskBuf.get() +
1540 80 : nFirstPx * GDALGetDataTypeSizeBytes(
1541 80 : m_maskDataType)
1542 749 : : nullptr,
1543 749 : m_pabyCoverageBuf.get() +
1544 1498 : iRow * oGeomWindow.nXSize *
1545 749 : GDALGetDataTypeSizeBytes(
1546 749 : m_coverageDataType),
1547 168 : m_padfX ? m_padfX.get() + nCoverageXOff
1548 749 : : nullptr,
1549 168 : m_padfY ? m_padfY.get() + nCoverageYOff + iRow
1550 749 : : nullptr,
1551 749 : oGeomWindow.nXSize, 1);
1552 : }
1553 : }
1554 : }
1555 : }
1556 :
1557 110 : if (pfnProgress != nullptr)
1558 : {
1559 0 : ++iWindow;
1560 0 : pfnProgress(static_cast<double>(iWindow) /
1561 0 : static_cast<double>(nIterCount),
1562 : "", pProgressData);
1563 : }
1564 : }
1565 :
1566 47 : OGRLayer *poDstLayer = GetOutputLayer(false);
1567 47 : if (!poDstLayer)
1568 0 : return false;
1569 :
1570 120 : for (size_t iFeature = 0; iFeature < features.size(); iFeature++)
1571 : {
1572 : auto poDstFeature =
1573 73 : std::make_unique<OGRFeature>(poDstLayer->GetLayerDefn());
1574 73 : poDstFeature->SetFrom(features[iFeature].get());
1575 152 : for (int iBand : m_options.bands)
1576 : {
1577 79 : SetStatFields(*poDstFeature, iBand, statsMap[iBand][iFeature]);
1578 : }
1579 73 : if (poDstLayer->CreateFeature(poDstFeature.get()) != OGRERR_NONE)
1580 : {
1581 0 : return false;
1582 : }
1583 : }
1584 :
1585 47 : return true;
1586 : #endif
1587 : }
1588 :
1589 76 : bool ProcessVectorZonesByFeature(GDALProgressFunc pfnProgress,
1590 : void *pProgressData)
1591 : {
1592 76 : if (!Init())
1593 : {
1594 17 : return false;
1595 : }
1596 :
1597 59 : OGREnvelope oGeomExtent;
1598 : GDALRasterWindow oWindow;
1599 :
1600 59 : std::unique_ptr<GDALDataset> poAlignedWeightsDS;
1601 : // Align the weighting dataset to the values.
1602 59 : if (m_weights)
1603 : {
1604 35 : bool resampled = false;
1605 35 : poAlignedWeightsDS = GetVRT(*m_weights, m_src, resampled);
1606 35 : if (!poAlignedWeightsDS)
1607 : {
1608 0 : return false;
1609 : }
1610 35 : if (resampled)
1611 : {
1612 1 : CPLError(CE_Warning, CPLE_AppDefined,
1613 : "Resampled weights to match source raster using "
1614 : "average resampling.");
1615 : }
1616 : }
1617 :
1618 59 : size_t nBufSize = 0;
1619 :
1620 59 : OGRLayer *poSrcLayer = std::get<OGRLayer *>(m_zones);
1621 59 : OGRLayer *poDstLayer = GetOutputLayer(false);
1622 59 : if (!poDstLayer)
1623 0 : return false;
1624 59 : size_t i = 0;
1625 59 : auto nFeatures = poSrcLayer->GetFeatureCount();
1626 : GDALRasterWindow oRasterWindow;
1627 59 : oRasterWindow.nXOff = 0;
1628 59 : oRasterWindow.nYOff = 0;
1629 59 : oRasterWindow.nXSize = m_src.GetRasterXSize();
1630 59 : oRasterWindow.nYSize = m_src.GetRasterYSize();
1631 59 : const OGREnvelope oRasterExtent = ToEnvelope(oRasterWindow);
1632 :
1633 170 : for (const auto &poFeature : *poSrcLayer)
1634 : {
1635 112 : const auto *poGeom = poFeature->GetGeometryRef();
1636 :
1637 112 : oWindow.nXSize = 0;
1638 112 : oWindow.nYSize = 0;
1639 112 : if (poGeom == nullptr || poGeom->IsEmpty())
1640 : {
1641 : // do nothing
1642 : }
1643 103 : else if (poGeom->getDimension() != 2)
1644 : {
1645 1 : CPLError(CE_Failure, CPLE_AppDefined,
1646 : "Non-polygonal geometry encountered.");
1647 1 : return false;
1648 : }
1649 : else
1650 : {
1651 102 : poGeom->getEnvelope(&oGeomExtent);
1652 102 : if (oGeomExtent.Intersects(oRasterExtent))
1653 : {
1654 94 : oGeomExtent.Intersect(oRasterExtent);
1655 94 : if (!m_srcInvGT.Apply(oGeomExtent, oWindow))
1656 : {
1657 0 : return false;
1658 : }
1659 94 : oWindow.nXOff =
1660 94 : std::max(oWindow.nXOff, oRasterWindow.nXOff);
1661 94 : oWindow.nYOff =
1662 94 : std::max(oWindow.nYOff, oRasterWindow.nYOff);
1663 94 : oWindow.nXSize =
1664 188 : std::min(oWindow.nXSize, oRasterWindow.nXOff +
1665 188 : oRasterWindow.nXSize -
1666 94 : oWindow.nXOff);
1667 94 : oWindow.nYSize =
1668 188 : std::min(oWindow.nYSize, oRasterWindow.nYOff +
1669 94 : oRasterWindow.nYSize -
1670 94 : oWindow.nYOff);
1671 : }
1672 : }
1673 :
1674 : std::unique_ptr<OGRFeature> poDstFeature(
1675 111 : OGRFeature::CreateFeature(poDstLayer->GetLayerDefn()));
1676 111 : poDstFeature->SetFrom(poFeature.get());
1677 :
1678 111 : if (oWindow.nXSize == 0 || oWindow.nYSize == 0)
1679 : {
1680 34 : const gdal::RasterStats<double> empty(CreateStats());
1681 34 : for (int iBand : m_options.bands)
1682 : {
1683 17 : SetStatFields(*poDstFeature, iBand, empty);
1684 17 : }
1685 : }
1686 : else
1687 : {
1688 : // Calculate how many rows of raster data we can read in at
1689 : // a time while remaining within maxCells.
1690 94 : const int nRowsPerChunk = std::min(
1691 : oWindow.nYSize,
1692 188 : std::max(1, static_cast<int>(
1693 94 : m_maxCells /
1694 94 : static_cast<size_t>(oWindow.nXSize))));
1695 :
1696 94 : const size_t nWindowSize = static_cast<size_t>(oWindow.nXSize) *
1697 94 : static_cast<size_t>(nRowsPerChunk);
1698 :
1699 94 : if (nBufSize < nWindowSize)
1700 : {
1701 61 : bool bAllocSuccess = true;
1702 61 : Realloc(m_pabyValuesBuf, nWindowSize,
1703 61 : GDALGetDataTypeSizeBytes(m_workingDataType),
1704 : bAllocSuccess);
1705 61 : Realloc(m_pabyCoverageBuf, nWindowSize,
1706 61 : GDALGetDataTypeSizeBytes(m_coverageDataType),
1707 : bAllocSuccess);
1708 61 : Realloc(m_pabyMaskBuf, nWindowSize,
1709 61 : GDALGetDataTypeSizeBytes(m_maskDataType),
1710 : bAllocSuccess);
1711 :
1712 61 : if (m_stats_options.store_xy)
1713 : {
1714 9 : Realloc(m_padfX, oWindow.nXSize,
1715 9 : GDALGetDataTypeSizeBytes(GDT_Float64),
1716 : bAllocSuccess);
1717 9 : Realloc(m_padfY, oWindow.nYSize,
1718 9 : GDALGetDataTypeSizeBytes(GDT_Float64),
1719 : bAllocSuccess);
1720 : }
1721 :
1722 61 : if (m_weights != nullptr)
1723 : {
1724 35 : Realloc(m_padfWeightsBuf, nWindowSize,
1725 35 : GDALGetDataTypeSizeBytes(GDT_Float64),
1726 : bAllocSuccess);
1727 35 : Realloc(m_pabyWeightsMaskBuf, nWindowSize,
1728 35 : GDALGetDataTypeSizeBytes(m_maskDataType),
1729 : bAllocSuccess);
1730 : }
1731 61 : if (!bAllocSuccess)
1732 : {
1733 0 : return false;
1734 : }
1735 :
1736 61 : nBufSize = nWindowSize;
1737 : }
1738 :
1739 94 : if (m_padfX && m_padfY)
1740 : {
1741 12 : CalculateCellCenters(oWindow, m_srcGT, m_padfX.get(),
1742 : m_padfY.get());
1743 : }
1744 :
1745 94 : std::vector<gdal::RasterStats<double>> aoStats;
1746 94 : aoStats.resize(m_options.bands.size(), CreateStats());
1747 :
1748 94 : for (int nYOff = oWindow.nYOff;
1749 197 : nYOff < oWindow.nYOff + oWindow.nYSize;
1750 103 : nYOff += nRowsPerChunk)
1751 : {
1752 : GDALRasterWindow oSubWindow;
1753 103 : oSubWindow.nXOff = oWindow.nXOff;
1754 103 : oSubWindow.nXSize = oWindow.nXSize;
1755 103 : oSubWindow.nYOff = nYOff;
1756 103 : oSubWindow.nYSize = std::min(
1757 103 : nRowsPerChunk, oWindow.nYOff + oWindow.nYSize - nYOff);
1758 :
1759 103 : const auto nCoverageXOff = oSubWindow.nXOff - oWindow.nXOff;
1760 103 : const auto nCoverageYOff = oSubWindow.nYOff - oWindow.nYOff;
1761 :
1762 : const OGREnvelope oSnappedGeomExtent =
1763 103 : ToEnvelope(oSubWindow);
1764 :
1765 103 : if (!CalculateCoverage(poGeom, oSnappedGeomExtent,
1766 : oSubWindow.nXSize, oSubWindow.nYSize,
1767 : m_pabyCoverageBuf.get()))
1768 : {
1769 0 : return false;
1770 : }
1771 :
1772 103 : if (m_weights != nullptr)
1773 : {
1774 : GDALRasterBand *poWeightsBand =
1775 35 : poAlignedWeightsDS->GetRasterBand(
1776 : m_options.weights_band);
1777 :
1778 35 : if (!ReadWindow(*poWeightsBand, oSubWindow,
1779 : reinterpret_cast<GByte *>(
1780 35 : m_padfWeightsBuf.get()),
1781 : GDT_Float64))
1782 : {
1783 0 : return false;
1784 : }
1785 35 : if (!ReadWindow(*poWeightsBand->GetMaskBand(),
1786 : oSubWindow, m_pabyWeightsMaskBuf.get(),
1787 : GDT_UInt8))
1788 : {
1789 0 : return false;
1790 : }
1791 : }
1792 :
1793 215 : for (size_t iBandInd = 0; iBandInd < m_options.bands.size();
1794 : iBandInd++)
1795 : {
1796 : GDALRasterBand *poBand =
1797 112 : m_src.GetRasterBand(m_options.bands[iBandInd]);
1798 :
1799 112 : if (!ReadWindow(*poBand, oSubWindow,
1800 : m_pabyValuesBuf.get(),
1801 112 : m_workingDataType))
1802 : {
1803 0 : return false;
1804 : }
1805 112 : if (!ReadWindow(*poBand->GetMaskBand(), oSubWindow,
1806 112 : m_pabyMaskBuf.get(), m_maskDataType))
1807 : {
1808 0 : return false;
1809 : }
1810 :
1811 336 : UpdateStats(
1812 112 : aoStats[iBandInd], m_pabyValuesBuf.get(),
1813 112 : m_pabyMaskBuf.get(), m_padfWeightsBuf.get(),
1814 112 : m_pabyWeightsMaskBuf.get(), m_pabyCoverageBuf.get(),
1815 127 : m_padfX ? m_padfX.get() + nCoverageXOff : nullptr,
1816 15 : m_padfY ? m_padfY.get() + nCoverageYOff : nullptr,
1817 112 : oSubWindow.nXSize, oSubWindow.nYSize);
1818 : }
1819 : }
1820 :
1821 194 : for (size_t iBandInd = 0; iBandInd < m_options.bands.size();
1822 : iBandInd++)
1823 : {
1824 100 : SetStatFields(*poDstFeature, m_options.bands[iBandInd],
1825 100 : aoStats[iBandInd]);
1826 : }
1827 : }
1828 :
1829 111 : if (poDstLayer->CreateFeature(poDstFeature.get()) != OGRERR_NONE)
1830 : {
1831 0 : return false;
1832 : }
1833 :
1834 111 : if (pfnProgress)
1835 : {
1836 0 : pfnProgress(static_cast<double>(i + 1) /
1837 0 : static_cast<double>(nFeatures),
1838 : "", pProgressData);
1839 : }
1840 111 : i++;
1841 : }
1842 :
1843 58 : return true;
1844 : }
1845 :
1846 861 : void UpdateStats(gdal::RasterStats<double> &stats, const GByte *pabyValues,
1847 : const GByte *pabyMask, const double *padfWeights,
1848 : const GByte *pabyWeightsMask, const GByte *pabyCoverage,
1849 : const double *pdfX, const double *pdfY, size_t nX,
1850 : size_t nY) const
1851 : {
1852 861 : if (m_coverageDataType == GDT_Float32)
1853 : {
1854 312 : stats.process(reinterpret_cast<const double *>(pabyValues),
1855 : pabyMask, padfWeights, pabyWeightsMask,
1856 : reinterpret_cast<const float *>(pabyCoverage), pdfX,
1857 : pdfY, nX, nY);
1858 : }
1859 : else
1860 : {
1861 549 : stats.process(reinterpret_cast<const double *>(pabyValues),
1862 : pabyMask, padfWeights, pabyWeightsMask, pabyCoverage,
1863 : pdfX, pdfY, nX, nY);
1864 : }
1865 861 : }
1866 :
1867 213 : bool CalculateCoverage(const OGRGeometry *poGeom,
1868 : const OGREnvelope &oSnappedGeomExtent, int nXSize,
1869 : int nYSize, GByte *pabyCoverageBuf) const
1870 : {
1871 : #if GEOS_GRID_INTERSECTION_AVAILABLE
1872 213 : if (m_options.pixels == GDALZonalStatsOptions::FRACTIONAL)
1873 : {
1874 83 : std::memset(pabyCoverageBuf, 0,
1875 83 : static_cast<size_t>(nXSize) * nYSize *
1876 83 : GDALGetDataTypeSizeBytes(GDT_Float32));
1877 : GEOSGeometry *poGeosGeom =
1878 83 : poGeom->exportToGEOS(m_geosContext, true);
1879 83 : if (!poGeosGeom)
1880 : {
1881 0 : CPLError(CE_Failure, CPLE_AppDefined,
1882 : "Failed to convert geometry to GEOS.");
1883 0 : return false;
1884 : }
1885 :
1886 166 : const bool bRet = GEOSGridIntersectionFractions_r(
1887 83 : m_geosContext, poGeosGeom, oSnappedGeomExtent.MinX,
1888 83 : oSnappedGeomExtent.MinY, oSnappedGeomExtent.MaxX,
1889 83 : oSnappedGeomExtent.MaxY, nXSize, nYSize,
1890 83 : reinterpret_cast<float *>(pabyCoverageBuf));
1891 83 : if (!bRet)
1892 : {
1893 0 : CPLError(CE_Failure, CPLE_AppDefined,
1894 : "Failed to calculate pixel intersection fractions.");
1895 : }
1896 83 : GEOSGeom_destroy_r(m_geosContext, poGeosGeom);
1897 :
1898 83 : return bRet;
1899 : }
1900 : else
1901 : #endif
1902 : {
1903 130 : GDALGeoTransform oCoverageGT;
1904 130 : oCoverageGT.xorig = oSnappedGeomExtent.MinX;
1905 130 : oCoverageGT.xscale = m_srcGT.xscale;
1906 130 : oCoverageGT.xrot = 0;
1907 :
1908 130 : oCoverageGT.yorig = m_srcGT.yscale < 0 ? oSnappedGeomExtent.MaxY
1909 : : oSnappedGeomExtent.MinY;
1910 130 : oCoverageGT.yscale = m_srcGT.yscale;
1911 130 : oCoverageGT.yrot = 0;
1912 :
1913 : // Create a memory dataset that wraps the coverage buffer so that
1914 : // we can invoke GDALRasterize
1915 : std::unique_ptr<MEMDataset> poMemDS(MEMDataset::Create(
1916 260 : "", nXSize, nYSize, 0, m_coverageDataType, nullptr));
1917 130 : poMemDS->SetGeoTransform(oCoverageGT);
1918 130 : constexpr double dfBurnValue = 255.0;
1919 130 : constexpr int nBand = 1;
1920 :
1921 : MEMRasterBand *poCoverageBand =
1922 130 : new MEMRasterBand(poMemDS.get(), 1, pabyCoverageBuf,
1923 130 : m_coverageDataType, 0, 0, false, nullptr);
1924 130 : poMemDS->AddMEMBand(poCoverageBand);
1925 130 : poCoverageBand->Fill(0);
1926 :
1927 130 : CPLStringList aosOptions;
1928 130 : if (m_options.pixels == GDALZonalStatsOptions::ALL_TOUCHED)
1929 : {
1930 35 : aosOptions.AddString("ALL_TOUCHED=1");
1931 : }
1932 :
1933 : OGRGeometryH hGeom =
1934 130 : OGRGeometry::ToHandle(const_cast<OGRGeometry *>(poGeom));
1935 :
1936 130 : const auto eErr = GDALRasterizeGeometries(
1937 130 : GDALDataset::ToHandle(poMemDS.get()), 1, &nBand, 1, &hGeom,
1938 130 : nullptr, nullptr, &dfBurnValue, aosOptions.List(), nullptr,
1939 : nullptr);
1940 :
1941 130 : return eErr == CE_None;
1942 : }
1943 : }
1944 :
1945 : #ifdef HAVE_GEOS
1946 174 : GEOSGeometry *CreateGEOSEnvelope(const OGREnvelope &oEnv) const
1947 : {
1948 174 : GEOSCoordSequence *seq = GEOSCoordSeq_create_r(m_geosContext, 2, 2);
1949 174 : if (seq == nullptr)
1950 : {
1951 0 : return nullptr;
1952 : }
1953 174 : GEOSCoordSeq_setXY_r(m_geosContext, seq, 0, oEnv.MinX, oEnv.MinY);
1954 174 : GEOSCoordSeq_setXY_r(m_geosContext, seq, 1, oEnv.MaxX, oEnv.MaxY);
1955 174 : return GEOSGeom_createLineString_r(m_geosContext, seq);
1956 : }
1957 : #endif
1958 :
1959 : CPL_DISALLOW_COPY_ASSIGN(GDALZonalStatsImpl)
1960 :
1961 : GDALDataset &m_src;
1962 : GDALDataset *m_weights;
1963 : GDALDataset &m_dst;
1964 : const BandOrLayer m_zones;
1965 :
1966 : const GDALDataType m_coverageDataType;
1967 : const GDALDataType m_workingDataType = GDT_Float64;
1968 : const GDALDataType m_maskDataType = GDT_UInt8;
1969 : static constexpr GDALDataType m_zonesDataType = GDT_Float64;
1970 :
1971 : GDALGeoTransform m_srcGT{};
1972 : GDALGeoTransform m_srcInvGT{};
1973 :
1974 : GDALZonalStatsOptions m_options{};
1975 : gdal::RasterStatsOptions m_stats_options{};
1976 :
1977 : size_t m_maxCells{0};
1978 :
1979 : static constexpr auto NUM_STATS = Stat::INVALID + 1;
1980 : std::map<int, std::array<int, NUM_STATS>> m_statFields{};
1981 :
1982 : std::unique_ptr<GByte, VSIFreeReleaser> m_pabyCoverageBuf{};
1983 : std::unique_ptr<GByte, VSIFreeReleaser> m_pabyMaskBuf{};
1984 : std::unique_ptr<GByte, VSIFreeReleaser> m_pabyValuesBuf{};
1985 : std::unique_ptr<double, VSIFreeReleaser> m_padfWeightsBuf{};
1986 : std::unique_ptr<GByte, VSIFreeReleaser> m_pabyWeightsMaskBuf{};
1987 : std::unique_ptr<double, VSIFreeReleaser> m_padfX{};
1988 : std::unique_ptr<double, VSIFreeReleaser> m_padfY{};
1989 :
1990 : #ifdef HAVE_GEOS
1991 : GEOSContextHandle_t m_geosContext{nullptr};
1992 : #endif
1993 : };
1994 :
1995 147 : static CPLErr GDALZonalStats(GDALDataset &srcDataset, GDALDataset *poWeights,
1996 : GDALDataset &zonesDataset, GDALDataset &dstDataset,
1997 : const GDALZonalStatsOptions &options,
1998 : GDALProgressFunc pfnProgress, void *pProgressData)
1999 : {
2000 147 : int nZonesBand = options.zones_band;
2001 294 : std::string osZonesLayer = options.zones_layer;
2002 :
2003 147 : if (nZonesBand < 1 && osZonesLayer.empty())
2004 : {
2005 145 : if (zonesDataset.GetRasterCount() + zonesDataset.GetLayerCount() > 1)
2006 : {
2007 0 : CPLError(CE_Failure, CPLE_AppDefined,
2008 : "Zones dataset has more than one band or layer. Use "
2009 : "the --zone-band or --zone-layer argument to specify "
2010 : "which should be used.");
2011 0 : return CE_Failure;
2012 : }
2013 145 : if (zonesDataset.GetRasterCount() > 0)
2014 : {
2015 19 : nZonesBand = 1;
2016 : }
2017 126 : else if (zonesDataset.GetLayerCount() > 0)
2018 : {
2019 125 : osZonesLayer = zonesDataset.GetLayer(0)->GetName();
2020 : }
2021 : else
2022 : {
2023 1 : CPLError(CE_Failure, CPLE_AppDefined,
2024 : "Zones dataset has no band or layer.");
2025 1 : return CE_Failure;
2026 : }
2027 : }
2028 :
2029 146 : GDALZonalStatsImpl::BandOrLayer poZones;
2030 :
2031 146 : if (nZonesBand > 0)
2032 : {
2033 20 : if (nZonesBand > zonesDataset.GetRasterCount())
2034 : {
2035 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid zones band: %d",
2036 : nZonesBand);
2037 1 : return CE_Failure;
2038 : }
2039 19 : GDALRasterBand *poZonesBand = zonesDataset.GetRasterBand(nZonesBand);
2040 19 : if (poZonesBand == nullptr)
2041 : {
2042 0 : CPLError(CE_Failure, CPLE_AppDefined,
2043 : "Specified zones band %d not found", nZonesBand);
2044 0 : return CE_Failure;
2045 : }
2046 19 : poZones = poZonesBand;
2047 : }
2048 : else
2049 : {
2050 : OGRLayer *poZonesLayer =
2051 126 : zonesDataset.GetLayerByName(osZonesLayer.c_str());
2052 126 : if (poZonesLayer == nullptr)
2053 : {
2054 1 : CPLError(CE_Failure, CPLE_AppDefined,
2055 : "Specified zones layer '%s' not found",
2056 : options.zones_layer.c_str());
2057 1 : return CE_Failure;
2058 : }
2059 125 : poZones = poZonesLayer;
2060 : }
2061 :
2062 144 : GDALZonalStatsImpl alg(srcDataset, dstDataset, poWeights, poZones, options);
2063 144 : return alg.Process(pfnProgress, pProgressData) ? CE_None : CE_Failure;
2064 : }
2065 :
2066 : /** Compute statistics of raster values within defined zones
2067 : *
2068 : * @param hSrcDS raster dataset containing values to be summarized
2069 : * @param hWeightsDS optional raster dataset containing weights
2070 : * @param hZonesDS raster or vector dataset containing zones across which values will be summarized
2071 : * @param hOutDS dataset to which output layer will be written
2072 : * @param papszOptions list of options
2073 : * BANDS: a comma-separated list of band indices to be processed from the
2074 : * source dataset. If not present, all bands will be processed.
2075 : * INCLUDE_FIELDS: a comma-separated list of field names from the zones
2076 : * dataset to be included in output features.
2077 : * PIXEL_INTERSECTION: controls which pixels are included in calculations:
2078 : * - DEFAULT: use default options to GDALRasterize
2079 : * - ALL_TOUCHED: use ALL_TOUCHED option of GDALRasterize
2080 : * - FRACTIONAL: calculate fraction of each pixel that is covered
2081 : * by the zone. Requires the GEOS library, version >= 3.14.
2082 : * RASTER_CHUNK_SIZE_BYTES: sets a maximum amount of raster data to read
2083 : * into memory at a single time (from a single source)
2084 : * STATS: comma-separated list of stats. The following stats are supported:
2085 : * - center_x
2086 : * - center_y
2087 : * - count
2088 : * - coverage
2089 : * - frac
2090 : * - max
2091 : * - max_center_x
2092 : * - max_center_y
2093 : * - mean
2094 : * - min
2095 : * - min_center_x
2096 : * - min_center_y
2097 : * - minority
2098 : * - mode
2099 : * - stdev
2100 : * - sum
2101 : * - unique
2102 : * - values
2103 : * - variance
2104 : * - weighted_frac
2105 : * - mean
2106 : * - weighted_sum
2107 : * - weighted_stdev
2108 : * - weighted_variance
2109 : * - weights
2110 : * STRATEGY: determine how to perform processing with vector zones:
2111 : * - FEATURE_SEQUENTIAL: iterate over zones, finding raster pixels
2112 : * that intersect with each, calculating stats, and writing output
2113 : * to hOutDS.
2114 : * - RASTER_SEQUENTIAL: iterate over chunks of the raster, finding
2115 : * zones that intersect with each chunk and updating stats.
2116 : * Features are written to hOutDS after all processing has been
2117 : * completed.
2118 : * WEIGHTS_BAND: the band to read from WeightsDS
2119 : * ZONES_BAND: the band to read from hZonesDS, if hZonesDS is a raster
2120 : * ZONES_LAYER: the layer to read from hZonesDS, if hZonesDS is a vector
2121 : * LCO_{key}: layer creation option {key}
2122 : *
2123 : * @param pfnProgress optional progress reporting callback
2124 : * @param pProgressArg optional data for progress callback
2125 : * @return CE_Failure if an error occurred, CE_None otherwise
2126 : */
2127 147 : CPLErr GDALZonalStats(GDALDatasetH hSrcDS, GDALDatasetH hWeightsDS,
2128 : GDALDatasetH hZonesDS, GDALDatasetH hOutDS,
2129 : CSLConstList papszOptions, GDALProgressFunc pfnProgress,
2130 : void *pProgressArg)
2131 : {
2132 147 : VALIDATE_POINTER1(hSrcDS, __func__, CE_Failure);
2133 147 : VALIDATE_POINTER1(hZonesDS, __func__, CE_Failure);
2134 147 : VALIDATE_POINTER1(hOutDS, __func__, CE_Failure);
2135 :
2136 294 : GDALZonalStatsOptions sOptions;
2137 147 : if (papszOptions)
2138 : {
2139 147 : if (auto eErr = sOptions.Init(papszOptions); eErr != CE_None)
2140 : {
2141 0 : return eErr;
2142 : }
2143 : }
2144 :
2145 294 : return GDALZonalStats(
2146 147 : *GDALDataset::FromHandle(hSrcDS), GDALDataset::FromHandle(hWeightsDS),
2147 147 : *GDALDataset::FromHandle(hZonesDS), *GDALDataset::FromHandle(hOutDS),
2148 147 : sOptions, pfnProgress, pProgressArg);
2149 : }
|
