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