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