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