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