Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Vector rasterization.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "gdal_alg.h"
16 : #include "gdal_alg_priv.h"
17 :
18 : #include <climits>
19 : #include <cstddef>
20 : #include <cstdlib>
21 : #include <cstring>
22 : #include <cfloat>
23 : #include <limits>
24 : #include <vector>
25 : #include <algorithm>
26 :
27 : #include "cpl_conv.h"
28 : #include "cpl_error.h"
29 : #include "cpl_progress.h"
30 : #include "cpl_string.h"
31 : #include "cpl_vsi.h"
32 : #include "gdal.h"
33 : #include "gdal_priv.h"
34 : #include "gdal_priv_templates.hpp"
35 : #include "ogr_api.h"
36 : #include "ogr_core.h"
37 : #include "ogr_feature.h"
38 : #include "ogr_geometry.h"
39 : #include "ogr_spatialref.h"
40 : #include "ogrsf_frmts.h"
41 :
42 0 : template <typename T> static inline T SaturatedAddSigned(T a, T b)
43 : {
44 0 : if (a > 0 && b > 0 && a > std::numeric_limits<T>::max() - b)
45 : {
46 0 : return std::numeric_limits<T>::max();
47 : }
48 0 : else if (a < 0 && b < 0 && a < std::numeric_limits<T>::min() - b)
49 : {
50 0 : return std::numeric_limits<T>::min();
51 : }
52 : else
53 : {
54 0 : return a + b;
55 : }
56 : }
57 :
58 : /************************************************************************/
59 : /* MakeKey() */
60 : /************************************************************************/
61 :
62 818 : inline uint64_t MakeKey(int y, int x)
63 : {
64 818 : return (static_cast<uint64_t>(y) << 32) | static_cast<uint64_t>(x);
65 : }
66 :
67 : /************************************************************************/
68 : /* gvBurnScanlineBasic() */
69 : /************************************************************************/
70 : template <typename T>
71 15686 : static inline void gvBurnScanlineBasic(GDALRasterizeInfo *psInfo, int nY,
72 : int nXStart, int nXEnd, double dfVariant)
73 :
74 : {
75 32362 : for (int iBand = 0; iBand < psInfo->nBands; iBand++)
76 : {
77 16676 : const double burnValue =
78 16676 : (psInfo->burnValues.double_values[iBand] +
79 16676 : ((psInfo->eBurnValueSource == GBV_UserBurnValue) ? 0 : dfVariant));
80 :
81 16676 : unsigned char *pabyInsert =
82 16676 : psInfo->pabyChunkBuf + iBand * psInfo->nBandSpace +
83 16676 : nY * psInfo->nLineSpace + nXStart * psInfo->nPixelSpace;
84 16676 : if (psInfo->eMergeAlg == GRMA_Add)
85 : {
86 280 : if (psInfo->poSetVisitedPoints)
87 : {
88 280 : CPLAssert(!psInfo->bFillSetVisitedPoints);
89 280 : uint64_t nKey = MakeKey(nY, nXStart);
90 280 : auto &oSetVisitedPoints = *(psInfo->poSetVisitedPoints);
91 3791 : for (int nX = nXStart; nX <= nXEnd; ++nX)
92 : {
93 3511 : if (oSetVisitedPoints.find(nKey) == oSetVisitedPoints.end())
94 : {
95 3379 : double dfVal = static_cast<double>(
96 3379 : *reinterpret_cast<T *>(pabyInsert)) +
97 : burnValue;
98 3379 : GDALCopyWord(dfVal, *reinterpret_cast<T *>(pabyInsert));
99 : }
100 3511 : pabyInsert += psInfo->nPixelSpace;
101 3511 : ++nKey;
102 : }
103 : }
104 : else
105 : {
106 0 : for (int nX = nXStart; nX <= nXEnd; ++nX)
107 : {
108 0 : double dfVal = static_cast<double>(
109 0 : *reinterpret_cast<T *>(pabyInsert)) +
110 : burnValue;
111 0 : GDALCopyWord(dfVal, *reinterpret_cast<T *>(pabyInsert));
112 0 : pabyInsert += psInfo->nPixelSpace;
113 : }
114 : }
115 : }
116 : else
117 : {
118 : T nVal;
119 16396 : GDALCopyWord(burnValue, nVal);
120 11800192 : for (int nX = nXStart; nX <= nXEnd; ++nX)
121 : {
122 11783849 : *reinterpret_cast<T *>(pabyInsert) = nVal;
123 11783849 : pabyInsert += psInfo->nPixelSpace;
124 : }
125 : }
126 : }
127 15686 : }
128 :
129 3 : static inline void gvBurnScanlineInt64UserBurnValue(GDALRasterizeInfo *psInfo,
130 : int nY, int nXStart,
131 : int nXEnd)
132 :
133 : {
134 6 : for (int iBand = 0; iBand < psInfo->nBands; iBand++)
135 : {
136 3 : const std::int64_t burnValue = psInfo->burnValues.int64_values[iBand];
137 :
138 3 : unsigned char *pabyInsert =
139 3 : psInfo->pabyChunkBuf + iBand * psInfo->nBandSpace +
140 3 : nY * psInfo->nLineSpace + nXStart * psInfo->nPixelSpace;
141 3 : if (psInfo->eMergeAlg == GRMA_Add)
142 : {
143 0 : if (psInfo->poSetVisitedPoints)
144 : {
145 0 : CPLAssert(!psInfo->bFillSetVisitedPoints);
146 0 : uint64_t nKey = MakeKey(nY, nXStart);
147 0 : auto &oSetVisitedPoints = *(psInfo->poSetVisitedPoints);
148 0 : for (int nX = nXStart; nX <= nXEnd; ++nX)
149 : {
150 0 : if (oSetVisitedPoints.find(nKey) == oSetVisitedPoints.end())
151 : {
152 0 : *reinterpret_cast<std::int64_t *>(pabyInsert) =
153 0 : SaturatedAddSigned(
154 : *reinterpret_cast<std::int64_t *>(pabyInsert),
155 : burnValue);
156 : }
157 0 : pabyInsert += psInfo->nPixelSpace;
158 0 : ++nKey;
159 : }
160 : }
161 : else
162 : {
163 0 : for (int nX = nXStart; nX <= nXEnd; ++nX)
164 : {
165 0 : *reinterpret_cast<std::int64_t *>(pabyInsert) =
166 0 : SaturatedAddSigned(
167 : *reinterpret_cast<std::int64_t *>(pabyInsert),
168 : burnValue);
169 0 : pabyInsert += psInfo->nPixelSpace;
170 : }
171 : }
172 : }
173 : else
174 : {
175 9 : for (int nX = nXStart; nX <= nXEnd; ++nX)
176 : {
177 6 : *reinterpret_cast<std::int64_t *>(pabyInsert) = burnValue;
178 6 : pabyInsert += psInfo->nPixelSpace;
179 : }
180 : }
181 : }
182 3 : }
183 :
184 : /************************************************************************/
185 : /* gvBurnScanline() */
186 : /************************************************************************/
187 15868 : static void gvBurnScanline(GDALRasterizeInfo *psInfo, int nY, int nXStart,
188 : int nXEnd, double dfVariant)
189 :
190 : {
191 15868 : if (nXStart > nXEnd)
192 179 : return;
193 :
194 15689 : CPLAssert(nY >= 0 && nY < psInfo->nYSize);
195 15689 : CPLAssert(nXStart < psInfo->nXSize);
196 15689 : CPLAssert(nXEnd >= 0);
197 :
198 15689 : if (nXStart < 0)
199 6370 : nXStart = 0;
200 15689 : if (nXEnd >= psInfo->nXSize)
201 337 : nXEnd = psInfo->nXSize - 1;
202 :
203 15689 : if (psInfo->eBurnValueType == GDT_Int64)
204 : {
205 3 : if (psInfo->eType == GDT_Int64 &&
206 3 : psInfo->eBurnValueSource == GBV_UserBurnValue)
207 : {
208 3 : gvBurnScanlineInt64UserBurnValue(psInfo, nY, nXStart, nXEnd);
209 : }
210 : else
211 : {
212 0 : CPLAssert(false);
213 : }
214 3 : return;
215 : }
216 :
217 15686 : switch (psInfo->eType)
218 : {
219 9039 : case GDT_UInt8:
220 9039 : gvBurnScanlineBasic<GByte>(psInfo, nY, nXStart, nXEnd, dfVariant);
221 9039 : break;
222 0 : case GDT_Int8:
223 0 : gvBurnScanlineBasic<GInt8>(psInfo, nY, nXStart, nXEnd, dfVariant);
224 0 : break;
225 329 : case GDT_Int16:
226 329 : gvBurnScanlineBasic<GInt16>(psInfo, nY, nXStart, nXEnd, dfVariant);
227 329 : break;
228 0 : case GDT_UInt16:
229 0 : gvBurnScanlineBasic<GUInt16>(psInfo, nY, nXStart, nXEnd, dfVariant);
230 0 : break;
231 0 : case GDT_Int32:
232 0 : gvBurnScanlineBasic<GInt32>(psInfo, nY, nXStart, nXEnd, dfVariant);
233 0 : break;
234 0 : case GDT_UInt32:
235 0 : gvBurnScanlineBasic<GUInt32>(psInfo, nY, nXStart, nXEnd, dfVariant);
236 0 : break;
237 0 : case GDT_Int64:
238 0 : gvBurnScanlineBasic<std::int64_t>(psInfo, nY, nXStart, nXEnd,
239 : dfVariant);
240 0 : break;
241 0 : case GDT_UInt64:
242 0 : gvBurnScanlineBasic<std::uint64_t>(psInfo, nY, nXStart, nXEnd,
243 : dfVariant);
244 0 : break;
245 0 : case GDT_Float16:
246 0 : gvBurnScanlineBasic<GFloat16>(psInfo, nY, nXStart, nXEnd,
247 : dfVariant);
248 0 : break;
249 263 : case GDT_Float32:
250 263 : gvBurnScanlineBasic<float>(psInfo, nY, nXStart, nXEnd, dfVariant);
251 263 : break;
252 6055 : case GDT_Float64:
253 6055 : gvBurnScanlineBasic<double>(psInfo, nY, nXStart, nXEnd, dfVariant);
254 6055 : break;
255 0 : case GDT_CInt16:
256 : case GDT_CInt32:
257 : case GDT_CFloat16:
258 : case GDT_CFloat32:
259 : case GDT_CFloat64:
260 : case GDT_Unknown:
261 : case GDT_TypeCount:
262 0 : CPLAssert(false);
263 : break;
264 : }
265 : }
266 :
267 : /************************************************************************/
268 : /* gvBurnPointBasic() */
269 : /************************************************************************/
270 : template <typename T>
271 45375 : static inline void gvBurnPointBasic(GDALRasterizeInfo *psInfo, int nY, int nX,
272 : double dfVariant)
273 :
274 : {
275 95427 : for (int iBand = 0; iBand < psInfo->nBands; iBand++)
276 : {
277 50052 : double burnValue =
278 50052 : (psInfo->burnValues.double_values[iBand] +
279 50052 : ((psInfo->eBurnValueSource == GBV_UserBurnValue) ? 0 : dfVariant));
280 50052 : unsigned char *pbyInsert =
281 50052 : psInfo->pabyChunkBuf + iBand * psInfo->nBandSpace +
282 50052 : nY * psInfo->nLineSpace + nX * psInfo->nPixelSpace;
283 :
284 50052 : T *pbyPixel = reinterpret_cast<T *>(pbyInsert);
285 50052 : if (psInfo->eMergeAlg == GRMA_Add)
286 1020 : burnValue += static_cast<double>(*pbyPixel);
287 50052 : GDALCopyWord(burnValue, *pbyPixel);
288 : }
289 45375 : }
290 :
291 0 : static inline void gvBurnPointInt64UserBurnValue(GDALRasterizeInfo *psInfo,
292 : int nY, int nX)
293 :
294 : {
295 0 : for (int iBand = 0; iBand < psInfo->nBands; iBand++)
296 : {
297 0 : std::int64_t burnValue = psInfo->burnValues.int64_values[iBand];
298 0 : unsigned char *pbyInsert =
299 0 : psInfo->pabyChunkBuf + iBand * psInfo->nBandSpace +
300 0 : nY * psInfo->nLineSpace + nX * psInfo->nPixelSpace;
301 :
302 0 : std::int64_t *pbyPixel = reinterpret_cast<std::int64_t *>(pbyInsert);
303 0 : if (psInfo->eMergeAlg == GRMA_Add)
304 : {
305 0 : burnValue = SaturatedAddSigned(burnValue, *pbyPixel);
306 : }
307 0 : *pbyPixel = burnValue;
308 : }
309 0 : }
310 :
311 : /************************************************************************/
312 : /* gvBurnPoint() */
313 : /************************************************************************/
314 45387 : static void gvBurnPoint(GDALRasterizeInfo *psInfo, int nY, int nX,
315 : double dfVariant)
316 :
317 : {
318 :
319 45387 : CPLAssert(nY >= 0 && nY < psInfo->nYSize);
320 45387 : CPLAssert(nX >= 0 && nX < psInfo->nXSize);
321 :
322 45387 : if (psInfo->poSetVisitedPoints)
323 : {
324 538 : const uint64_t nKey = MakeKey(nY, nX);
325 538 : if (psInfo->poSetVisitedPoints->find(nKey) ==
326 1076 : psInfo->poSetVisitedPoints->end())
327 : {
328 526 : if (psInfo->bFillSetVisitedPoints)
329 526 : psInfo->poSetVisitedPoints->insert(nKey);
330 : }
331 : else
332 : {
333 12 : return;
334 : }
335 : }
336 :
337 45375 : if (psInfo->eBurnValueType == GDT_Int64)
338 : {
339 0 : if (psInfo->eType == GDT_Int64 &&
340 0 : psInfo->eBurnValueSource == GBV_UserBurnValue)
341 : {
342 0 : gvBurnPointInt64UserBurnValue(psInfo, nY, nX);
343 : }
344 : else
345 : {
346 0 : CPLAssert(false);
347 : }
348 0 : return;
349 : }
350 :
351 45375 : switch (psInfo->eType)
352 : {
353 18119 : case GDT_UInt8:
354 18119 : gvBurnPointBasic<GByte>(psInfo, nY, nX, dfVariant);
355 18119 : break;
356 0 : case GDT_Int8:
357 0 : gvBurnPointBasic<GInt8>(psInfo, nY, nX, dfVariant);
358 0 : break;
359 0 : case GDT_Int16:
360 0 : gvBurnPointBasic<GInt16>(psInfo, nY, nX, dfVariant);
361 0 : break;
362 0 : case GDT_UInt16:
363 0 : gvBurnPointBasic<GUInt16>(psInfo, nY, nX, dfVariant);
364 0 : break;
365 0 : case GDT_Int32:
366 0 : gvBurnPointBasic<GInt32>(psInfo, nY, nX, dfVariant);
367 0 : break;
368 0 : case GDT_UInt32:
369 0 : gvBurnPointBasic<GUInt32>(psInfo, nY, nX, dfVariant);
370 0 : break;
371 0 : case GDT_Int64:
372 0 : gvBurnPointBasic<std::int64_t>(psInfo, nY, nX, dfVariant);
373 0 : break;
374 0 : case GDT_UInt64:
375 0 : gvBurnPointBasic<std::uint64_t>(psInfo, nY, nX, dfVariant);
376 0 : break;
377 0 : case GDT_Float16:
378 0 : gvBurnPointBasic<GFloat16>(psInfo, nY, nX, dfVariant);
379 0 : break;
380 732 : case GDT_Float32:
381 732 : gvBurnPointBasic<float>(psInfo, nY, nX, dfVariant);
382 732 : break;
383 26524 : case GDT_Float64:
384 26524 : gvBurnPointBasic<double>(psInfo, nY, nX, dfVariant);
385 26524 : break;
386 0 : case GDT_CInt16:
387 : case GDT_CInt32:
388 : case GDT_CFloat16:
389 : case GDT_CFloat32:
390 : case GDT_CFloat64:
391 : case GDT_Unknown:
392 : case GDT_TypeCount:
393 0 : CPLAssert(false);
394 : }
395 : }
396 :
397 : /************************************************************************/
398 : /* GDALCollectRingsFromGeometry() */
399 : /************************************************************************/
400 :
401 2896 : static void GDALCollectRingsFromGeometry(const OGRGeometry *poShape,
402 : std::vector<double> &aPointX,
403 : std::vector<double> &aPointY,
404 : std::vector<double> &aPointVariant,
405 : std::vector<int> &aPartSize,
406 : GDALBurnValueSrc eBurnValueSrc)
407 :
408 : {
409 2896 : if (poShape == nullptr || poShape->IsEmpty())
410 0 : return;
411 :
412 2896 : const OGRwkbGeometryType eFlatType = wkbFlatten(poShape->getGeometryType());
413 :
414 2896 : if (eFlatType == wkbPoint)
415 : {
416 5 : const auto poPoint = poShape->toPoint();
417 :
418 5 : aPointX.push_back(poPoint->getX());
419 5 : aPointY.push_back(poPoint->getY());
420 5 : aPartSize.push_back(1);
421 5 : if (eBurnValueSrc != GBV_UserBurnValue)
422 : {
423 : // TODO(schwehr): Why not have the option for M r18164?
424 : // switch( eBurnValueSrc )
425 : // {
426 : // case GBV_Z:*/
427 0 : aPointVariant.push_back(poPoint->getZ());
428 : // break;
429 : // case GBV_M:
430 : // aPointVariant.reserve( nNewCount );
431 : // aPointVariant.push_back( poPoint->getM() );
432 : }
433 : }
434 2891 : else if (EQUAL(poShape->getGeometryName(), "LINEARRING"))
435 : {
436 508 : const auto poRing = poShape->toLinearRing();
437 508 : const int nCount = poRing->getNumPoints();
438 508 : const size_t nNewCount = aPointX.size() + static_cast<size_t>(nCount);
439 :
440 508 : aPointX.reserve(nNewCount);
441 508 : aPointY.reserve(nNewCount);
442 508 : if (eBurnValueSrc != GBV_UserBurnValue)
443 7 : aPointVariant.reserve(nNewCount);
444 508 : if (poRing->isClockwise())
445 : {
446 7916 : for (int i = 0; i < nCount; i++)
447 : {
448 7592 : aPointX.push_back(poRing->getX(i));
449 7592 : aPointY.push_back(poRing->getY(i));
450 7592 : if (eBurnValueSrc != GBV_UserBurnValue)
451 : {
452 : /*switch( eBurnValueSrc )
453 : {
454 : case GBV_Z:*/
455 19 : aPointVariant.push_back(poRing->getZ(i));
456 : /*break;
457 : case GBV_M:
458 : aPointVariant.push_back( poRing->getM(i) );
459 : }*/
460 : }
461 : }
462 : }
463 : else
464 : {
465 6092 : for (int i = nCount - 1; i >= 0; i--)
466 : {
467 5908 : aPointX.push_back(poRing->getX(i));
468 5908 : aPointY.push_back(poRing->getY(i));
469 5908 : if (eBurnValueSrc != GBV_UserBurnValue)
470 : {
471 : /*switch( eBurnValueSrc )
472 : {
473 : case GBV_Z:*/
474 15 : aPointVariant.push_back(poRing->getZ(i));
475 : /*break;
476 : case GBV_M:
477 : aPointVariant.push_back( poRing->getM(i) );
478 : }*/
479 : }
480 : }
481 : }
482 508 : aPartSize.push_back(nCount);
483 : }
484 2383 : else if (eFlatType == wkbLineString)
485 : {
486 1886 : const auto poLine = poShape->toLineString();
487 1886 : const int nCount = poLine->getNumPoints();
488 1886 : const size_t nNewCount = aPointX.size() + static_cast<size_t>(nCount);
489 :
490 1886 : aPointX.reserve(nNewCount);
491 1886 : aPointY.reserve(nNewCount);
492 1886 : if (eBurnValueSrc != GBV_UserBurnValue)
493 1856 : aPointVariant.reserve(nNewCount);
494 50769 : for (int i = nCount - 1; i >= 0; i--)
495 : {
496 48883 : aPointX.push_back(poLine->getX(i));
497 48883 : aPointY.push_back(poLine->getY(i));
498 48883 : if (eBurnValueSrc != GBV_UserBurnValue)
499 : {
500 : /*switch( eBurnValueSrc )
501 : {
502 : case GBV_Z:*/
503 48614 : aPointVariant.push_back(poLine->getZ(i));
504 : /*break;
505 : case GBV_M:
506 : aPointVariant.push_back( poLine->getM(i) );
507 : }*/
508 : }
509 : }
510 1886 : aPartSize.push_back(nCount);
511 : }
512 497 : else if (eFlatType == wkbPolygon)
513 : {
514 497 : const auto poPolygon = poShape->toPolygon();
515 :
516 497 : GDALCollectRingsFromGeometry(poPolygon->getExteriorRing(), aPointX,
517 : aPointY, aPointVariant, aPartSize,
518 : eBurnValueSrc);
519 :
520 508 : for (int i = 0; i < poPolygon->getNumInteriorRings(); i++)
521 11 : GDALCollectRingsFromGeometry(poPolygon->getInteriorRing(i), aPointX,
522 : aPointY, aPointVariant, aPartSize,
523 : eBurnValueSrc);
524 : }
525 0 : else if (eFlatType == wkbMultiPoint || eFlatType == wkbMultiLineString ||
526 0 : eFlatType == wkbMultiPolygon || eFlatType == wkbGeometryCollection)
527 : {
528 0 : const auto poGC = poShape->toGeometryCollection();
529 0 : for (int i = 0; i < poGC->getNumGeometries(); i++)
530 0 : GDALCollectRingsFromGeometry(poGC->getGeometryRef(i), aPointX,
531 : aPointY, aPointVariant, aPartSize,
532 0 : eBurnValueSrc);
533 : }
534 : else
535 : {
536 0 : CPLDebug("GDAL", "Rasterizer ignoring non-polygonal geometry.");
537 : }
538 : }
539 :
540 : /************************************************************************
541 : * gv_rasterize_one_shape()
542 : *
543 : * @param pabyChunkBuf buffer to which values will be burned
544 : * @param nXOff chunk column offset from left edge of raster
545 : * @param nYOff chunk scanline offset from top of raster
546 : * @param nXSize number of columns in chunk
547 : * @param nYSize number of rows in chunk
548 : * @param nBands number of bands in chunk
549 : * @param eType data type of pabyChunkBuf
550 : * @param nPixelSpace number of bytes between adjacent pixels in chunk
551 : * (0 to calculate automatically)
552 : * @param nLineSpace number of bytes between adjacent scanlines in chunk
553 : * (0 to calculate automatically)
554 : * @param nBandSpace number of bytes between adjacent bands in chunk
555 : * (0 to calculate automatically)
556 : * @param bAllTouched burn value to all touched pixels?
557 : * @param poShape geometry to rasterize, in original coordinates.
558 : * Since GDAL 3.14, curved geometries will be linearized before
559 : * rasterization. (In previous versions, they are ignored.)
560 : * @param eBurnValueType type of value to be burned (must be Float64 or Int64)
561 : * @param padfBurnValues array of nBands values to burn (Float64), or nullptr
562 : * @param panBurnValues array of nBands values to burn (Int64), or nullptr
563 : * @param eBurnValueSrc whether to burn values from padfBurnValues /
564 : * panBurnValues, or from the Z or M values of poShape
565 : * @param eMergeAlg whether the burn value should replace or be added to the
566 : * existing values
567 : * @param pfnTransformer transformer from CRS of geometry to pixel/line
568 : * coordinates of raster
569 : * @param pTransformArg arguments to pass to pfnTransformer
570 : ************************************************************************/
571 2417 : static void gv_rasterize_one_shape(
572 : unsigned char *pabyChunkBuf, int nXOff, int nYOff, int nXSize, int nYSize,
573 : int nBands, GDALDataType eType, int nPixelSpace, GSpacing nLineSpace,
574 : GSpacing nBandSpace, int bAllTouched, const OGRGeometry *poShape,
575 : GDALDataType eBurnValueType, const double *padfBurnValues,
576 : const int64_t *panBurnValues, GDALBurnValueSrc eBurnValueSrc,
577 : GDALRasterMergeAlg eMergeAlg, GDALTransformerFunc pfnTransformer,
578 : void *pTransformArg)
579 :
580 : {
581 0 : std::unique_ptr<OGRGeometry> poLinearizedShape;
582 :
583 2417 : if (poShape == nullptr || poShape->IsEmpty())
584 1 : return;
585 :
586 2416 : if (poShape->hasCurveGeometry())
587 : {
588 5 : poLinearizedShape.reset(poShape->getLinearGeometry());
589 5 : if (!poLinearizedShape)
590 : {
591 0 : CPLError(CE_Failure, CPLE_AppDefined,
592 : "Failed to linearize geometry");
593 0 : return;
594 : }
595 5 : poShape = poLinearizedShape.get();
596 : }
597 :
598 2416 : const auto eGeomType = wkbFlatten(poShape->getGeometryType());
599 :
600 2416 : if ((eGeomType == wkbMultiLineString || eGeomType == wkbMultiPolygon ||
601 28 : eGeomType == wkbGeometryCollection) &&
602 : eMergeAlg == GRMA_Replace)
603 : {
604 : // Speed optimization: in replace mode, we can rasterize each part of
605 : // a geometry collection separately.
606 28 : const auto poGC = poShape->toGeometryCollection();
607 65 : for (const auto poPart : *poGC)
608 : {
609 37 : gv_rasterize_one_shape(
610 : pabyChunkBuf, nXOff, nYOff, nXSize, nYSize, nBands, eType,
611 : nPixelSpace, nLineSpace, nBandSpace, bAllTouched, poPart,
612 : eBurnValueType, padfBurnValues, panBurnValues, eBurnValueSrc,
613 : eMergeAlg, pfnTransformer, pTransformArg);
614 : }
615 28 : return;
616 : }
617 :
618 2388 : if (nPixelSpace == 0)
619 : {
620 2388 : nPixelSpace = GDALGetDataTypeSizeBytes(eType);
621 : }
622 2388 : if (nLineSpace == 0)
623 : {
624 2388 : nLineSpace = static_cast<GSpacing>(nXSize) * nPixelSpace;
625 : }
626 2388 : if (nBandSpace == 0)
627 : {
628 2388 : nBandSpace = nYSize * nLineSpace;
629 : }
630 :
631 : GDALRasterizeInfo sInfo;
632 2388 : sInfo.nXSize = nXSize;
633 2388 : sInfo.nYSize = nYSize;
634 2388 : sInfo.nBands = nBands;
635 2388 : sInfo.pabyChunkBuf = pabyChunkBuf;
636 2388 : sInfo.eType = eType;
637 2388 : sInfo.nPixelSpace = nPixelSpace;
638 2388 : sInfo.nLineSpace = nLineSpace;
639 2388 : sInfo.nBandSpace = nBandSpace;
640 2388 : sInfo.eBurnValueType = eBurnValueType;
641 2388 : if (eBurnValueType == GDT_Float64)
642 2387 : sInfo.burnValues.double_values = padfBurnValues;
643 1 : else if (eBurnValueType == GDT_Int64)
644 1 : sInfo.burnValues.int64_values = panBurnValues;
645 : else
646 : {
647 0 : CPLAssert(false);
648 : }
649 2388 : sInfo.eBurnValueSource = eBurnValueSrc;
650 2388 : sInfo.eMergeAlg = eMergeAlg;
651 2388 : sInfo.bFillSetVisitedPoints = false;
652 2388 : sInfo.poSetVisitedPoints = nullptr;
653 :
654 : /* -------------------------------------------------------------------- */
655 : /* Transform polygon geometries into a set of rings and a part */
656 : /* size list. */
657 : /* -------------------------------------------------------------------- */
658 : std::vector<double>
659 4776 : aPointX; // coordinate X values from all rings/components
660 : std::vector<double>
661 4776 : aPointY; // coordinate Y values from all rings/components
662 4776 : std::vector<double> aPointVariant; // coordinate Z values
663 4776 : std::vector<int> aPartSize; // number of X/Y/(Z) values associated with
664 : // each ring/component
665 :
666 2388 : GDALCollectRingsFromGeometry(poShape, aPointX, aPointY, aPointVariant,
667 : aPartSize, eBurnValueSrc);
668 :
669 : /* -------------------------------------------------------------------- */
670 : /* Transform points if needed. */
671 : /* -------------------------------------------------------------------- */
672 2388 : if (pfnTransformer != nullptr)
673 : {
674 : int *panSuccess =
675 2388 : static_cast<int *>(CPLCalloc(sizeof(int), aPointX.size()));
676 :
677 : // TODO: We need to add all appropriate error checking at some point.
678 2388 : pfnTransformer(pTransformArg, FALSE, static_cast<int>(aPointX.size()),
679 : aPointX.data(), aPointY.data(), nullptr, panSuccess);
680 2388 : CPLFree(panSuccess);
681 : }
682 :
683 : /* -------------------------------------------------------------------- */
684 : /* Shift to account for the buffer offset of this buffer. */
685 : /* -------------------------------------------------------------------- */
686 64776 : for (unsigned int i = 0; i < aPointX.size(); i++)
687 62388 : aPointX[i] -= nXOff;
688 64776 : for (unsigned int i = 0; i < aPointY.size(); i++)
689 62388 : aPointY[i] -= nYOff;
690 :
691 : /* -------------------------------------------------------------------- */
692 : /* Perform the rasterization. */
693 : /* According to the C++ Standard/23.2.4, elements of a vector are */
694 : /* stored in continuous memory block. */
695 : /* -------------------------------------------------------------------- */
696 :
697 2388 : switch (eGeomType)
698 : {
699 5 : case wkbPoint:
700 : case wkbMultiPoint:
701 10 : GDALdllImagePoint(
702 5 : sInfo.nXSize, nYSize, static_cast<int>(aPartSize.size()),
703 5 : aPartSize.data(), aPointX.data(), aPointY.data(),
704 : (eBurnValueSrc == GBV_UserBurnValue) ? nullptr
705 0 : : aPointVariant.data(),
706 : gvBurnPoint, &sInfo);
707 5 : break;
708 1886 : case wkbLineString:
709 : case wkbMultiLineString:
710 : {
711 1886 : if (eMergeAlg == GRMA_Add)
712 : {
713 10 : sInfo.bFillSetVisitedPoints = true;
714 10 : sInfo.poSetVisitedPoints = new std::set<uint64_t>();
715 : }
716 1886 : if (bAllTouched)
717 20 : GDALdllImageLineAllTouched(
718 10 : sInfo.nXSize, nYSize, static_cast<int>(aPartSize.size()),
719 10 : aPartSize.data(), aPointX.data(), aPointY.data(),
720 : (eBurnValueSrc == GBV_UserBurnValue) ? nullptr
721 0 : : aPointVariant.data(),
722 : gvBurnPoint, &sInfo, eMergeAlg == GRMA_Add, false);
723 : else
724 3752 : GDALdllImageLine(
725 1876 : sInfo.nXSize, nYSize, static_cast<int>(aPartSize.size()),
726 1876 : aPartSize.data(), aPointX.data(), aPointY.data(),
727 : (eBurnValueSrc == GBV_UserBurnValue) ? nullptr
728 1856 : : aPointVariant.data(),
729 : gvBurnPoint, &sInfo);
730 : }
731 1886 : break;
732 :
733 497 : default:
734 : {
735 497 : if (eMergeAlg == GRMA_Add)
736 : {
737 15 : sInfo.bFillSetVisitedPoints = true;
738 15 : sInfo.poSetVisitedPoints = new std::set<uint64_t>();
739 : }
740 497 : if (bAllTouched)
741 : {
742 : // Reverting the variants to the first value because the
743 : // polygon is filled using the variant from the first point of
744 : // the first segment. Should be removed when the code to full
745 : // polygons more appropriately is added.
746 142 : if (eBurnValueSrc == GBV_UserBurnValue)
747 : {
748 139 : GDALdllImageLineAllTouched(
749 : sInfo.nXSize, nYSize,
750 139 : static_cast<int>(aPartSize.size()), aPartSize.data(),
751 139 : aPointX.data(), aPointY.data(), nullptr, gvBurnPoint,
752 : &sInfo, eMergeAlg == GRMA_Add, true);
753 : }
754 : else
755 : {
756 6 : for (unsigned int i = 0, n = 0;
757 6 : i < static_cast<unsigned int>(aPartSize.size()); i++)
758 : {
759 17 : for (int j = 0; j < aPartSize[i]; j++)
760 14 : aPointVariant[n++] = aPointVariant[0];
761 : }
762 :
763 3 : GDALdllImageLineAllTouched(
764 : sInfo.nXSize, nYSize,
765 3 : static_cast<int>(aPartSize.size()), aPartSize.data(),
766 3 : aPointX.data(), aPointY.data(), aPointVariant.data(),
767 : gvBurnPoint, &sInfo, eMergeAlg == GRMA_Add, true);
768 : }
769 : }
770 497 : sInfo.bFillSetVisitedPoints = false;
771 994 : GDALdllImageFilledPolygon(
772 497 : sInfo.nXSize, nYSize, static_cast<int>(aPartSize.size()),
773 497 : aPartSize.data(), aPointX.data(), aPointY.data(),
774 : (eBurnValueSrc == GBV_UserBurnValue) ? nullptr
775 6 : : aPointVariant.data(),
776 : gvBurnScanline, &sInfo, eMergeAlg == GRMA_Add);
777 : }
778 497 : break;
779 : }
780 :
781 2388 : delete sInfo.poSetVisitedPoints;
782 : }
783 :
784 : /************************************************************************/
785 : /* GDALRasterizeOptions() */
786 : /* */
787 : /* Recognise a few rasterize options used by all three entry */
788 : /* points. */
789 : /************************************************************************/
790 :
791 351 : static CPLErr GDALRasterizeOptions(CSLConstList papszOptions, int *pbAllTouched,
792 : GDALBurnValueSrc *peBurnValueSource,
793 : GDALRasterMergeAlg *peMergeAlg,
794 : GDALRasterizeOptim *peOptim)
795 : {
796 351 : *pbAllTouched = CPLFetchBool(papszOptions, "ALL_TOUCHED", false);
797 :
798 351 : const char *pszOpt = CSLFetchNameValue(papszOptions, "BURN_VALUE_FROM");
799 351 : *peBurnValueSource = GBV_UserBurnValue;
800 351 : if (pszOpt)
801 : {
802 5 : if (EQUAL(pszOpt, "Z"))
803 : {
804 5 : *peBurnValueSource = GBV_Z;
805 : }
806 : // else if( EQUAL(pszOpt, "M"))
807 : // eBurnValueSource = GBV_M;
808 : else
809 : {
810 0 : CPLError(CE_Failure, CPLE_AppDefined,
811 : "Unrecognized value '%s' for BURN_VALUE_FROM.", pszOpt);
812 0 : return CE_Failure;
813 : }
814 : }
815 :
816 : /* -------------------------------------------------------------------- */
817 : /* MERGE_ALG=[REPLACE]/ADD */
818 : /* -------------------------------------------------------------------- */
819 351 : *peMergeAlg = GRMA_Replace;
820 351 : pszOpt = CSLFetchNameValue(papszOptions, "MERGE_ALG");
821 351 : if (pszOpt)
822 : {
823 18 : if (EQUAL(pszOpt, "ADD"))
824 : {
825 18 : *peMergeAlg = GRMA_Add;
826 : }
827 0 : else if (EQUAL(pszOpt, "REPLACE"))
828 : {
829 0 : *peMergeAlg = GRMA_Replace;
830 : }
831 : else
832 : {
833 0 : CPLError(CE_Failure, CPLE_AppDefined,
834 : "Unrecognized value '%s' for MERGE_ALG.", pszOpt);
835 0 : return CE_Failure;
836 : }
837 : }
838 :
839 : /* -------------------------------------------------------------------- */
840 : /* OPTIM=[AUTO]/RASTER/VECTOR */
841 : /* -------------------------------------------------------------------- */
842 351 : pszOpt = CSLFetchNameValue(papszOptions, "OPTIM");
843 351 : if (pszOpt)
844 : {
845 39 : if (peOptim)
846 : {
847 39 : *peOptim = GRO_Auto;
848 39 : if (EQUAL(pszOpt, "RASTER"))
849 : {
850 1 : *peOptim = GRO_Raster;
851 : }
852 38 : else if (EQUAL(pszOpt, "VECTOR"))
853 : {
854 1 : *peOptim = GRO_Vector;
855 : }
856 37 : else if (EQUAL(pszOpt, "AUTO"))
857 : {
858 37 : *peOptim = GRO_Auto;
859 : }
860 : else
861 : {
862 0 : CPLError(CE_Failure, CPLE_AppDefined,
863 : "Unrecognized value '%s' for OPTIM.", pszOpt);
864 0 : return CE_Failure;
865 : }
866 : }
867 : else
868 : {
869 0 : CPLError(CE_Warning, CPLE_NotSupported,
870 : "Option OPTIM is not supported by this function");
871 : }
872 : }
873 :
874 351 : return CE_None;
875 : }
876 :
877 : /************************************************************************/
878 : /* GDALRasterizeGeometries() */
879 : /************************************************************************/
880 :
881 : static CPLErr GDALRasterizeGeometriesInternal(
882 : GDALDatasetH hDS, int nBandCount, const int *panBandList, int nGeomCount,
883 : const OGRGeometryH *pahGeometries, GDALTransformerFunc pfnTransformer,
884 : void *pTransformArg, GDALDataType eBurnValueType,
885 : const double *padfGeomBurnValues, const int64_t *panGeomBurnValues,
886 : CSLConstList papszOptions, GDALProgressFunc pfnProgress,
887 : void *pProgressArg);
888 :
889 : /**
890 : * Burn geometries into raster.
891 : *
892 : * Rasterize a list of geometric objects into a raster dataset. The
893 : * geometries are passed as an array of OGRGeometryH handlers.
894 : *
895 : * If the geometries are in the georeferenced coordinates of the raster
896 : * dataset, then the pfnTransform may be passed in NULL and one will be
897 : * derived internally from the geotransform of the dataset. The transform
898 : * needs to transform the geometry locations into pixel/line coordinates
899 : * on the raster dataset.
900 : *
901 : * The output raster may be of any GDAL supported datatype. An explicit list
902 : * of burn values for each geometry for each band must be passed in.
903 : *
904 : * The papszOption list of options currently only supports one option. The
905 : * "ALL_TOUCHED" option may be enabled by setting it to "TRUE".
906 : *
907 : * @param hDS output data, must be opened in update mode.
908 : * @param nBandCount the number of bands to be updated.
909 : * @param panBandList the list of bands to be updated.
910 : * @param nGeomCount the number of geometries being passed in pahGeometries.
911 : * @param pahGeometries the array of geometries to burn in. Since GDAL 3.14,
912 : * curved geometries will be converted to linear types.
913 : * @param pfnTransformer transformation to apply to geometries to put into
914 : * pixel/line coordinates on raster. If NULL a geotransform based one will
915 : * be created internally.
916 : * @param pTransformArg callback data for transformer.
917 : * @param padfGeomBurnValues the array of values to burn into the raster.
918 : * There should be nBandCount values for each geometry.
919 : * @param papszOptions special options controlling rasterization
920 : * <ul>
921 : * <li>"ALL_TOUCHED": May be set to TRUE to set all pixels touched
922 : * by the line or polygons, not just those whose center is within the polygon
923 : * (behavior is unspecified when the polygon is just touching the pixel center)
924 : * or that are selected by Brezenham's line algorithm. Defaults to FALSE.</li>
925 : * <li>"BURN_VALUE_FROM": May be set to "Z" to use the Z values of the
926 : * geometries. dfBurnValue is added to this before burning.
927 : * Defaults to GDALBurnValueSrc.GBV_UserBurnValue in which case just the
928 : * dfBurnValue is burned. This is implemented only for points and lines for
929 : * now. The M value may be supported in the future.</li>
930 : * <li>"MERGE_ALG": May be REPLACE (the default) or ADD. REPLACE results in
931 : * overwriting of value, while ADD adds the new value to the existing raster,
932 : * suitable for heatmaps for instance.</li>
933 : * <li>"CHUNKYSIZE": The height in lines of the chunk to operate on.
934 : * The larger the chunk size the less times we need to make a pass through all
935 : * the shapes. If it is not set or set to zero the default chunk size will be
936 : * used. Default size will be estimated based on the GDAL cache buffer size
937 : * using formula: cache_size_bytes/scanline_size_bytes, so the chunk will
938 : * not exceed the cache. Not used in OPTIM=RASTER mode.</li>
939 : * <li>"OPTIM": May be set to "AUTO", "RASTER", "VECTOR". Force the algorithm
940 : * used (results are identical). The raster mode is used in most cases and
941 : * optimise read/write operations. The vector mode is useful with a decent
942 : * amount of input features and optimize the CPU use. That mode has to be used
943 : * with tiled images to be efficient. The auto mode (the default) will chose
944 : * the algorithm based on input and output properties.
945 : * </li>
946 : * </ul>
947 : * @param pfnProgress the progress function to report completion.
948 : * @param pProgressArg callback data for progress function.
949 : *
950 : * @return CE_None on success or CE_Failure on error.
951 : *
952 : * <strong>Example</strong><br>
953 : * GDALRasterizeGeometries rasterize output to MEM Dataset :<br>
954 : * @code
955 : * int nBufXSize = 1024;
956 : * int nBufYSize = 1024;
957 : * int nBandCount = 1;
958 : * GDALDataType eType = GDT_UInt8;
959 : * int nDataTypeSize = GDALGetDataTypeSizeBytes(eType);
960 : *
961 : * void* pData = CPLCalloc( nBufXSize*nBufYSize*nBandCount, nDataTypeSize );
962 : * char memdsetpath[1024];
963 : * sprintf(memdsetpath,"MEM:::DATAPOINTER=0x%p,PIXELS=%d,LINES=%d,"
964 : * "BANDS=%d,DATATYPE=%s,PIXELOFFSET=%d,LINEOFFSET=%d",
965 : * pData,nBufXSize,nBufYSize,nBandCount,GDALGetDataTypeName(eType),
966 : * nBandCount*nDataTypeSize, nBufXSize*nBandCount*nDataTypeSize );
967 : *
968 : * // Open Memory Dataset
969 : * GDALDatasetH hMemDset = GDALOpen(memdsetpath, GA_Update);
970 : * // or create it as follows
971 : * // GDALDriverH hMemDriver = GDALGetDriverByName("MEM");
972 : * // GDALDatasetH hMemDset = GDALCreate(hMemDriver, "", nBufXSize,
973 : * nBufYSize, nBandCount, eType, NULL);
974 : *
975 : * double adfGeoTransform[6];
976 : * // Assign GeoTransform parameters,Omitted here.
977 : *
978 : * GDALSetGeoTransform(hMemDset,adfGeoTransform);
979 : * GDALSetProjection(hMemDset,pszProjection); // Can not
980 : *
981 : * // Do something ...
982 : * // Need an array of OGRGeometry objects,The assumption here is pahGeoms
983 : *
984 : * int bandList[3] = { 1, 2, 3};
985 : * std::vector<double> geomBurnValue(nGeomCount*nBandCount,255.0);
986 : * CPLErr err = GDALRasterizeGeometries(
987 : * hMemDset, nBandCount, bandList, nGeomCount, pahGeoms, pfnTransformer,
988 : * pTransformArg, geomBurnValue.data(), papszOptions,
989 : * pfnProgress, pProgressArg);
990 : * if( err != CE_None )
991 : * {
992 : * // Do something ...
993 : * }
994 : * GDALClose(hMemDset);
995 : * CPLFree(pData);
996 : *@endcode
997 : */
998 :
999 302 : CPLErr GDALRasterizeGeometries(
1000 : GDALDatasetH hDS, int nBandCount, const int *panBandList, int nGeomCount,
1001 : const OGRGeometryH *pahGeometries, GDALTransformerFunc pfnTransformer,
1002 : void *pTransformArg, const double *padfGeomBurnValues,
1003 : CSLConstList papszOptions, GDALProgressFunc pfnProgress, void *pProgressArg)
1004 :
1005 : {
1006 302 : VALIDATE_POINTER1(hDS, "GDALRasterizeGeometries", CE_Failure);
1007 :
1008 302 : return GDALRasterizeGeometriesInternal(
1009 : hDS, nBandCount, panBandList, nGeomCount, pahGeometries, pfnTransformer,
1010 : pTransformArg, GDT_Float64, padfGeomBurnValues, nullptr, papszOptions,
1011 302 : pfnProgress, pProgressArg);
1012 : }
1013 :
1014 : /**
1015 : * Burn geometries into raster.
1016 : *
1017 : * Same as GDALRasterizeGeometries(), except that the burn values array is
1018 : * of type Int64. And the datatype of the output raster *must* be GDT_Int64.
1019 : *
1020 : * @since GDAL 3.5
1021 : */
1022 1 : CPLErr GDALRasterizeGeometriesInt64(
1023 : GDALDatasetH hDS, int nBandCount, const int *panBandList, int nGeomCount,
1024 : const OGRGeometryH *pahGeometries, GDALTransformerFunc pfnTransformer,
1025 : void *pTransformArg, const int64_t *panGeomBurnValues,
1026 : CSLConstList papszOptions, GDALProgressFunc pfnProgress, void *pProgressArg)
1027 :
1028 : {
1029 1 : VALIDATE_POINTER1(hDS, "GDALRasterizeGeometriesInt64", CE_Failure);
1030 :
1031 1 : return GDALRasterizeGeometriesInternal(
1032 : hDS, nBandCount, panBandList, nGeomCount, pahGeometries, pfnTransformer,
1033 : pTransformArg, GDT_Int64, nullptr, panGeomBurnValues, papszOptions,
1034 1 : pfnProgress, pProgressArg);
1035 : }
1036 :
1037 303 : static CPLErr GDALRasterizeGeometriesInternal(
1038 : GDALDatasetH hDS, int nBandCount, const int *panBandList, int nGeomCount,
1039 : const OGRGeometryH *pahGeometries, GDALTransformerFunc pfnTransformer,
1040 : void *pTransformArg, GDALDataType eBurnValueType,
1041 : const double *padfGeomBurnValues, const int64_t *panGeomBurnValues,
1042 : CSLConstList papszOptions, GDALProgressFunc pfnProgress, void *pProgressArg)
1043 :
1044 : {
1045 303 : if (pfnProgress == nullptr)
1046 239 : pfnProgress = GDALDummyProgress;
1047 :
1048 303 : GDALDataset *poDS = GDALDataset::FromHandle(hDS);
1049 : /* -------------------------------------------------------------------- */
1050 : /* Do some rudimentary arg checking. */
1051 : /* -------------------------------------------------------------------- */
1052 303 : if (nBandCount == 0 || nGeomCount == 0)
1053 : {
1054 0 : pfnProgress(1.0, "", pProgressArg);
1055 0 : return CE_None;
1056 : }
1057 :
1058 303 : if (eBurnValueType == GDT_Int64)
1059 : {
1060 2 : for (int i = 0; i < nBandCount; i++)
1061 : {
1062 1 : GDALRasterBand *poBand = poDS->GetRasterBand(panBandList[i]);
1063 1 : if (poBand == nullptr)
1064 0 : return CE_Failure;
1065 1 : if (poBand->GetRasterDataType() != GDT_Int64)
1066 : {
1067 0 : CPLError(CE_Failure, CPLE_NotSupported,
1068 : "GDALRasterizeGeometriesInt64() only supported on "
1069 : "Int64 raster");
1070 0 : return CE_Failure;
1071 : }
1072 : }
1073 : }
1074 :
1075 : // Prototype band.
1076 303 : GDALRasterBand *poBand = poDS->GetRasterBand(panBandList[0]);
1077 303 : if (poBand == nullptr)
1078 0 : return CE_Failure;
1079 :
1080 : /* -------------------------------------------------------------------- */
1081 : /* Options */
1082 : /* -------------------------------------------------------------------- */
1083 303 : int bAllTouched = FALSE;
1084 303 : GDALBurnValueSrc eBurnValueSource = GBV_UserBurnValue;
1085 303 : GDALRasterMergeAlg eMergeAlg = GRMA_Replace;
1086 303 : GDALRasterizeOptim eOptim = GRO_Auto;
1087 303 : if (GDALRasterizeOptions(papszOptions, &bAllTouched, &eBurnValueSource,
1088 303 : &eMergeAlg, &eOptim) == CE_Failure)
1089 : {
1090 0 : return CE_Failure;
1091 : }
1092 :
1093 : /* -------------------------------------------------------------------- */
1094 : /* If we have no transformer, assume the geometries are in file */
1095 : /* georeferenced coordinates, and create a transformer to */
1096 : /* convert that to pixel/line coordinates. */
1097 : /* */
1098 : /* We really just need to apply an affine transform, but for */
1099 : /* simplicity we use the more general GenImgProjTransformer. */
1100 : /* -------------------------------------------------------------------- */
1101 303 : bool bNeedToFreeTransformer = false;
1102 :
1103 303 : if (pfnTransformer == nullptr)
1104 : {
1105 271 : bNeedToFreeTransformer = true;
1106 :
1107 271 : CPLStringList aosTransformerOptions;
1108 271 : GDALGeoTransform gt;
1109 277 : if (poDS->GetGeoTransform(gt) != CE_None && poDS->GetGCPCount() == 0 &&
1110 6 : poDS->GetMetadata("RPC") == nullptr)
1111 : {
1112 5 : aosTransformerOptions.SetNameValue("DST_METHOD", "NO_GEOTRANSFORM");
1113 : }
1114 :
1115 271 : pTransformArg = GDALCreateGenImgProjTransformer2(
1116 271 : nullptr, hDS, aosTransformerOptions.List());
1117 :
1118 271 : pfnTransformer = GDALGenImgProjTransform;
1119 271 : if (pTransformArg == nullptr)
1120 : {
1121 0 : return CE_Failure;
1122 : }
1123 : }
1124 :
1125 : /* -------------------------------------------------------------------- */
1126 : /* Choice of optimisation in auto mode. Use vector optim : */
1127 : /* 1) if output is tiled */
1128 : /* 2) if large number of features is present (>10000) */
1129 : /* 3) if the nb of pixels > 50 * nb of features (not-too-small ft) */
1130 : /* -------------------------------------------------------------------- */
1131 : int nXBlockSize, nYBlockSize;
1132 303 : poBand->GetBlockSize(&nXBlockSize, &nYBlockSize);
1133 :
1134 303 : if (eOptim == GRO_Auto)
1135 : {
1136 301 : eOptim = GRO_Raster;
1137 : // TODO make more tests with various inputs/outputs to adjust the
1138 : // parameters
1139 301 : if (nYBlockSize > 1 && nGeomCount > 10000 &&
1140 0 : (poBand->GetXSize() * static_cast<long long>(poBand->GetYSize()) /
1141 0 : nGeomCount >
1142 : 50))
1143 : {
1144 0 : eOptim = GRO_Vector;
1145 0 : CPLDebug("GDAL", "The vector optim has been chosen automatically");
1146 : }
1147 : }
1148 :
1149 : /* -------------------------------------------------------------------- */
1150 : /* The original algorithm */
1151 : /* Optimized for raster writing */
1152 : /* (optimal on a small number of large vectors) */
1153 : /* -------------------------------------------------------------------- */
1154 : unsigned char *pabyChunkBuf;
1155 303 : CPLErr eErr = CE_None;
1156 303 : if (eOptim == GRO_Raster)
1157 : {
1158 : /* --------------------------------------------------------------------
1159 : */
1160 : /* Establish a chunksize to operate on. The larger the chunk */
1161 : /* size the less times we need to make a pass through all the */
1162 : /* shapes. */
1163 : /* --------------------------------------------------------------------
1164 : */
1165 : const GDALDataType eType =
1166 302 : GDALGetNonComplexDataType(poBand->GetRasterDataType());
1167 :
1168 604 : const uint64_t nScanlineBytes = static_cast<uint64_t>(nBandCount) *
1169 302 : poDS->GetRasterXSize() *
1170 302 : GDALGetDataTypeSizeBytes(eType);
1171 :
1172 : #if SIZEOF_VOIDP < 8
1173 : // Only on 32-bit systems and in pathological cases
1174 : if (nScanlineBytes > std::numeric_limits<size_t>::max())
1175 : {
1176 : CPLError(CE_Failure, CPLE_OutOfMemory, "Too big raster");
1177 : if (bNeedToFreeTransformer)
1178 : GDALDestroyTransformer(pTransformArg);
1179 : return CE_Failure;
1180 : }
1181 : #endif
1182 :
1183 : int nYChunkSize =
1184 302 : atoi(CSLFetchNameValueDef(papszOptions, "CHUNKYSIZE", "0"));
1185 302 : if (nYChunkSize <= 0)
1186 : {
1187 302 : const GIntBig nYChunkSize64 = GDALGetCacheMax64() / nScanlineBytes;
1188 302 : const int knIntMax = std::numeric_limits<int>::max();
1189 302 : nYChunkSize = nYChunkSize64 > knIntMax
1190 302 : ? knIntMax
1191 : : static_cast<int>(nYChunkSize64);
1192 : }
1193 :
1194 302 : if (nYChunkSize < 1)
1195 0 : nYChunkSize = 1;
1196 302 : if (nYChunkSize > poDS->GetRasterYSize())
1197 302 : nYChunkSize = poDS->GetRasterYSize();
1198 :
1199 302 : CPLDebug("GDAL", "Rasterizer operating on %d swaths of %d scanlines.",
1200 302 : DIV_ROUND_UP(poDS->GetRasterYSize(), nYChunkSize),
1201 : nYChunkSize);
1202 :
1203 302 : pabyChunkBuf = static_cast<unsigned char *>(VSI_MALLOC2_VERBOSE(
1204 : nYChunkSize, static_cast<size_t>(nScanlineBytes)));
1205 302 : if (pabyChunkBuf == nullptr)
1206 : {
1207 0 : if (bNeedToFreeTransformer)
1208 0 : GDALDestroyTransformer(pTransformArg);
1209 0 : return CE_Failure;
1210 : }
1211 :
1212 : /* ====================================================================
1213 : */
1214 : /* Loop over image in designated chunks. */
1215 : /* ====================================================================
1216 : */
1217 302 : pfnProgress(0.0, nullptr, pProgressArg);
1218 :
1219 604 : for (int iY = 0; iY < poDS->GetRasterYSize() && eErr == CE_None;
1220 302 : iY += nYChunkSize)
1221 : {
1222 302 : int nThisYChunkSize = nYChunkSize;
1223 302 : if (nThisYChunkSize + iY > poDS->GetRasterYSize())
1224 0 : nThisYChunkSize = poDS->GetRasterYSize() - iY;
1225 :
1226 302 : eErr = poDS->RasterIO(
1227 : GF_Read, 0, iY, poDS->GetRasterXSize(), nThisYChunkSize,
1228 : pabyChunkBuf, poDS->GetRasterXSize(), nThisYChunkSize, eType,
1229 : nBandCount, panBandList, 0, 0, 0, nullptr);
1230 302 : if (eErr != CE_None)
1231 0 : break;
1232 :
1233 2615 : for (int iShape = 0; iShape < nGeomCount; iShape++)
1234 : {
1235 4627 : gv_rasterize_one_shape(
1236 : pabyChunkBuf, 0, iY, poDS->GetRasterXSize(),
1237 : nThisYChunkSize, nBandCount, eType, 0, 0, 0, bAllTouched,
1238 2313 : OGRGeometry::FromHandle(pahGeometries[iShape]),
1239 : eBurnValueType,
1240 : padfGeomBurnValues
1241 2312 : ? padfGeomBurnValues +
1242 2312 : static_cast<size_t>(iShape) * nBandCount
1243 : : nullptr,
1244 : panGeomBurnValues
1245 1 : ? panGeomBurnValues +
1246 1 : static_cast<size_t>(iShape) * nBandCount
1247 : : nullptr,
1248 : eBurnValueSource, eMergeAlg, pfnTransformer, pTransformArg);
1249 : }
1250 :
1251 302 : eErr = poDS->RasterIO(
1252 : GF_Write, 0, iY, poDS->GetRasterXSize(), nThisYChunkSize,
1253 : pabyChunkBuf, poDS->GetRasterXSize(), nThisYChunkSize, eType,
1254 : nBandCount, panBandList, 0, 0, 0, nullptr);
1255 :
1256 302 : if (!pfnProgress((iY + nThisYChunkSize) /
1257 302 : static_cast<double>(poDS->GetRasterYSize()),
1258 : "", pProgressArg))
1259 : {
1260 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1261 0 : eErr = CE_Failure;
1262 : }
1263 : }
1264 : }
1265 : /* -------------------------------------------------------------------- */
1266 : /* The new algorithm */
1267 : /* Optimized to minimize the vector computation */
1268 : /* (optimal on a large number of vectors & tiled raster) */
1269 : /* -------------------------------------------------------------------- */
1270 : else
1271 : {
1272 : /* --------------------------------------------------------------------
1273 : */
1274 : /* Establish a chunksize to operate on. Its size is defined by */
1275 : /* the block size of the output file. */
1276 : /* --------------------------------------------------------------------
1277 : */
1278 1 : const int nXBlocks = DIV_ROUND_UP(poBand->GetXSize(), nXBlockSize);
1279 1 : const int nYBlocks = DIV_ROUND_UP(poBand->GetYSize(), nYBlockSize);
1280 :
1281 : const GDALDataType eType =
1282 1 : poBand->GetRasterDataType() == GDT_UInt8 ? GDT_UInt8 : GDT_Float64;
1283 :
1284 1 : const int nPixelSize = nBandCount * GDALGetDataTypeSizeBytes(eType);
1285 :
1286 : // rem: optimized for square blocks
1287 : const GIntBig nbMaxBlocks64 =
1288 1 : GDALGetCacheMax64() / nPixelSize / nYBlockSize / nXBlockSize;
1289 1 : const int knIntMax = std::numeric_limits<int>::max();
1290 : const int nbMaxBlocks = static_cast<int>(
1291 2 : std::min(static_cast<GIntBig>(knIntMax / nPixelSize / nYBlockSize /
1292 : nXBlockSize),
1293 1 : nbMaxBlocks64));
1294 1 : const int nbBlocksX = std::max(
1295 2 : 1,
1296 2 : std::min(static_cast<int>(sqrt(static_cast<double>(nbMaxBlocks))),
1297 1 : nXBlocks));
1298 : const int nbBlocksY =
1299 1 : std::max(1, std::min(nbMaxBlocks / nbBlocksX, nYBlocks));
1300 :
1301 1 : const uint64_t nChunkSize = static_cast<uint64_t>(nXBlockSize) *
1302 1 : nbBlocksX * nYBlockSize * nbBlocksY;
1303 :
1304 : #if SIZEOF_VOIDP < 8
1305 : // Only on 32-bit systems and in pathological cases
1306 : if (nChunkSize > std::numeric_limits<size_t>::max())
1307 : {
1308 : CPLError(CE_Failure, CPLE_OutOfMemory, "Too big raster");
1309 : if (bNeedToFreeTransformer)
1310 : GDALDestroyTransformer(pTransformArg);
1311 : return CE_Failure;
1312 : }
1313 : #endif
1314 :
1315 : pabyChunkBuf = static_cast<unsigned char *>(
1316 1 : VSI_MALLOC2_VERBOSE(nPixelSize, static_cast<size_t>(nChunkSize)));
1317 1 : if (pabyChunkBuf == nullptr)
1318 : {
1319 0 : if (bNeedToFreeTransformer)
1320 0 : GDALDestroyTransformer(pTransformArg);
1321 0 : return CE_Failure;
1322 : }
1323 :
1324 1 : OGREnvelope sRasterEnvelope;
1325 1 : sRasterEnvelope.MinX = 0;
1326 1 : sRasterEnvelope.MinY = 0;
1327 1 : sRasterEnvelope.MaxX = poDS->GetRasterXSize();
1328 1 : sRasterEnvelope.MaxY = poDS->GetRasterYSize();
1329 :
1330 : /* --------------------------------------------------------------------
1331 : */
1332 : /* loop over the vectorial geometries */
1333 : /* --------------------------------------------------------------------
1334 : */
1335 1 : pfnProgress(0.0, nullptr, pProgressArg);
1336 3 : for (int iShape = 0; iShape < nGeomCount; iShape++)
1337 : {
1338 :
1339 : const OGRGeometry *poGeometry =
1340 2 : OGRGeometry::FromHandle(pahGeometries[iShape]);
1341 2 : if (poGeometry == nullptr || poGeometry->IsEmpty())
1342 0 : continue;
1343 : /* ------------------------------------------------------------ */
1344 : /* get the envelope of the geometry and transform it to */
1345 : /* pixels coordinates. */
1346 : /* ------------------------------------------------------------ */
1347 2 : OGREnvelope sGeomEnvelope;
1348 2 : poGeometry->getEnvelope(&sGeomEnvelope);
1349 2 : if (pfnTransformer != nullptr)
1350 : {
1351 2 : int anSuccessTransform[2] = {0};
1352 : double apCorners[4];
1353 2 : apCorners[0] = sGeomEnvelope.MinX;
1354 2 : apCorners[1] = sGeomEnvelope.MaxX;
1355 2 : apCorners[2] = sGeomEnvelope.MinY;
1356 2 : apCorners[3] = sGeomEnvelope.MaxY;
1357 :
1358 2 : if (!pfnTransformer(pTransformArg, FALSE, 2, &(apCorners[0]),
1359 : &(apCorners[2]), nullptr,
1360 2 : anSuccessTransform) ||
1361 2 : !anSuccessTransform[0] || !anSuccessTransform[1])
1362 : {
1363 0 : continue;
1364 : }
1365 2 : sGeomEnvelope.MinX = std::min(apCorners[0], apCorners[1]);
1366 2 : sGeomEnvelope.MaxX = std::max(apCorners[0], apCorners[1]);
1367 2 : sGeomEnvelope.MinY = std::min(apCorners[2], apCorners[3]);
1368 2 : sGeomEnvelope.MaxY = std::max(apCorners[2], apCorners[3]);
1369 : }
1370 2 : if (!sGeomEnvelope.Intersects(sRasterEnvelope))
1371 0 : continue;
1372 2 : sGeomEnvelope.Intersect(sRasterEnvelope);
1373 2 : CPLAssert(sGeomEnvelope.MinX >= 0 &&
1374 : sGeomEnvelope.MinX <= poDS->GetRasterXSize());
1375 2 : CPLAssert(sGeomEnvelope.MinY >= 0 &&
1376 : sGeomEnvelope.MinY <= poDS->GetRasterYSize());
1377 2 : CPLAssert(sGeomEnvelope.MaxX >= 0 &&
1378 : sGeomEnvelope.MaxX <= poDS->GetRasterXSize());
1379 2 : CPLAssert(sGeomEnvelope.MaxY >= 0 &&
1380 : sGeomEnvelope.MaxY <= poDS->GetRasterYSize());
1381 2 : const int minBlockX = int(sGeomEnvelope.MinX) / nXBlockSize;
1382 2 : const int minBlockY = int(sGeomEnvelope.MinY) / nYBlockSize;
1383 2 : const int maxBlockX = int(sGeomEnvelope.MaxX + 1) / nXBlockSize;
1384 2 : const int maxBlockY = int(sGeomEnvelope.MaxY + 1) / nYBlockSize;
1385 :
1386 : /* ------------------------------------------------------------ */
1387 : /* loop over the blocks concerned by the geometry */
1388 : /* (by packs of nbBlocksX x nbBlocksY) */
1389 : /* ------------------------------------------------------------ */
1390 :
1391 5 : for (int xB = minBlockX; xB <= maxBlockX; xB += nbBlocksX)
1392 : {
1393 6 : for (int yB = minBlockY; yB <= maxBlockY; yB += nbBlocksY)
1394 : {
1395 :
1396 : /* --------------------------------------------------------------------
1397 : */
1398 : /* ensure to stay in the image */
1399 : /* --------------------------------------------------------------------
1400 : */
1401 3 : int remSBX = std::min(maxBlockX - xB + 1, nbBlocksX);
1402 3 : int remSBY = std::min(maxBlockY - yB + 1, nbBlocksY);
1403 3 : int nThisXChunkSize = nXBlockSize * remSBX;
1404 3 : int nThisYChunkSize = nYBlockSize * remSBY;
1405 6 : if (xB * nXBlockSize + nThisXChunkSize >
1406 3 : poDS->GetRasterXSize())
1407 1 : nThisXChunkSize =
1408 1 : poDS->GetRasterXSize() - xB * nXBlockSize;
1409 6 : if (yB * nYBlockSize + nThisYChunkSize >
1410 3 : poDS->GetRasterYSize())
1411 2 : nThisYChunkSize =
1412 2 : poDS->GetRasterYSize() - yB * nYBlockSize;
1413 :
1414 : /* --------------------------------------------------------------------
1415 : */
1416 : /* read image / process buffer / write buffer */
1417 : /* --------------------------------------------------------------------
1418 : */
1419 3 : eErr = poDS->RasterIO(
1420 : GF_Read, xB * nXBlockSize, yB * nYBlockSize,
1421 : nThisXChunkSize, nThisYChunkSize, pabyChunkBuf,
1422 : nThisXChunkSize, nThisYChunkSize, eType, nBandCount,
1423 : panBandList, 0, 0, 0, nullptr);
1424 3 : if (eErr != CE_None)
1425 0 : break;
1426 :
1427 6 : gv_rasterize_one_shape(
1428 : pabyChunkBuf, xB * nXBlockSize, yB * nYBlockSize,
1429 : nThisXChunkSize, nThisYChunkSize, nBandCount, eType, 0,
1430 : 0, 0, bAllTouched,
1431 3 : OGRGeometry::FromHandle(pahGeometries[iShape]),
1432 : eBurnValueType,
1433 : padfGeomBurnValues
1434 3 : ? padfGeomBurnValues +
1435 3 : static_cast<size_t>(iShape) * nBandCount
1436 : : nullptr,
1437 : panGeomBurnValues
1438 0 : ? panGeomBurnValues +
1439 0 : static_cast<size_t>(iShape) * nBandCount
1440 : : nullptr,
1441 : eBurnValueSource, eMergeAlg, pfnTransformer,
1442 : pTransformArg);
1443 :
1444 3 : eErr = poDS->RasterIO(
1445 : GF_Write, xB * nXBlockSize, yB * nYBlockSize,
1446 : nThisXChunkSize, nThisYChunkSize, pabyChunkBuf,
1447 : nThisXChunkSize, nThisYChunkSize, eType, nBandCount,
1448 : panBandList, 0, 0, 0, nullptr);
1449 3 : if (eErr != CE_None)
1450 0 : break;
1451 : }
1452 : }
1453 :
1454 2 : if (!pfnProgress(iShape / static_cast<double>(nGeomCount), "",
1455 : pProgressArg))
1456 : {
1457 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1458 0 : eErr = CE_Failure;
1459 : }
1460 : }
1461 :
1462 1 : if (!pfnProgress(1., "", pProgressArg))
1463 : {
1464 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1465 0 : eErr = CE_Failure;
1466 : }
1467 : }
1468 :
1469 : /* -------------------------------------------------------------------- */
1470 : /* cleanup */
1471 : /* -------------------------------------------------------------------- */
1472 303 : VSIFree(pabyChunkBuf);
1473 :
1474 303 : if (bNeedToFreeTransformer)
1475 271 : GDALDestroyTransformer(pTransformArg);
1476 :
1477 303 : return eErr;
1478 : }
1479 :
1480 : /************************************************************************/
1481 : /* GDALRasterizeLayers() */
1482 : /************************************************************************/
1483 :
1484 : /**
1485 : * Burn geometries from the specified list of layers into raster.
1486 : *
1487 : * Rasterize all the geometric objects from a list of layers into a raster
1488 : * dataset. The layers are passed as an array of OGRLayerH handlers.
1489 : *
1490 : * If the geometries are in the georeferenced coordinates of the raster
1491 : * dataset, then the pfnTransform may be passed in NULL and one will be
1492 : * derived internally from the geotransform of the dataset. The transform
1493 : * needs to transform the geometry locations into pixel/line coordinates
1494 : * on the raster dataset.
1495 : *
1496 : * The output raster may be of any GDAL supported datatype. An explicit list
1497 : * of burn values for each layer for each band must be passed in.
1498 : *
1499 : * @param hDS output data, must be opened in update mode.
1500 : * @param nBandCount the number of bands to be updated.
1501 : * @param panBandList the list of bands to be updated.
1502 : * @param nLayerCount the number of layers being passed in pahLayers array.
1503 : * @param pahLayers the array of layers to burn in.
1504 : * @param pfnTransformer transformation to apply to geometries to put into
1505 : * pixel/line coordinates on raster. If NULL a geotransform based one will
1506 : * be created internally.
1507 : * @param pTransformArg callback data for transformer.
1508 : * @param padfLayerBurnValues the array of values to burn into the raster.
1509 : * There should be nBandCount values for each layer.
1510 : * @param papszOptions special options controlling rasterization:
1511 : * <ul>
1512 : * <li>"ATTRIBUTE": Identifies an attribute field on the features to be
1513 : * used for a burn in value. The value will be burned into all output
1514 : * bands. If specified, padfLayerBurnValues will not be used and can be a NULL
1515 : * pointer.</li>
1516 : * <li>"CHUNKYSIZE": The height in lines of the chunk to operate on.
1517 : * The larger the chunk size the less times we need to make a pass through all
1518 : * the shapes. If it is not set or set to zero the default chunk size will be
1519 : * used. Default size will be estimated based on the GDAL cache buffer size
1520 : * using formula: cache_size_bytes/scanline_size_bytes, so the chunk will
1521 : * not exceed the cache.</li>
1522 : * <li>"ALL_TOUCHED": May be set to TRUE to set all pixels touched
1523 : * by the line or polygons, not just those whose center is within the polygon
1524 : * (behavior is unspecified when the polygon is just touching the pixel center)
1525 : * or that are selected by Brezenham's line algorithm. Defaults to FALSE.
1526 : * <li>"BURN_VALUE_FROM": May be set to "Z" to use the Z values of the</li>
1527 : * geometries. The value from padfLayerBurnValues or the attribute field value
1528 : * is added to this before burning. In default case dfBurnValue is burned as it
1529 : * is. This is implemented properly only for points and lines for now. Polygons
1530 : * will be burned using the Z value from the first point. The M value may be
1531 : * supported in the future.</li>
1532 : * <li>"MERGE_ALG": May be REPLACE (the default) or ADD. REPLACE results in
1533 : * overwriting of value, while ADD adds the new value to the existing raster,
1534 : * suitable for heatmaps for instance.</li>
1535 : * </ul>
1536 : * @param pfnProgress the progress function to report completion.
1537 : * @param pProgressArg callback data for progress function.
1538 : *
1539 : * @return CE_None on success or CE_Failure on error.
1540 : */
1541 :
1542 48 : CPLErr GDALRasterizeLayers(GDALDatasetH hDS, int nBandCount, int *panBandList,
1543 : int nLayerCount, OGRLayerH *pahLayers,
1544 : GDALTransformerFunc pfnTransformer,
1545 : void *pTransformArg, double *padfLayerBurnValues,
1546 : char **papszOptions, GDALProgressFunc pfnProgress,
1547 : void *pProgressArg)
1548 :
1549 : {
1550 48 : VALIDATE_POINTER1(hDS, "GDALRasterizeLayers", CE_Failure);
1551 :
1552 48 : if (pfnProgress == nullptr)
1553 48 : pfnProgress = GDALDummyProgress;
1554 :
1555 : /* -------------------------------------------------------------------- */
1556 : /* Do some rudimentary arg checking. */
1557 : /* -------------------------------------------------------------------- */
1558 48 : if (nBandCount == 0 || nLayerCount == 0)
1559 0 : return CE_None;
1560 :
1561 48 : GDALDataset *poDS = GDALDataset::FromHandle(hDS);
1562 :
1563 : // Prototype band.
1564 48 : GDALRasterBand *poBand = poDS->GetRasterBand(panBandList[0]);
1565 48 : if (poBand == nullptr)
1566 0 : return CE_Failure;
1567 :
1568 : /* -------------------------------------------------------------------- */
1569 : /* Options */
1570 : /* -------------------------------------------------------------------- */
1571 48 : int bAllTouched = FALSE;
1572 48 : GDALBurnValueSrc eBurnValueSource = GBV_UserBurnValue;
1573 48 : GDALRasterMergeAlg eMergeAlg = GRMA_Replace;
1574 48 : if (GDALRasterizeOptions(papszOptions, &bAllTouched, &eBurnValueSource,
1575 48 : &eMergeAlg, nullptr) == CE_Failure)
1576 : {
1577 0 : return CE_Failure;
1578 : }
1579 :
1580 : /* -------------------------------------------------------------------- */
1581 : /* Establish a chunksize to operate on. The larger the chunk */
1582 : /* size the less times we need to make a pass through all the */
1583 : /* shapes. */
1584 : /* -------------------------------------------------------------------- */
1585 48 : const char *pszYChunkSize = CSLFetchNameValue(papszOptions, "CHUNKYSIZE");
1586 :
1587 48 : const GDALDataType eType = poBand->GetRasterDataType();
1588 :
1589 : const int nScanlineBytes =
1590 48 : nBandCount * poDS->GetRasterXSize() * GDALGetDataTypeSizeBytes(eType);
1591 :
1592 48 : int nYChunkSize = 0;
1593 48 : if (!(pszYChunkSize && ((nYChunkSize = atoi(pszYChunkSize))) != 0))
1594 : {
1595 48 : const GIntBig nYChunkSize64 = GDALGetCacheMax64() / nScanlineBytes;
1596 48 : nYChunkSize = static_cast<int>(
1597 48 : std::min<GIntBig>(nYChunkSize64, std::numeric_limits<int>::max()));
1598 : }
1599 :
1600 48 : if (nYChunkSize < 1)
1601 0 : nYChunkSize = 1;
1602 48 : if (nYChunkSize > poDS->GetRasterYSize())
1603 48 : nYChunkSize = poDS->GetRasterYSize();
1604 :
1605 48 : CPLDebug("GDAL", "Rasterizer operating on %d swaths of %d scanlines.",
1606 48 : DIV_ROUND_UP(poDS->GetRasterYSize(), nYChunkSize), nYChunkSize);
1607 : unsigned char *pabyChunkBuf = static_cast<unsigned char *>(
1608 48 : VSI_MALLOC2_VERBOSE(nYChunkSize, nScanlineBytes));
1609 48 : if (pabyChunkBuf == nullptr)
1610 : {
1611 0 : return CE_Failure;
1612 : }
1613 :
1614 : /* -------------------------------------------------------------------- */
1615 : /* Read the image once for all layers if user requested to render */
1616 : /* the whole raster in single chunk. */
1617 : /* -------------------------------------------------------------------- */
1618 48 : if (nYChunkSize == poDS->GetRasterYSize())
1619 : {
1620 48 : if (poDS->RasterIO(GF_Read, 0, 0, poDS->GetRasterXSize(), nYChunkSize,
1621 : pabyChunkBuf, poDS->GetRasterXSize(), nYChunkSize,
1622 : eType, nBandCount, panBandList, 0, 0, 0,
1623 48 : nullptr) != CE_None)
1624 : {
1625 0 : CPLFree(pabyChunkBuf);
1626 0 : return CE_Failure;
1627 : }
1628 : }
1629 :
1630 : /* ==================================================================== */
1631 : /* Read the specified layers transforming and rasterizing */
1632 : /* geometries. */
1633 : /* ==================================================================== */
1634 48 : CPLErr eErr = CE_None;
1635 48 : const char *pszBurnAttribute = CSLFetchNameValue(papszOptions, "ATTRIBUTE");
1636 :
1637 48 : pfnProgress(0.0, nullptr, pProgressArg);
1638 :
1639 96 : for (int iLayer = 0; iLayer < nLayerCount; iLayer++)
1640 : {
1641 48 : OGRLayer *poLayer = reinterpret_cast<OGRLayer *>(pahLayers[iLayer]);
1642 :
1643 48 : if (!poLayer)
1644 : {
1645 0 : CPLError(CE_Warning, CPLE_AppDefined,
1646 : "Layer element number %d is NULL, skipping.", iLayer);
1647 0 : continue;
1648 : }
1649 :
1650 : /* --------------------------------------------------------------------
1651 : */
1652 : /* If the layer does not contain any features just skip it. */
1653 : /* Do not force the feature count, so if driver doesn't know */
1654 : /* exact number of features, go down the normal way. */
1655 : /* --------------------------------------------------------------------
1656 : */
1657 48 : if (poLayer->GetFeatureCount(FALSE) == 0)
1658 0 : continue;
1659 :
1660 48 : int iBurnField = -1;
1661 48 : double *padfBurnValues = nullptr;
1662 :
1663 48 : if (pszBurnAttribute)
1664 : {
1665 : iBurnField =
1666 1 : poLayer->GetLayerDefn()->GetFieldIndex(pszBurnAttribute);
1667 1 : if (iBurnField == -1)
1668 : {
1669 0 : CPLError(CE_Warning, CPLE_AppDefined,
1670 : "Failed to find field %s on layer %s, skipping.",
1671 0 : pszBurnAttribute, poLayer->GetLayerDefn()->GetName());
1672 0 : continue;
1673 : }
1674 : }
1675 : else
1676 : {
1677 47 : padfBurnValues = padfLayerBurnValues + iLayer * nBandCount;
1678 : }
1679 :
1680 : /* --------------------------------------------------------------------
1681 : */
1682 : /* If we have no transformer, create the one from input file */
1683 : /* projection. Note that each layer can be georefernced */
1684 : /* separately. */
1685 : /* --------------------------------------------------------------------
1686 : */
1687 48 : bool bNeedToFreeTransformer = false;
1688 :
1689 48 : if (pfnTransformer == nullptr)
1690 : {
1691 48 : char *pszProjection = nullptr;
1692 48 : bNeedToFreeTransformer = true;
1693 :
1694 48 : const OGRSpatialReference *poSRS = poLayer->GetSpatialRef();
1695 48 : if (!poSRS)
1696 : {
1697 9 : if (poDS->GetSpatialRef() != nullptr ||
1698 16 : poDS->GetGCPSpatialRef() != nullptr ||
1699 7 : poDS->GetMetadata("RPC") != nullptr)
1700 : {
1701 2 : CPLError(
1702 : CE_Warning, CPLE_AppDefined,
1703 : "Failed to fetch spatial reference on layer %s "
1704 : "to build transformer, assuming matching coordinate "
1705 : "systems.",
1706 2 : poLayer->GetLayerDefn()->GetName());
1707 : }
1708 : }
1709 : else
1710 : {
1711 39 : poSRS->exportToWkt(&pszProjection);
1712 : }
1713 :
1714 48 : CPLStringList aosTransformerOptions;
1715 48 : if (pszProjection != nullptr)
1716 39 : aosTransformerOptions.SetNameValue("SRC_SRS", pszProjection);
1717 48 : GDALGeoTransform gt;
1718 48 : if (poDS->GetGeoTransform(gt) != CE_None &&
1719 48 : poDS->GetGCPCount() == 0 && poDS->GetMetadata("RPC") == nullptr)
1720 : {
1721 : aosTransformerOptions.SetNameValue("DST_METHOD",
1722 0 : "NO_GEOTRANSFORM");
1723 : }
1724 :
1725 48 : pTransformArg = GDALCreateGenImgProjTransformer2(
1726 48 : nullptr, hDS, aosTransformerOptions.List());
1727 48 : pfnTransformer = GDALGenImgProjTransform;
1728 :
1729 48 : CPLFree(pszProjection);
1730 48 : if (pTransformArg == nullptr)
1731 : {
1732 0 : CPLFree(pabyChunkBuf);
1733 0 : return CE_Failure;
1734 : }
1735 : }
1736 :
1737 48 : poLayer->ResetReading();
1738 :
1739 : /* --------------------------------------------------------------------
1740 : */
1741 : /* Loop over image in designated chunks. */
1742 : /* --------------------------------------------------------------------
1743 : */
1744 :
1745 : double *padfAttrValues = static_cast<double *>(
1746 48 : VSI_MALLOC_VERBOSE(sizeof(double) * nBandCount));
1747 48 : if (padfAttrValues == nullptr)
1748 0 : eErr = CE_Failure;
1749 :
1750 96 : for (int iY = 0; iY < poDS->GetRasterYSize() && eErr == CE_None;
1751 48 : iY += nYChunkSize)
1752 : {
1753 48 : int nThisYChunkSize = nYChunkSize;
1754 48 : if (nThisYChunkSize + iY > poDS->GetRasterYSize())
1755 0 : nThisYChunkSize = poDS->GetRasterYSize() - iY;
1756 :
1757 : // Only re-read image if not a single chunk is being rendered.
1758 48 : if (nYChunkSize < poDS->GetRasterYSize())
1759 : {
1760 0 : eErr = poDS->RasterIO(
1761 : GF_Read, 0, iY, poDS->GetRasterXSize(), nThisYChunkSize,
1762 : pabyChunkBuf, poDS->GetRasterXSize(), nThisYChunkSize,
1763 : eType, nBandCount, panBandList, 0, 0, 0, nullptr);
1764 0 : if (eErr != CE_None)
1765 0 : break;
1766 : }
1767 :
1768 112 : for (auto &poFeat : poLayer)
1769 : {
1770 64 : OGRGeometry *poGeom = poFeat->GetGeometryRef();
1771 :
1772 64 : if (pszBurnAttribute)
1773 : {
1774 : const double dfAttrValue =
1775 5 : poFeat->GetFieldAsDouble(iBurnField);
1776 20 : for (int iBand = 0; iBand < nBandCount; iBand++)
1777 15 : padfAttrValues[iBand] = dfAttrValue;
1778 :
1779 5 : padfBurnValues = padfAttrValues;
1780 : }
1781 :
1782 64 : gv_rasterize_one_shape(
1783 : pabyChunkBuf, 0, iY, poDS->GetRasterXSize(),
1784 : nThisYChunkSize, nBandCount, eType, 0, 0, 0, bAllTouched,
1785 : poGeom, GDT_Float64, padfBurnValues, nullptr,
1786 : eBurnValueSource, eMergeAlg, pfnTransformer, pTransformArg);
1787 : }
1788 :
1789 : // Only write image if not a single chunk is being rendered.
1790 48 : if (nYChunkSize < poDS->GetRasterYSize())
1791 : {
1792 0 : eErr = poDS->RasterIO(
1793 : GF_Write, 0, iY, poDS->GetRasterXSize(), nThisYChunkSize,
1794 : pabyChunkBuf, poDS->GetRasterXSize(), nThisYChunkSize,
1795 : eType, nBandCount, panBandList, 0, 0, 0, nullptr);
1796 : }
1797 :
1798 48 : poLayer->ResetReading();
1799 :
1800 48 : if (!pfnProgress((iY + nThisYChunkSize) /
1801 48 : static_cast<double>(poDS->GetRasterYSize()),
1802 : "", pProgressArg))
1803 : {
1804 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1805 0 : eErr = CE_Failure;
1806 : }
1807 : }
1808 :
1809 48 : VSIFree(padfAttrValues);
1810 :
1811 48 : if (bNeedToFreeTransformer)
1812 : {
1813 48 : GDALDestroyTransformer(pTransformArg);
1814 48 : pTransformArg = nullptr;
1815 48 : pfnTransformer = nullptr;
1816 : }
1817 : }
1818 :
1819 : /* -------------------------------------------------------------------- */
1820 : /* Write out the image once for all layers if user requested */
1821 : /* to render the whole raster in single chunk. */
1822 : /* -------------------------------------------------------------------- */
1823 48 : if (eErr == CE_None && nYChunkSize == poDS->GetRasterYSize())
1824 : {
1825 : eErr =
1826 48 : poDS->RasterIO(GF_Write, 0, 0, poDS->GetRasterXSize(), nYChunkSize,
1827 : pabyChunkBuf, poDS->GetRasterXSize(), nYChunkSize,
1828 : eType, nBandCount, panBandList, 0, 0, 0, nullptr);
1829 : }
1830 :
1831 : /* -------------------------------------------------------------------- */
1832 : /* cleanup */
1833 : /* -------------------------------------------------------------------- */
1834 48 : VSIFree(pabyChunkBuf);
1835 :
1836 48 : return eErr;
1837 : }
1838 :
1839 : /************************************************************************/
1840 : /* GDALRasterizeLayersBuf() */
1841 : /************************************************************************/
1842 :
1843 : /**
1844 : * Burn geometries from the specified list of layer into raster.
1845 : *
1846 : * Rasterize all the geometric objects from a list of layers into supplied
1847 : * raster buffer. The layers are passed as an array of OGRLayerH handlers.
1848 : *
1849 : * If the geometries are in the georeferenced coordinates of the raster
1850 : * dataset, then the pfnTransform may be passed in NULL and one will be
1851 : * derived internally from the geotransform of the dataset. The transform
1852 : * needs to transform the geometry locations into pixel/line coordinates
1853 : * of the target raster.
1854 : *
1855 : * The output raster may be of any GDAL supported datatype(non complex).
1856 : *
1857 : * @param pData pointer to the output data array.
1858 : *
1859 : * @param nBufXSize width of the output data array in pixels.
1860 : *
1861 : * @param nBufYSize height of the output data array in pixels.
1862 : *
1863 : * @param eBufType data type of the output data array.
1864 : *
1865 : * @param nPixelSpace The byte offset from the start of one pixel value in
1866 : * pData to the start of the next pixel value within a scanline. If defaulted
1867 : * (0) the size of the datatype eBufType is used.
1868 : *
1869 : * @param nLineSpace The byte offset from the start of one scanline in
1870 : * pData to the start of the next. If defaulted the size of the datatype
1871 : * eBufType * nBufXSize is used.
1872 : *
1873 : * @param nLayerCount the number of layers being passed in pahLayers array.
1874 : *
1875 : * @param pahLayers the array of layers to burn in.
1876 : *
1877 : * @param pszDstProjection WKT defining the coordinate system of the target
1878 : * raster.
1879 : *
1880 : * @param padfDstGeoTransform geotransformation matrix of the target raster.
1881 : *
1882 : * @param pfnTransformer transformation to apply to geometries to put into
1883 : * pixel/line coordinates on raster. If NULL a geotransform based one will
1884 : * be created internally.
1885 : *
1886 : * @param pTransformArg callback data for transformer.
1887 : *
1888 : * @param dfBurnValue the value to burn into the raster.
1889 : *
1890 : * @param papszOptions special options controlling rasterization:
1891 : * <ul>
1892 : * <li>"ATTRIBUTE": Identifies an attribute field on the features to be
1893 : * used for a burn in value. The value will be burned into all output
1894 : * bands. If specified, padfLayerBurnValues will not be used and can be a NULL
1895 : * pointer.</li>
1896 : * <li>"ALL_TOUCHED": May be set to TRUE to set all pixels touched
1897 : * by the line or polygons, not just those whose center is within the polygon
1898 : * (behavior is unspecified when the polygon is just touching the pixel center)
1899 : * or that are selected by Brezenham's line algorithm. Defaults to FALSE.</li>
1900 : * <li>"BURN_VALUE_FROM": May be set to "Z" to use
1901 : * the Z values of the geometries. dfBurnValue or the attribute field value is
1902 : * added to this before burning. In default case dfBurnValue is burned as it
1903 : * is. This is implemented properly only for points and lines for now. Polygons
1904 : * will be burned using the Z value from the first point. The M value may
1905 : * be supported in the future.</li>
1906 : * <li>"MERGE_ALG": May be REPLACE (the default) or ADD. REPLACE
1907 : * results in overwriting of value, while ADD adds the new value to the
1908 : * existing raster, suitable for heatmaps for instance.</li>
1909 : * </ul>
1910 : *
1911 : * @param pfnProgress the progress function to report completion.
1912 : *
1913 : * @param pProgressArg callback data for progress function.
1914 : *
1915 : *
1916 : * @return CE_None on success or CE_Failure on error.
1917 : */
1918 :
1919 0 : CPLErr GDALRasterizeLayersBuf(
1920 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
1921 : int nPixelSpace, int nLineSpace, int nLayerCount, OGRLayerH *pahLayers,
1922 : const char *pszDstProjection, double *padfDstGeoTransform,
1923 : GDALTransformerFunc pfnTransformer, void *pTransformArg, double dfBurnValue,
1924 : char **papszOptions, GDALProgressFunc pfnProgress, void *pProgressArg)
1925 :
1926 : {
1927 : /* -------------------------------------------------------------------- */
1928 : /* check eType, Avoid not supporting data types */
1929 : /* -------------------------------------------------------------------- */
1930 0 : if (GDALDataTypeIsComplex(eBufType) || eBufType <= GDT_Unknown ||
1931 0 : eBufType >= GDT_TypeCount)
1932 : {
1933 0 : CPLError(CE_Failure, CPLE_NotSupported,
1934 : "GDALRasterizeLayersBuf(): unsupported data type of eBufType");
1935 0 : return CE_Failure;
1936 : }
1937 :
1938 : /* -------------------------------------------------------------------- */
1939 : /* If pixel and line spaceing are defaulted assign reasonable */
1940 : /* value assuming a packed buffer. */
1941 : /* -------------------------------------------------------------------- */
1942 0 : int nTypeSizeBytes = GDALGetDataTypeSizeBytes(eBufType);
1943 0 : if (nPixelSpace == 0)
1944 : {
1945 0 : nPixelSpace = nTypeSizeBytes;
1946 : }
1947 0 : if (nPixelSpace < nTypeSizeBytes)
1948 : {
1949 0 : CPLError(CE_Failure, CPLE_NotSupported,
1950 : "GDALRasterizeLayersBuf(): unsupported value of nPixelSpace");
1951 0 : return CE_Failure;
1952 : }
1953 :
1954 0 : if (nLineSpace == 0)
1955 : {
1956 0 : nLineSpace = nPixelSpace * nBufXSize;
1957 : }
1958 0 : if (nLineSpace < nPixelSpace * nBufXSize)
1959 : {
1960 0 : CPLError(CE_Failure, CPLE_NotSupported,
1961 : "GDALRasterizeLayersBuf(): unsupported value of nLineSpace");
1962 0 : return CE_Failure;
1963 : }
1964 :
1965 0 : if (pfnProgress == nullptr)
1966 0 : pfnProgress = GDALDummyProgress;
1967 :
1968 : /* -------------------------------------------------------------------- */
1969 : /* Do some rudimentary arg checking. */
1970 : /* -------------------------------------------------------------------- */
1971 0 : if (nLayerCount == 0)
1972 0 : return CE_None;
1973 :
1974 : /* -------------------------------------------------------------------- */
1975 : /* Options */
1976 : /* -------------------------------------------------------------------- */
1977 0 : int bAllTouched = FALSE;
1978 0 : GDALBurnValueSrc eBurnValueSource = GBV_UserBurnValue;
1979 0 : GDALRasterMergeAlg eMergeAlg = GRMA_Replace;
1980 0 : if (GDALRasterizeOptions(papszOptions, &bAllTouched, &eBurnValueSource,
1981 0 : &eMergeAlg, nullptr) == CE_Failure)
1982 : {
1983 0 : return CE_Failure;
1984 : }
1985 :
1986 : /* ==================================================================== */
1987 : /* Read the specified layers transforming and rasterizing */
1988 : /* geometries. */
1989 : /* ==================================================================== */
1990 0 : CPLErr eErr = CE_None;
1991 0 : const char *pszBurnAttribute = CSLFetchNameValue(papszOptions, "ATTRIBUTE");
1992 :
1993 0 : pfnProgress(0.0, nullptr, pProgressArg);
1994 :
1995 0 : for (int iLayer = 0; iLayer < nLayerCount; iLayer++)
1996 : {
1997 0 : OGRLayer *poLayer = reinterpret_cast<OGRLayer *>(pahLayers[iLayer]);
1998 :
1999 0 : if (!poLayer)
2000 : {
2001 0 : CPLError(CE_Warning, CPLE_AppDefined,
2002 : "Layer element number %d is NULL, skipping.", iLayer);
2003 0 : continue;
2004 : }
2005 :
2006 : /* --------------------------------------------------------------------
2007 : */
2008 : /* If the layer does not contain any features just skip it. */
2009 : /* Do not force the feature count, so if driver doesn't know */
2010 : /* exact number of features, go down the normal way. */
2011 : /* --------------------------------------------------------------------
2012 : */
2013 0 : if (poLayer->GetFeatureCount(FALSE) == 0)
2014 0 : continue;
2015 :
2016 0 : int iBurnField = -1;
2017 0 : if (pszBurnAttribute)
2018 : {
2019 : iBurnField =
2020 0 : poLayer->GetLayerDefn()->GetFieldIndex(pszBurnAttribute);
2021 0 : if (iBurnField == -1)
2022 : {
2023 0 : CPLError(CE_Warning, CPLE_AppDefined,
2024 : "Failed to find field %s on layer %s, skipping.",
2025 0 : pszBurnAttribute, poLayer->GetLayerDefn()->GetName());
2026 0 : continue;
2027 : }
2028 : }
2029 :
2030 : /* --------------------------------------------------------------------
2031 : */
2032 : /* If we have no transformer, create the one from input file */
2033 : /* projection. Note that each layer can be georefernced */
2034 : /* separately. */
2035 : /* --------------------------------------------------------------------
2036 : */
2037 0 : bool bNeedToFreeTransformer = false;
2038 :
2039 0 : if (pfnTransformer == nullptr)
2040 : {
2041 0 : char *pszProjection = nullptr;
2042 0 : bNeedToFreeTransformer = true;
2043 :
2044 0 : const OGRSpatialReference *poSRS = poLayer->GetSpatialRef();
2045 0 : if (!poSRS)
2046 : {
2047 0 : CPLError(CE_Warning, CPLE_AppDefined,
2048 : "Failed to fetch spatial reference on layer %s "
2049 : "to build transformer, assuming matching coordinate "
2050 : "systems.",
2051 0 : poLayer->GetLayerDefn()->GetName());
2052 : }
2053 : else
2054 : {
2055 0 : poSRS->exportToWkt(&pszProjection);
2056 : }
2057 :
2058 0 : pTransformArg = GDALCreateGenImgProjTransformer3(
2059 : pszProjection, nullptr, pszDstProjection, padfDstGeoTransform);
2060 0 : pfnTransformer = GDALGenImgProjTransform;
2061 :
2062 0 : CPLFree(pszProjection);
2063 : }
2064 :
2065 0 : for (auto &poFeat : poLayer)
2066 : {
2067 0 : OGRGeometry *poGeom = poFeat->GetGeometryRef();
2068 :
2069 0 : if (pszBurnAttribute)
2070 0 : dfBurnValue = poFeat->GetFieldAsDouble(iBurnField);
2071 :
2072 0 : gv_rasterize_one_shape(
2073 : static_cast<unsigned char *>(pData), 0, 0, nBufXSize, nBufYSize,
2074 : 1, eBufType, nPixelSpace, nLineSpace, 0, bAllTouched, poGeom,
2075 : GDT_Float64, &dfBurnValue, nullptr, eBurnValueSource, eMergeAlg,
2076 : pfnTransformer, pTransformArg);
2077 : }
2078 :
2079 0 : poLayer->ResetReading();
2080 :
2081 0 : if (!pfnProgress(1, "", pProgressArg))
2082 : {
2083 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
2084 0 : eErr = CE_Failure;
2085 : }
2086 :
2087 0 : if (bNeedToFreeTransformer)
2088 : {
2089 0 : GDALDestroyTransformer(pTransformArg);
2090 0 : pTransformArg = nullptr;
2091 0 : pfnTransformer = nullptr;
2092 : }
2093 : }
2094 :
2095 0 : return eErr;
2096 : }
|
