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