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