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