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