Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Vector rasterization.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "gdal_alg.h"
16 : #include "gdal_alg_priv.h"
17 :
18 : #include <climits>
19 : #include <cstddef>
20 : #include <cstdlib>
21 : #include <cstring>
22 : #include <cfloat>
23 : #include <limits>
24 : #include <vector>
25 : #include <algorithm>
26 :
27 : #include "cpl_conv.h"
28 : #include "cpl_error.h"
29 : #include "cpl_progress.h"
30 : #include "cpl_string.h"
31 : #include "cpl_vsi.h"
32 : #include "gdal.h"
33 : #include "gdal_priv.h"
34 : #include "gdal_priv_templates.hpp"
35 : #include "ogr_api.h"
36 : #include "ogr_core.h"
37 : #include "ogr_feature.h"
38 : #include "ogr_geometry.h"
39 : #include "ogr_spatialref.h"
40 : #include "ogrsf_frmts.h"
41 :
42 0 : template <typename T> static inline T SaturatedAddSigned(T a, T b)
43 : {
44 0 : if (a > 0 && b > 0 && a > std::numeric_limits<T>::max() - b)
45 : {
46 0 : return std::numeric_limits<T>::max();
47 : }
48 0 : else if (a < 0 && b < 0 && a < std::numeric_limits<T>::min() - b)
49 : {
50 0 : return std::numeric_limits<T>::min();
51 : }
52 : else
53 : {
54 0 : return a + b;
55 : }
56 : }
57 :
58 : /************************************************************************/
59 : /* MakeKey() */
60 : /************************************************************************/
61 :
62 818 : inline uint64_t MakeKey(int y, int x)
63 : {
64 818 : return (static_cast<uint64_t>(y) << 32) | static_cast<uint64_t>(x);
65 : }
66 :
67 : /************************************************************************/
68 : /* gvBurnScanlineBasic() */
69 : /************************************************************************/
70 : template <typename T>
71 10416 : static inline void gvBurnScanlineBasic(GDALRasterizeInfo *psInfo, int nY,
72 : int nXStart, int nXEnd, double dfVariant)
73 :
74 : {
75 21726 : for (int iBand = 0; iBand < psInfo->nBands; iBand++)
76 : {
77 11310 : const double burnValue =
78 11310 : (psInfo->burnValues.double_values[iBand] +
79 11310 : ((psInfo->eBurnValueSource == GBV_UserBurnValue) ? 0 : dfVariant));
80 :
81 11310 : unsigned char *pabyInsert =
82 11310 : psInfo->pabyChunkBuf + iBand * psInfo->nBandSpace +
83 11310 : nY * psInfo->nLineSpace + nXStart * psInfo->nPixelSpace;
84 11310 : if (psInfo->eMergeAlg == GRMA_Add)
85 : {
86 280 : if (psInfo->poSetVisitedPoints)
87 : {
88 280 : CPLAssert(!psInfo->bFillSetVisitedPoints);
89 280 : uint64_t nKey = MakeKey(nY, nXStart);
90 280 : auto &oSetVisitedPoints = *(psInfo->poSetVisitedPoints);
91 3791 : for (int nX = nXStart; nX <= nXEnd; ++nX)
92 : {
93 3511 : if (oSetVisitedPoints.find(nKey) == oSetVisitedPoints.end())
94 : {
95 3379 : double dfVal = static_cast<double>(
96 3379 : *reinterpret_cast<T *>(pabyInsert)) +
97 : burnValue;
98 3379 : GDALCopyWord(dfVal, *reinterpret_cast<T *>(pabyInsert));
99 : }
100 3511 : pabyInsert += psInfo->nPixelSpace;
101 3511 : ++nKey;
102 : }
103 : }
104 : else
105 : {
106 0 : for (int nX = nXStart; nX <= nXEnd; ++nX)
107 : {
108 0 : double dfVal = static_cast<double>(
109 0 : *reinterpret_cast<T *>(pabyInsert)) +
110 : burnValue;
111 0 : GDALCopyWord(dfVal, *reinterpret_cast<T *>(pabyInsert));
112 0 : pabyInsert += psInfo->nPixelSpace;
113 : }
114 : }
115 : }
116 : else
117 : {
118 : T nVal;
119 11030 : GDALCopyWord(burnValue, nVal);
120 11645094 : for (int nX = nXStart; nX <= nXEnd; ++nX)
121 : {
122 11634137 : *reinterpret_cast<T *>(pabyInsert) = nVal;
123 11634137 : pabyInsert += psInfo->nPixelSpace;
124 : }
125 : }
126 : }
127 10416 : }
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 10461 : static void gvBurnScanline(GDALRasterizeInfo *psInfo, int nY, int nXStart,
188 : int nXEnd, double dfVariant)
189 :
190 : {
191 10461 : if (nXStart > nXEnd)
192 42 : return;
193 :
194 10419 : CPLAssert(nY >= 0 && nY < psInfo->nYSize);
195 10419 : CPLAssert(nXStart < psInfo->nXSize);
196 10419 : CPLAssert(nXEnd >= 0);
197 :
198 10419 : if (nXStart < 0)
199 6158 : nXStart = 0;
200 10419 : if (nXEnd >= psInfo->nXSize)
201 125 : nXEnd = psInfo->nXSize - 1;
202 :
203 10419 : 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 10416 : switch (psInfo->eType)
218 : {
219 7859 : case GDT_Byte:
220 7859 : gvBurnScanlineBasic<GByte>(psInfo, nY, nXStart, nXEnd, dfVariant);
221 7859 : 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 2294 : case GDT_Float64:
253 2294 : gvBurnScanlineBasic<double>(psInfo, nY, nXStart, nXEnd, dfVariant);
254 2294 : 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 43282 : static inline void gvBurnPointBasic(GDALRasterizeInfo *psInfo, int nY, int nX,
272 : double dfVariant)
273 :
274 : {
275 90749 : for (int iBand = 0; iBand < psInfo->nBands; iBand++)
276 : {
277 47467 : double burnValue =
278 47467 : (psInfo->burnValues.double_values[iBand] +
279 47467 : ((psInfo->eBurnValueSource == GBV_UserBurnValue) ? 0 : dfVariant));
280 47467 : unsigned char *pbyInsert =
281 47467 : psInfo->pabyChunkBuf + iBand * psInfo->nBandSpace +
282 47467 : nY * psInfo->nLineSpace + nX * psInfo->nPixelSpace;
283 :
284 47467 : T *pbyPixel = reinterpret_cast<T *>(pbyInsert);
285 47467 : if (psInfo->eMergeAlg == GRMA_Add)
286 1020 : burnValue += static_cast<double>(*pbyPixel);
287 47467 : GDALCopyWord(burnValue, *pbyPixel);
288 : }
289 43282 : }
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 43294 : static void gvBurnPoint(GDALRasterizeInfo *psInfo, int nY, int nX,
315 : double dfVariant)
316 :
317 : {
318 :
319 43294 : CPLAssert(nY >= 0 && nY < psInfo->nYSize);
320 43294 : CPLAssert(nX >= 0 && nX < psInfo->nXSize);
321 :
322 43294 : if (psInfo->poSetVisitedPoints)
323 : {
324 538 : const uint64_t nKey = MakeKey(nY, nX);
325 538 : if (psInfo->poSetVisitedPoints->find(nKey) ==
326 1076 : psInfo->poSetVisitedPoints->end())
327 : {
328 526 : if (psInfo->bFillSetVisitedPoints)
329 526 : psInfo->poSetVisitedPoints->insert(nKey);
330 : }
331 : else
332 : {
333 12 : return;
334 : }
335 : }
336 :
337 43282 : 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 43282 : switch (psInfo->eType)
352 : {
353 16272 : case GDT_Byte:
354 16272 : gvBurnPointBasic<GByte>(psInfo, nY, nX, dfVariant);
355 16272 : 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 26278 : case GDT_Float64:
384 26278 : gvBurnPointBasic<double>(psInfo, nY, nX, dfVariant);
385 26278 : 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 2254 : 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 2254 : if (poShape == nullptr || poShape->IsEmpty())
410 0 : return;
411 :
412 2254 : const OGRwkbGeometryType eFlatType = wkbFlatten(poShape->getGeometryType());
413 :
414 2254 : 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 2249 : else if (EQUAL(poShape->getGeometryName(), "LINEARRING"))
435 : {
436 190 : const auto poRing = poShape->toLinearRing();
437 190 : const int nCount = poRing->getNumPoints();
438 190 : const size_t nNewCount = aPointX.size() + static_cast<size_t>(nCount);
439 :
440 190 : aPointX.reserve(nNewCount);
441 190 : aPointY.reserve(nNewCount);
442 190 : if (eBurnValueSrc != GBV_UserBurnValue)
443 7 : aPointVariant.reserve(nNewCount);
444 190 : if (poRing->isClockwise())
445 : {
446 2338 : for (int i = 0; i < nCount; i++)
447 : {
448 2235 : aPointX.push_back(poRing->getX(i));
449 2235 : aPointY.push_back(poRing->getY(i));
450 2235 : if (eBurnValueSrc != GBV_UserBurnValue)
451 : {
452 : /*switch( eBurnValueSrc )
453 : {
454 : case GBV_Z:*/
455 19 : aPointVariant.push_back(poRing->getZ(i));
456 : /*break;
457 : case GBV_M:
458 : aPointVariant.push_back( poRing->getM(i) );
459 : }*/
460 : }
461 : }
462 : }
463 : else
464 : {
465 5545 : for (int i = nCount - 1; i >= 0; i--)
466 : {
467 5458 : aPointX.push_back(poRing->getX(i));
468 5458 : aPointY.push_back(poRing->getY(i));
469 5458 : if (eBurnValueSrc != GBV_UserBurnValue)
470 : {
471 : /*switch( eBurnValueSrc )
472 : {
473 : case GBV_Z:*/
474 15 : aPointVariant.push_back(poRing->getZ(i));
475 : /*break;
476 : case GBV_M:
477 : aPointVariant.push_back( poRing->getM(i) );
478 : }*/
479 : }
480 : }
481 : }
482 190 : aPartSize.push_back(nCount);
483 : }
484 2059 : 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 179 : else if (eFlatType == wkbPolygon)
513 : {
514 179 : const auto poPolygon = poShape->toPolygon();
515 :
516 179 : GDALCollectRingsFromGeometry(poPolygon->getExteriorRing(), aPointX,
517 : aPointY, aPointVariant, aPartSize,
518 : eBurnValueSrc);
519 :
520 190 : for (int i = 0; i < poPolygon->getNumInteriorRings(); i++)
521 11 : GDALCollectRingsFromGeometry(poPolygon->getInteriorRing(i), aPointX,
522 : aPointY, aPointVariant, aPartSize,
523 : eBurnValueSrc);
524 : }
525 0 : else if (eFlatType == wkbMultiPoint || eFlatType == wkbMultiLineString ||
526 0 : eFlatType == wkbMultiPolygon || eFlatType == wkbGeometryCollection)
527 : {
528 0 : const auto poGC = poShape->toGeometryCollection();
529 0 : for (int i = 0; i < poGC->getNumGeometries(); i++)
530 0 : GDALCollectRingsFromGeometry(poGC->getGeometryRef(i), aPointX,
531 : aPointY, aPointVariant, aPartSize,
532 0 : eBurnValueSrc);
533 : }
534 : else
535 : {
536 0 : CPLDebug("GDAL", "Rasterizer ignoring non-polygonal geometry.");
537 : }
538 : }
539 :
540 : /************************************************************************
541 : * gv_rasterize_one_shape()
542 : *
543 : * @param pabyChunkBuf buffer to which values will be burned
544 : * @param nXOff chunk column offset from left edge of raster
545 : * @param nYOff chunk scanline offset from top of raster
546 : * @param nXSize number of columns in chunk
547 : * @param nYSize number of rows in chunk
548 : * @param nBands number of bands in chunk
549 : * @param eType data type of pabyChunkBuf
550 : * @param nPixelSpace number of bytes between adjacent pixels in chunk
551 : * (0 to calculate automatically)
552 : * @param nLineSpace number of bytes between adjacent scanlines in chunk
553 : * (0 to calculate automatically)
554 : * @param nBandSpace number of bytes between adjacent bands in chunk
555 : * (0 to calculate automatically)
556 : * @param bAllTouched burn value to all touched pixels?
557 : * @param poShape geometry to rasterize, in original coordinates
558 : * @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 2091 : 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 2091 : if (poShape == nullptr || poShape->IsEmpty())
580 27 : return;
581 2090 : const auto eGeomType = wkbFlatten(poShape->getGeometryType());
582 :
583 2090 : 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 2064 : if (nPixelSpace == 0)
602 : {
603 2064 : nPixelSpace = GDALGetDataTypeSizeBytes(eType);
604 : }
605 2064 : if (nLineSpace == 0)
606 : {
607 2064 : nLineSpace = static_cast<GSpacing>(nXSize) * nPixelSpace;
608 : }
609 2064 : if (nBandSpace == 0)
610 : {
611 2064 : nBandSpace = nYSize * nLineSpace;
612 : }
613 :
614 : GDALRasterizeInfo sInfo;
615 2064 : sInfo.nXSize = nXSize;
616 2064 : sInfo.nYSize = nYSize;
617 2064 : sInfo.nBands = nBands;
618 2064 : sInfo.pabyChunkBuf = pabyChunkBuf;
619 2064 : sInfo.eType = eType;
620 2064 : sInfo.nPixelSpace = nPixelSpace;
621 2064 : sInfo.nLineSpace = nLineSpace;
622 2064 : sInfo.nBandSpace = nBandSpace;
623 2064 : sInfo.eBurnValueType = eBurnValueType;
624 2064 : if (eBurnValueType == GDT_Float64)
625 2063 : 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 2064 : sInfo.eBurnValueSource = eBurnValueSrc;
633 2064 : sInfo.eMergeAlg = eMergeAlg;
634 2064 : sInfo.bFillSetVisitedPoints = false;
635 2064 : 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 4128 : aPointX; // coordinate X values from all rings/components
643 : std::vector<double>
644 4128 : aPointY; // coordinate Y values from all rings/components
645 4128 : std::vector<double> aPointVariant; // coordinate Z values
646 4128 : std::vector<int> aPartSize; // number of X/Y/(Z) values associated with
647 : // each ring/component
648 :
649 2064 : GDALCollectRingsFromGeometry(poShape, aPointX, aPointY, aPointVariant,
650 : aPartSize, eBurnValueSrc);
651 :
652 : /* -------------------------------------------------------------------- */
653 : /* Transform points if needed. */
654 : /* -------------------------------------------------------------------- */
655 2064 : if (pfnTransformer != nullptr)
656 : {
657 : int *panSuccess =
658 2064 : static_cast<int *>(CPLCalloc(sizeof(int), aPointX.size()));
659 :
660 : // TODO: We need to add all appropriate error checking at some point.
661 2064 : pfnTransformer(pTransformArg, FALSE, static_cast<int>(aPointX.size()),
662 : aPointX.data(), aPointY.data(), nullptr, panSuccess);
663 2064 : CPLFree(panSuccess);
664 : }
665 :
666 : /* -------------------------------------------------------------------- */
667 : /* Shift to account for the buffer offset of this buffer. */
668 : /* -------------------------------------------------------------------- */
669 58455 : for (unsigned int i = 0; i < aPointX.size(); i++)
670 56391 : aPointX[i] -= nXOff;
671 58455 : for (unsigned int i = 0; i < aPointY.size(); i++)
672 56391 : 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 2064 : 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 179 : default:
717 : {
718 179 : if (eMergeAlg == GRMA_Add)
719 : {
720 15 : sInfo.bFillSetVisitedPoints = true;
721 15 : sInfo.poSetVisitedPoints = new std::set<uint64_t>();
722 : }
723 179 : 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 83 : if (eBurnValueSrc == GBV_UserBurnValue)
730 : {
731 80 : GDALdllImageLineAllTouched(
732 : sInfo.nXSize, nYSize,
733 80 : static_cast<int>(aPartSize.size()), aPartSize.data(),
734 80 : aPointX.data(), aPointY.data(), nullptr, gvBurnPoint,
735 : &sInfo, eMergeAlg == GRMA_Add, true);
736 : }
737 : else
738 : {
739 6 : for (unsigned int i = 0, n = 0;
740 6 : i < static_cast<unsigned int>(aPartSize.size()); i++)
741 : {
742 17 : for (int j = 0; j < aPartSize[i]; j++)
743 14 : aPointVariant[n++] = aPointVariant[0];
744 : }
745 :
746 3 : GDALdllImageLineAllTouched(
747 : sInfo.nXSize, nYSize,
748 3 : static_cast<int>(aPartSize.size()), aPartSize.data(),
749 3 : aPointX.data(), aPointY.data(), aPointVariant.data(),
750 : gvBurnPoint, &sInfo, eMergeAlg == GRMA_Add, true);
751 : }
752 : }
753 179 : sInfo.bFillSetVisitedPoints = false;
754 358 : GDALdllImageFilledPolygon(
755 179 : sInfo.nXSize, nYSize, static_cast<int>(aPartSize.size()),
756 179 : aPartSize.data(), aPointX.data(), aPointY.data(),
757 : (eBurnValueSrc == GBV_UserBurnValue) ? nullptr
758 6 : : aPointVariant.data(),
759 : gvBurnScanline, &sInfo, eMergeAlg == GRMA_Add);
760 : }
761 179 : break;
762 : }
763 :
764 2064 : 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 120 : static CPLErr GDALRasterizeOptions(CSLConstList papszOptions, int *pbAllTouched,
775 : GDALBurnValueSrc *peBurnValueSource,
776 : GDALRasterMergeAlg *peMergeAlg,
777 : GDALRasterizeOptim *peOptim)
778 : {
779 120 : *pbAllTouched = CPLFetchBool(papszOptions, "ALL_TOUCHED", false);
780 :
781 120 : const char *pszOpt = CSLFetchNameValue(papszOptions, "BURN_VALUE_FROM");
782 120 : *peBurnValueSource = GBV_UserBurnValue;
783 120 : if (pszOpt)
784 : {
785 5 : if (EQUAL(pszOpt, "Z"))
786 : {
787 5 : *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 120 : *peMergeAlg = GRMA_Replace;
803 120 : pszOpt = CSLFetchNameValue(papszOptions, "MERGE_ALG");
804 120 : if (pszOpt)
805 : {
806 18 : if (EQUAL(pszOpt, "ADD"))
807 : {
808 18 : *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 120 : pszOpt = CSLFetchNameValue(papszOptions, "OPTIM");
826 120 : if (pszOpt)
827 : {
828 25 : if (peOptim)
829 : {
830 25 : *peOptim = GRO_Auto;
831 25 : if (EQUAL(pszOpt, "RASTER"))
832 : {
833 1 : *peOptim = GRO_Raster;
834 : }
835 24 : else if (EQUAL(pszOpt, "VECTOR"))
836 : {
837 1 : *peOptim = GRO_Vector;
838 : }
839 23 : else if (EQUAL(pszOpt, "AUTO"))
840 : {
841 23 : *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 120 : 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 77 : 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 77 : VALIDATE_POINTER1(hDS, "GDALRasterizeGeometries", CE_Failure);
989 :
990 77 : return GDALRasterizeGeometriesInternal(
991 : hDS, nBandCount, panBandList, nGeomCount, pahGeometries, pfnTransformer,
992 : pTransformArg, GDT_Float64, padfGeomBurnValues, nullptr, papszOptions,
993 77 : 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 78 : 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 78 : if (pfnProgress == nullptr)
1028 30 : pfnProgress = GDALDummyProgress;
1029 :
1030 78 : GDALDataset *poDS = GDALDataset::FromHandle(hDS);
1031 : /* -------------------------------------------------------------------- */
1032 : /* Do some rudimentary arg checking. */
1033 : /* -------------------------------------------------------------------- */
1034 78 : if (nBandCount == 0 || nGeomCount == 0)
1035 : {
1036 0 : pfnProgress(1.0, "", pProgressArg);
1037 0 : return CE_None;
1038 : }
1039 :
1040 78 : 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 78 : GDALRasterBand *poBand = poDS->GetRasterBand(panBandList[0]);
1059 78 : if (poBand == nullptr)
1060 0 : return CE_Failure;
1061 :
1062 : /* -------------------------------------------------------------------- */
1063 : /* Options */
1064 : /* -------------------------------------------------------------------- */
1065 78 : int bAllTouched = FALSE;
1066 78 : GDALBurnValueSrc eBurnValueSource = GBV_UserBurnValue;
1067 78 : GDALRasterMergeAlg eMergeAlg = GRMA_Replace;
1068 78 : GDALRasterizeOptim eOptim = GRO_Auto;
1069 78 : if (GDALRasterizeOptions(papszOptions, &bAllTouched, &eBurnValueSource,
1070 78 : &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 78 : bool bNeedToFreeTransformer = false;
1084 :
1085 78 : if (pfnTransformer == nullptr)
1086 : {
1087 47 : bNeedToFreeTransformer = true;
1088 :
1089 47 : char **papszTransformerOptions = nullptr;
1090 47 : double adfGeoTransform[6] = {0.0};
1091 47 : if (poDS->GetGeoTransform(adfGeoTransform) != CE_None &&
1092 47 : poDS->GetGCPCount() == 0 && poDS->GetMetadata("RPC") == nullptr)
1093 : {
1094 3 : papszTransformerOptions = CSLSetNameValue(
1095 : papszTransformerOptions, "DST_METHOD", "NO_GEOTRANSFORM");
1096 : }
1097 :
1098 47 : pTransformArg = GDALCreateGenImgProjTransformer2(
1099 : nullptr, hDS, papszTransformerOptions);
1100 47 : CSLDestroy(papszTransformerOptions);
1101 :
1102 47 : pfnTransformer = GDALGenImgProjTransform;
1103 47 : 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 78 : poBand->GetBlockSize(&nXBlockSize, &nYBlockSize);
1117 :
1118 78 : if (eOptim == GRO_Auto)
1119 : {
1120 76 : eOptim = GRO_Raster;
1121 : // TODO make more tests with various inputs/outputs to adjust the
1122 : // parameters
1123 76 : 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 78 : CPLErr eErr = CE_None;
1140 78 : 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 77 : GDALGetNonComplexDataType(poBand->GetRasterDataType());
1151 :
1152 154 : const uint64_t nScanlineBytes = static_cast<uint64_t>(nBandCount) *
1153 77 : poDS->GetRasterXSize() *
1154 77 : 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 77 : atoi(CSLFetchNameValueDef(papszOptions, "CHUNKYSIZE", "0"));
1169 77 : if (nYChunkSize <= 0)
1170 : {
1171 77 : const GIntBig nYChunkSize64 = GDALGetCacheMax64() / nScanlineBytes;
1172 77 : const int knIntMax = std::numeric_limits<int>::max();
1173 77 : nYChunkSize = nYChunkSize64 > knIntMax
1174 77 : ? knIntMax
1175 : : static_cast<int>(nYChunkSize64);
1176 : }
1177 :
1178 77 : if (nYChunkSize < 1)
1179 0 : nYChunkSize = 1;
1180 77 : if (nYChunkSize > poDS->GetRasterYSize())
1181 77 : nYChunkSize = poDS->GetRasterYSize();
1182 :
1183 77 : CPLDebug("GDAL", "Rasterizer operating on %d swaths of %d scanlines.",
1184 77 : DIV_ROUND_UP(poDS->GetRasterYSize(), nYChunkSize),
1185 : nYChunkSize);
1186 :
1187 77 : pabyChunkBuf = static_cast<unsigned char *>(VSI_MALLOC2_VERBOSE(
1188 : nYChunkSize, static_cast<size_t>(nScanlineBytes)));
1189 77 : 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 77 : pfnProgress(0.0, nullptr, pProgressArg);
1202 :
1203 154 : for (int iY = 0; iY < poDS->GetRasterYSize() && eErr == CE_None;
1204 77 : iY += nYChunkSize)
1205 : {
1206 77 : int nThisYChunkSize = nYChunkSize;
1207 77 : if (nThisYChunkSize + iY > poDS->GetRasterYSize())
1208 0 : nThisYChunkSize = poDS->GetRasterYSize() - iY;
1209 :
1210 77 : 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 77 : if (eErr != CE_None)
1215 0 : break;
1216 :
1217 2074 : for (int iShape = 0; iShape < nGeomCount; iShape++)
1218 : {
1219 3995 : gv_rasterize_one_shape(
1220 : pabyChunkBuf, 0, iY, poDS->GetRasterXSize(),
1221 : nThisYChunkSize, nBandCount, eType, 0, 0, 0, bAllTouched,
1222 1997 : OGRGeometry::FromHandle(pahGeometries[iShape]),
1223 : eBurnValueType,
1224 : padfGeomBurnValues
1225 1996 : ? padfGeomBurnValues +
1226 1996 : 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 77 : 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 77 : if (!pfnProgress((iY + nThisYChunkSize) /
1241 77 : 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 1 : const int nXBlocks = DIV_ROUND_UP(poBand->GetXSize(), nXBlockSize);
1263 1 : const int nYBlocks = DIV_ROUND_UP(poBand->GetYSize(), nYBlockSize);
1264 :
1265 : const GDALDataType eType =
1266 1 : poBand->GetRasterDataType() == GDT_Byte ? GDT_Byte : GDT_Float64;
1267 :
1268 1 : const int nPixelSize = nBandCount * GDALGetDataTypeSizeBytes(eType);
1269 :
1270 : // rem: optimized for square blocks
1271 : const GIntBig nbMaxBlocks64 =
1272 1 : GDALGetCacheMax64() / nPixelSize / nYBlockSize / nXBlockSize;
1273 1 : const int knIntMax = std::numeric_limits<int>::max();
1274 : const int nbMaxBlocks = static_cast<int>(
1275 2 : std::min(static_cast<GIntBig>(knIntMax / nPixelSize / nYBlockSize /
1276 : nXBlockSize),
1277 1 : nbMaxBlocks64));
1278 1 : const int nbBlocksX = std::max(
1279 2 : 1,
1280 2 : std::min(static_cast<int>(sqrt(static_cast<double>(nbMaxBlocks))),
1281 1 : nXBlocks));
1282 : const int nbBlocksY =
1283 1 : std::max(1, std::min(nbMaxBlocks / nbBlocksX, nYBlocks));
1284 :
1285 1 : const uint64_t nChunkSize = static_cast<uint64_t>(nXBlockSize) *
1286 1 : nbBlocksX * nYBlockSize * nbBlocksY;
1287 :
1288 : #if SIZEOF_VOIDP < 8
1289 : // Only on 32-bit systems and in pathological cases
1290 : if (nChunkSize > std::numeric_limits<size_t>::max())
1291 : {
1292 : CPLError(CE_Failure, CPLE_OutOfMemory, "Too big raster");
1293 : if (bNeedToFreeTransformer)
1294 : GDALDestroyTransformer(pTransformArg);
1295 : return CE_Failure;
1296 : }
1297 : #endif
1298 :
1299 : pabyChunkBuf = static_cast<unsigned char *>(
1300 1 : VSI_MALLOC2_VERBOSE(nPixelSize, static_cast<size_t>(nChunkSize)));
1301 1 : if (pabyChunkBuf == nullptr)
1302 : {
1303 0 : if (bNeedToFreeTransformer)
1304 0 : GDALDestroyTransformer(pTransformArg);
1305 0 : return CE_Failure;
1306 : }
1307 :
1308 1 : OGREnvelope sRasterEnvelope;
1309 1 : sRasterEnvelope.MinX = 0;
1310 1 : sRasterEnvelope.MinY = 0;
1311 1 : sRasterEnvelope.MaxX = poDS->GetRasterXSize();
1312 1 : sRasterEnvelope.MaxY = poDS->GetRasterYSize();
1313 :
1314 : /* --------------------------------------------------------------------
1315 : */
1316 : /* loop over the vectorial geometries */
1317 : /* --------------------------------------------------------------------
1318 : */
1319 1 : pfnProgress(0.0, nullptr, pProgressArg);
1320 3 : for (int iShape = 0; iShape < nGeomCount; iShape++)
1321 : {
1322 :
1323 : const OGRGeometry *poGeometry =
1324 2 : OGRGeometry::FromHandle(pahGeometries[iShape]);
1325 2 : if (poGeometry == nullptr || poGeometry->IsEmpty())
1326 0 : continue;
1327 : /* ------------------------------------------------------------ */
1328 : /* get the envelope of the geometry and transform it to */
1329 : /* pixels coordinates. */
1330 : /* ------------------------------------------------------------ */
1331 2 : OGREnvelope sGeomEnvelope;
1332 2 : poGeometry->getEnvelope(&sGeomEnvelope);
1333 2 : if (pfnTransformer != nullptr)
1334 : {
1335 2 : int anSuccessTransform[2] = {0};
1336 : double apCorners[4];
1337 2 : apCorners[0] = sGeomEnvelope.MinX;
1338 2 : apCorners[1] = sGeomEnvelope.MaxX;
1339 2 : apCorners[2] = sGeomEnvelope.MinY;
1340 2 : apCorners[3] = sGeomEnvelope.MaxY;
1341 :
1342 2 : if (!pfnTransformer(pTransformArg, FALSE, 2, &(apCorners[0]),
1343 : &(apCorners[2]), nullptr,
1344 2 : anSuccessTransform) ||
1345 2 : !anSuccessTransform[0] || !anSuccessTransform[1])
1346 : {
1347 0 : continue;
1348 : }
1349 2 : sGeomEnvelope.MinX = std::min(apCorners[0], apCorners[1]);
1350 2 : sGeomEnvelope.MaxX = std::max(apCorners[0], apCorners[1]);
1351 2 : sGeomEnvelope.MinY = std::min(apCorners[2], apCorners[3]);
1352 2 : sGeomEnvelope.MaxY = std::max(apCorners[2], apCorners[3]);
1353 : }
1354 2 : if (!sGeomEnvelope.Intersects(sRasterEnvelope))
1355 0 : continue;
1356 2 : sGeomEnvelope.Intersect(sRasterEnvelope);
1357 2 : CPLAssert(sGeomEnvelope.MinX >= 0 &&
1358 : sGeomEnvelope.MinX <= poDS->GetRasterXSize());
1359 2 : CPLAssert(sGeomEnvelope.MinY >= 0 &&
1360 : sGeomEnvelope.MinY <= poDS->GetRasterYSize());
1361 2 : CPLAssert(sGeomEnvelope.MaxX >= 0 &&
1362 : sGeomEnvelope.MaxX <= poDS->GetRasterXSize());
1363 2 : CPLAssert(sGeomEnvelope.MaxY >= 0 &&
1364 : sGeomEnvelope.MaxY <= poDS->GetRasterYSize());
1365 2 : const int minBlockX = int(sGeomEnvelope.MinX) / nXBlockSize;
1366 2 : const int minBlockY = int(sGeomEnvelope.MinY) / nYBlockSize;
1367 2 : const int maxBlockX = int(sGeomEnvelope.MaxX + 1) / nXBlockSize;
1368 2 : const int maxBlockY = int(sGeomEnvelope.MaxY + 1) / nYBlockSize;
1369 :
1370 : /* ------------------------------------------------------------ */
1371 : /* loop over the blocks concerned by the geometry */
1372 : /* (by packs of nbBlocksX x nbBlocksY) */
1373 : /* ------------------------------------------------------------ */
1374 :
1375 5 : for (int xB = minBlockX; xB <= maxBlockX; xB += nbBlocksX)
1376 : {
1377 6 : for (int yB = minBlockY; yB <= maxBlockY; yB += nbBlocksY)
1378 : {
1379 :
1380 : /* --------------------------------------------------------------------
1381 : */
1382 : /* ensure to stay in the image */
1383 : /* --------------------------------------------------------------------
1384 : */
1385 3 : int remSBX = std::min(maxBlockX - xB + 1, nbBlocksX);
1386 3 : int remSBY = std::min(maxBlockY - yB + 1, nbBlocksY);
1387 3 : int nThisXChunkSize = nXBlockSize * remSBX;
1388 3 : int nThisYChunkSize = nYBlockSize * remSBY;
1389 6 : if (xB * nXBlockSize + nThisXChunkSize >
1390 3 : poDS->GetRasterXSize())
1391 1 : nThisXChunkSize =
1392 1 : poDS->GetRasterXSize() - xB * nXBlockSize;
1393 6 : if (yB * nYBlockSize + nThisYChunkSize >
1394 3 : poDS->GetRasterYSize())
1395 2 : nThisYChunkSize =
1396 2 : poDS->GetRasterYSize() - yB * nYBlockSize;
1397 :
1398 : /* --------------------------------------------------------------------
1399 : */
1400 : /* read image / process buffer / write buffer */
1401 : /* --------------------------------------------------------------------
1402 : */
1403 3 : eErr = poDS->RasterIO(
1404 : GF_Read, xB * nXBlockSize, yB * nYBlockSize,
1405 : nThisXChunkSize, nThisYChunkSize, pabyChunkBuf,
1406 : nThisXChunkSize, nThisYChunkSize, eType, nBandCount,
1407 : panBandList, 0, 0, 0, nullptr);
1408 3 : if (eErr != CE_None)
1409 0 : break;
1410 :
1411 6 : gv_rasterize_one_shape(
1412 : pabyChunkBuf, xB * nXBlockSize, yB * nYBlockSize,
1413 : nThisXChunkSize, nThisYChunkSize, nBandCount, eType, 0,
1414 : 0, 0, bAllTouched,
1415 3 : OGRGeometry::FromHandle(pahGeometries[iShape]),
1416 : eBurnValueType,
1417 : padfGeomBurnValues
1418 3 : ? padfGeomBurnValues +
1419 3 : static_cast<size_t>(iShape) * nBandCount
1420 : : nullptr,
1421 : panGeomBurnValues
1422 0 : ? panGeomBurnValues +
1423 0 : static_cast<size_t>(iShape) * nBandCount
1424 : : nullptr,
1425 : eBurnValueSource, eMergeAlg, pfnTransformer,
1426 : pTransformArg);
1427 :
1428 3 : eErr = poDS->RasterIO(
1429 : GF_Write, xB * nXBlockSize, yB * nYBlockSize,
1430 : nThisXChunkSize, nThisYChunkSize, pabyChunkBuf,
1431 : nThisXChunkSize, nThisYChunkSize, eType, nBandCount,
1432 : panBandList, 0, 0, 0, nullptr);
1433 3 : if (eErr != CE_None)
1434 0 : break;
1435 : }
1436 : }
1437 :
1438 2 : if (!pfnProgress(iShape / static_cast<double>(nGeomCount), "",
1439 : pProgressArg))
1440 : {
1441 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1442 0 : eErr = CE_Failure;
1443 : }
1444 : }
1445 :
1446 1 : if (!pfnProgress(1., "", pProgressArg))
1447 : {
1448 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1449 0 : eErr = CE_Failure;
1450 : }
1451 : }
1452 :
1453 : /* -------------------------------------------------------------------- */
1454 : /* cleanup */
1455 : /* -------------------------------------------------------------------- */
1456 78 : VSIFree(pabyChunkBuf);
1457 :
1458 78 : if (bNeedToFreeTransformer)
1459 47 : GDALDestroyTransformer(pTransformArg);
1460 :
1461 78 : return eErr;
1462 : }
1463 :
1464 : /************************************************************************/
1465 : /* GDALRasterizeLayers() */
1466 : /************************************************************************/
1467 :
1468 : /**
1469 : * Burn geometries from the specified list of layers into raster.
1470 : *
1471 : * Rasterize all the geometric objects from a list of layers into a raster
1472 : * dataset. The layers are passed as an array of OGRLayerH handlers.
1473 : *
1474 : * If the geometries are in the georeferenced coordinates of the raster
1475 : * dataset, then the pfnTransform may be passed in NULL and one will be
1476 : * derived internally from the geotransform of the dataset. The transform
1477 : * needs to transform the geometry locations into pixel/line coordinates
1478 : * on the raster dataset.
1479 : *
1480 : * The output raster may be of any GDAL supported datatype. An explicit list
1481 : * of burn values for each layer for each band must be passed in.
1482 : *
1483 : * @param hDS output data, must be opened in update mode.
1484 : * @param nBandCount the number of bands to be updated.
1485 : * @param panBandList the list of bands to be updated.
1486 : * @param nLayerCount the number of layers being passed in pahLayers array.
1487 : * @param pahLayers the array of layers to burn in.
1488 : * @param pfnTransformer transformation to apply to geometries to put into
1489 : * pixel/line coordinates on raster. If NULL a geotransform based one will
1490 : * be created internally.
1491 : * @param pTransformArg callback data for transformer.
1492 : * @param padfLayerBurnValues the array of values to burn into the raster.
1493 : * There should be nBandCount values for each layer.
1494 : * @param papszOptions special options controlling rasterization:
1495 : * <ul>
1496 : * <li>"ATTRIBUTE": Identifies an attribute field on the features to be
1497 : * used for a burn in value. The value will be burned into all output
1498 : * bands. If specified, padfLayerBurnValues will not be used and can be a NULL
1499 : * pointer.</li>
1500 : * <li>"CHUNKYSIZE": The height in lines of the chunk to operate on.
1501 : * The larger the chunk size the less times we need to make a pass through all
1502 : * the shapes. If it is not set or set to zero the default chunk size will be
1503 : * used. Default size will be estimated based on the GDAL cache buffer size
1504 : * using formula: cache_size_bytes/scanline_size_bytes, so the chunk will
1505 : * not exceed the cache.</li>
1506 : * <li>"ALL_TOUCHED": May be set to TRUE to set all pixels touched
1507 : * by the line or polygons, not just those whose center is within the polygon
1508 : * (behavior is unspecified when the polygon is just touching the pixel center)
1509 : * or that are selected by Brezenham's line algorithm. Defaults to FALSE.
1510 : * <li>"BURN_VALUE_FROM": May be set to "Z" to use the Z values of the</li>
1511 : * geometries. The value from padfLayerBurnValues or the attribute field value
1512 : * is added to this before burning. In default case dfBurnValue is burned as it
1513 : * is. This is implemented properly only for points and lines for now. Polygons
1514 : * will be burned using the Z value from the first point. The M value may be
1515 : * supported in the future.</li>
1516 : * <li>"MERGE_ALG": May be REPLACE (the default) or ADD. REPLACE results in
1517 : * overwriting of value, while ADD adds the new value to the existing raster,
1518 : * suitable for heatmaps for instance.</li>
1519 : * </ul>
1520 : * @param pfnProgress the progress function to report completion.
1521 : * @param pProgressArg callback data for progress function.
1522 : *
1523 : * @return CE_None on success or CE_Failure on error.
1524 : */
1525 :
1526 42 : CPLErr GDALRasterizeLayers(GDALDatasetH hDS, int nBandCount, int *panBandList,
1527 : int nLayerCount, OGRLayerH *pahLayers,
1528 : GDALTransformerFunc pfnTransformer,
1529 : void *pTransformArg, double *padfLayerBurnValues,
1530 : char **papszOptions, GDALProgressFunc pfnProgress,
1531 : void *pProgressArg)
1532 :
1533 : {
1534 42 : VALIDATE_POINTER1(hDS, "GDALRasterizeLayers", CE_Failure);
1535 :
1536 42 : if (pfnProgress == nullptr)
1537 42 : pfnProgress = GDALDummyProgress;
1538 :
1539 : /* -------------------------------------------------------------------- */
1540 : /* Do some rudimentary arg checking. */
1541 : /* -------------------------------------------------------------------- */
1542 42 : if (nBandCount == 0 || nLayerCount == 0)
1543 0 : return CE_None;
1544 :
1545 42 : GDALDataset *poDS = GDALDataset::FromHandle(hDS);
1546 :
1547 : // Prototype band.
1548 42 : GDALRasterBand *poBand = poDS->GetRasterBand(panBandList[0]);
1549 42 : if (poBand == nullptr)
1550 0 : return CE_Failure;
1551 :
1552 : /* -------------------------------------------------------------------- */
1553 : /* Options */
1554 : /* -------------------------------------------------------------------- */
1555 42 : int bAllTouched = FALSE;
1556 42 : GDALBurnValueSrc eBurnValueSource = GBV_UserBurnValue;
1557 42 : GDALRasterMergeAlg eMergeAlg = GRMA_Replace;
1558 42 : if (GDALRasterizeOptions(papszOptions, &bAllTouched, &eBurnValueSource,
1559 42 : &eMergeAlg, nullptr) == CE_Failure)
1560 : {
1561 0 : return CE_Failure;
1562 : }
1563 :
1564 : /* -------------------------------------------------------------------- */
1565 : /* Establish a chunksize to operate on. The larger the chunk */
1566 : /* size the less times we need to make a pass through all the */
1567 : /* shapes. */
1568 : /* -------------------------------------------------------------------- */
1569 42 : const char *pszYChunkSize = CSLFetchNameValue(papszOptions, "CHUNKYSIZE");
1570 :
1571 42 : const GDALDataType eType = poBand->GetRasterDataType();
1572 :
1573 : const int nScanlineBytes =
1574 42 : nBandCount * poDS->GetRasterXSize() * GDALGetDataTypeSizeBytes(eType);
1575 :
1576 42 : int nYChunkSize = 0;
1577 42 : if (!(pszYChunkSize && ((nYChunkSize = atoi(pszYChunkSize))) != 0))
1578 : {
1579 42 : const GIntBig nYChunkSize64 = GDALGetCacheMax64() / nScanlineBytes;
1580 42 : nYChunkSize = static_cast<int>(
1581 42 : std::min<GIntBig>(nYChunkSize64, std::numeric_limits<int>::max()));
1582 : }
1583 :
1584 42 : if (nYChunkSize < 1)
1585 0 : nYChunkSize = 1;
1586 42 : if (nYChunkSize > poDS->GetRasterYSize())
1587 42 : nYChunkSize = poDS->GetRasterYSize();
1588 :
1589 42 : CPLDebug("GDAL", "Rasterizer operating on %d swaths of %d scanlines.",
1590 42 : DIV_ROUND_UP(poDS->GetRasterYSize(), nYChunkSize), nYChunkSize);
1591 : unsigned char *pabyChunkBuf = static_cast<unsigned char *>(
1592 42 : VSI_MALLOC2_VERBOSE(nYChunkSize, nScanlineBytes));
1593 42 : if (pabyChunkBuf == nullptr)
1594 : {
1595 0 : return CE_Failure;
1596 : }
1597 :
1598 : /* -------------------------------------------------------------------- */
1599 : /* Read the image once for all layers if user requested to render */
1600 : /* the whole raster in single chunk. */
1601 : /* -------------------------------------------------------------------- */
1602 42 : if (nYChunkSize == poDS->GetRasterYSize())
1603 : {
1604 42 : if (poDS->RasterIO(GF_Read, 0, 0, poDS->GetRasterXSize(), nYChunkSize,
1605 : pabyChunkBuf, poDS->GetRasterXSize(), nYChunkSize,
1606 : eType, nBandCount, panBandList, 0, 0, 0,
1607 42 : nullptr) != CE_None)
1608 : {
1609 0 : CPLFree(pabyChunkBuf);
1610 0 : return CE_Failure;
1611 : }
1612 : }
1613 :
1614 : /* ==================================================================== */
1615 : /* Read the specified layers transforming and rasterizing */
1616 : /* geometries. */
1617 : /* ==================================================================== */
1618 42 : CPLErr eErr = CE_None;
1619 42 : const char *pszBurnAttribute = CSLFetchNameValue(papszOptions, "ATTRIBUTE");
1620 :
1621 42 : pfnProgress(0.0, nullptr, pProgressArg);
1622 :
1623 84 : for (int iLayer = 0; iLayer < nLayerCount; iLayer++)
1624 : {
1625 42 : OGRLayer *poLayer = reinterpret_cast<OGRLayer *>(pahLayers[iLayer]);
1626 :
1627 42 : if (!poLayer)
1628 : {
1629 0 : CPLError(CE_Warning, CPLE_AppDefined,
1630 : "Layer element number %d is NULL, skipping.", iLayer);
1631 0 : continue;
1632 : }
1633 :
1634 : /* --------------------------------------------------------------------
1635 : */
1636 : /* If the layer does not contain any features just skip it. */
1637 : /* Do not force the feature count, so if driver doesn't know */
1638 : /* exact number of features, go down the normal way. */
1639 : /* --------------------------------------------------------------------
1640 : */
1641 42 : if (poLayer->GetFeatureCount(FALSE) == 0)
1642 0 : continue;
1643 :
1644 42 : int iBurnField = -1;
1645 42 : double *padfBurnValues = nullptr;
1646 :
1647 42 : if (pszBurnAttribute)
1648 : {
1649 : iBurnField =
1650 1 : poLayer->GetLayerDefn()->GetFieldIndex(pszBurnAttribute);
1651 1 : if (iBurnField == -1)
1652 : {
1653 0 : CPLError(CE_Warning, CPLE_AppDefined,
1654 : "Failed to find field %s on layer %s, skipping.",
1655 0 : pszBurnAttribute, poLayer->GetLayerDefn()->GetName());
1656 0 : continue;
1657 : }
1658 : }
1659 : else
1660 : {
1661 41 : padfBurnValues = padfLayerBurnValues + iLayer * nBandCount;
1662 : }
1663 :
1664 : /* --------------------------------------------------------------------
1665 : */
1666 : /* If we have no transformer, create the one from input file */
1667 : /* projection. Note that each layer can be georefernced */
1668 : /* separately. */
1669 : /* --------------------------------------------------------------------
1670 : */
1671 42 : bool bNeedToFreeTransformer = false;
1672 :
1673 42 : if (pfnTransformer == nullptr)
1674 : {
1675 42 : char *pszProjection = nullptr;
1676 42 : bNeedToFreeTransformer = true;
1677 :
1678 42 : OGRSpatialReference *poSRS = poLayer->GetSpatialRef();
1679 42 : if (!poSRS)
1680 : {
1681 3 : if (poDS->GetSpatialRef() != nullptr ||
1682 4 : poDS->GetGCPSpatialRef() != nullptr ||
1683 1 : poDS->GetMetadata("RPC") != nullptr)
1684 : {
1685 2 : CPLError(
1686 : CE_Warning, CPLE_AppDefined,
1687 : "Failed to fetch spatial reference on layer %s "
1688 : "to build transformer, assuming matching coordinate "
1689 : "systems.",
1690 2 : poLayer->GetLayerDefn()->GetName());
1691 : }
1692 : }
1693 : else
1694 : {
1695 39 : poSRS->exportToWkt(&pszProjection);
1696 : }
1697 :
1698 42 : char **papszTransformerOptions = nullptr;
1699 42 : if (pszProjection != nullptr)
1700 39 : papszTransformerOptions = CSLSetNameValue(
1701 : papszTransformerOptions, "SRC_SRS", pszProjection);
1702 42 : double adfGeoTransform[6] = {};
1703 42 : if (poDS->GetGeoTransform(adfGeoTransform) != CE_None &&
1704 42 : poDS->GetGCPCount() == 0 && poDS->GetMetadata("RPC") == nullptr)
1705 : {
1706 0 : papszTransformerOptions = CSLSetNameValue(
1707 : papszTransformerOptions, "DST_METHOD", "NO_GEOTRANSFORM");
1708 : }
1709 :
1710 42 : pTransformArg = GDALCreateGenImgProjTransformer2(
1711 : nullptr, hDS, papszTransformerOptions);
1712 42 : pfnTransformer = GDALGenImgProjTransform;
1713 :
1714 42 : CPLFree(pszProjection);
1715 42 : CSLDestroy(papszTransformerOptions);
1716 42 : if (pTransformArg == nullptr)
1717 : {
1718 0 : CPLFree(pabyChunkBuf);
1719 0 : return CE_Failure;
1720 : }
1721 : }
1722 :
1723 42 : poLayer->ResetReading();
1724 :
1725 : /* --------------------------------------------------------------------
1726 : */
1727 : /* Loop over image in designated chunks. */
1728 : /* --------------------------------------------------------------------
1729 : */
1730 :
1731 : double *padfAttrValues = static_cast<double *>(
1732 42 : VSI_MALLOC_VERBOSE(sizeof(double) * nBandCount));
1733 42 : if (padfAttrValues == nullptr)
1734 0 : eErr = CE_Failure;
1735 :
1736 84 : for (int iY = 0; iY < poDS->GetRasterYSize() && eErr == CE_None;
1737 42 : iY += nYChunkSize)
1738 : {
1739 42 : int nThisYChunkSize = nYChunkSize;
1740 42 : if (nThisYChunkSize + iY > poDS->GetRasterYSize())
1741 0 : nThisYChunkSize = poDS->GetRasterYSize() - iY;
1742 :
1743 : // Only re-read image if not a single chunk is being rendered.
1744 42 : if (nYChunkSize < poDS->GetRasterYSize())
1745 : {
1746 0 : eErr = poDS->RasterIO(
1747 : GF_Read, 0, iY, poDS->GetRasterXSize(), nThisYChunkSize,
1748 : pabyChunkBuf, poDS->GetRasterXSize(), nThisYChunkSize,
1749 : eType, nBandCount, panBandList, 0, 0, 0, nullptr);
1750 0 : if (eErr != CE_None)
1751 0 : break;
1752 : }
1753 :
1754 100 : for (auto &poFeat : poLayer)
1755 : {
1756 58 : OGRGeometry *poGeom = poFeat->GetGeometryRef();
1757 :
1758 58 : if (pszBurnAttribute)
1759 : {
1760 : const double dfAttrValue =
1761 5 : poFeat->GetFieldAsDouble(iBurnField);
1762 20 : for (int iBand = 0; iBand < nBandCount; iBand++)
1763 15 : padfAttrValues[iBand] = dfAttrValue;
1764 :
1765 5 : padfBurnValues = padfAttrValues;
1766 : }
1767 :
1768 58 : gv_rasterize_one_shape(
1769 : pabyChunkBuf, 0, iY, poDS->GetRasterXSize(),
1770 : nThisYChunkSize, nBandCount, eType, 0, 0, 0, bAllTouched,
1771 : poGeom, GDT_Float64, padfBurnValues, nullptr,
1772 : eBurnValueSource, eMergeAlg, pfnTransformer, pTransformArg);
1773 : }
1774 :
1775 : // Only write image if not a single chunk is being rendered.
1776 42 : if (nYChunkSize < poDS->GetRasterYSize())
1777 : {
1778 0 : eErr = poDS->RasterIO(
1779 : GF_Write, 0, iY, poDS->GetRasterXSize(), nThisYChunkSize,
1780 : pabyChunkBuf, poDS->GetRasterXSize(), nThisYChunkSize,
1781 : eType, nBandCount, panBandList, 0, 0, 0, nullptr);
1782 : }
1783 :
1784 42 : poLayer->ResetReading();
1785 :
1786 42 : if (!pfnProgress((iY + nThisYChunkSize) /
1787 42 : static_cast<double>(poDS->GetRasterYSize()),
1788 : "", pProgressArg))
1789 : {
1790 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1791 0 : eErr = CE_Failure;
1792 : }
1793 : }
1794 :
1795 42 : VSIFree(padfAttrValues);
1796 :
1797 42 : if (bNeedToFreeTransformer)
1798 : {
1799 42 : GDALDestroyTransformer(pTransformArg);
1800 42 : pTransformArg = nullptr;
1801 42 : pfnTransformer = nullptr;
1802 : }
1803 : }
1804 :
1805 : /* -------------------------------------------------------------------- */
1806 : /* Write out the image once for all layers if user requested */
1807 : /* to render the whole raster in single chunk. */
1808 : /* -------------------------------------------------------------------- */
1809 42 : if (eErr == CE_None && nYChunkSize == poDS->GetRasterYSize())
1810 : {
1811 : eErr =
1812 42 : poDS->RasterIO(GF_Write, 0, 0, poDS->GetRasterXSize(), nYChunkSize,
1813 : pabyChunkBuf, poDS->GetRasterXSize(), nYChunkSize,
1814 : eType, nBandCount, panBandList, 0, 0, 0, nullptr);
1815 : }
1816 :
1817 : /* -------------------------------------------------------------------- */
1818 : /* cleanup */
1819 : /* -------------------------------------------------------------------- */
1820 42 : VSIFree(pabyChunkBuf);
1821 :
1822 42 : return eErr;
1823 : }
1824 :
1825 : /************************************************************************/
1826 : /* GDALRasterizeLayersBuf() */
1827 : /************************************************************************/
1828 :
1829 : /**
1830 : * Burn geometries from the specified list of layer into raster.
1831 : *
1832 : * Rasterize all the geometric objects from a list of layers into supplied
1833 : * raster buffer. The layers are passed as an array of OGRLayerH handlers.
1834 : *
1835 : * If the geometries are in the georeferenced coordinates of the raster
1836 : * dataset, then the pfnTransform may be passed in NULL and one will be
1837 : * derived internally from the geotransform of the dataset. The transform
1838 : * needs to transform the geometry locations into pixel/line coordinates
1839 : * of the target raster.
1840 : *
1841 : * The output raster may be of any GDAL supported datatype(non complex).
1842 : *
1843 : * @param pData pointer to the output data array.
1844 : *
1845 : * @param nBufXSize width of the output data array in pixels.
1846 : *
1847 : * @param nBufYSize height of the output data array in pixels.
1848 : *
1849 : * @param eBufType data type of the output data array.
1850 : *
1851 : * @param nPixelSpace The byte offset from the start of one pixel value in
1852 : * pData to the start of the next pixel value within a scanline. If defaulted
1853 : * (0) the size of the datatype eBufType is used.
1854 : *
1855 : * @param nLineSpace The byte offset from the start of one scanline in
1856 : * pData to the start of the next. If defaulted the size of the datatype
1857 : * eBufType * nBufXSize is used.
1858 : *
1859 : * @param nLayerCount the number of layers being passed in pahLayers array.
1860 : *
1861 : * @param pahLayers the array of layers to burn in.
1862 : *
1863 : * @param pszDstProjection WKT defining the coordinate system of the target
1864 : * raster.
1865 : *
1866 : * @param padfDstGeoTransform geotransformation matrix of the target raster.
1867 : *
1868 : * @param pfnTransformer transformation to apply to geometries to put into
1869 : * pixel/line coordinates on raster. If NULL a geotransform based one will
1870 : * be created internally.
1871 : *
1872 : * @param pTransformArg callback data for transformer.
1873 : *
1874 : * @param dfBurnValue the value to burn into the raster.
1875 : *
1876 : * @param papszOptions special options controlling rasterization:
1877 : * <ul>
1878 : * <li>"ATTRIBUTE": Identifies an attribute field on the features to be
1879 : * used for a burn in value. The value will be burned into all output
1880 : * bands. If specified, padfLayerBurnValues will not be used and can be a NULL
1881 : * pointer.</li>
1882 : * <li>"ALL_TOUCHED": May be set to TRUE to set all pixels touched
1883 : * by the line or polygons, not just those whose center is within the polygon
1884 : * (behavior is unspecified when the polygon is just touching the pixel center)
1885 : * or that are selected by Brezenham's line algorithm. Defaults to FALSE.</li>
1886 : * <li>"BURN_VALUE_FROM": May be set to "Z" to use
1887 : * the Z values of the geometries. dfBurnValue or the attribute field value is
1888 : * added to this before burning. In default case dfBurnValue is burned as it
1889 : * is. This is implemented properly only for points and lines for now. Polygons
1890 : * will be burned using the Z value from the first point. The M value may
1891 : * be supported in the future.</li>
1892 : * <li>"MERGE_ALG": May be REPLACE (the default) or ADD. REPLACE
1893 : * results in overwriting of value, while ADD adds the new value to the
1894 : * existing raster, suitable for heatmaps for instance.</li>
1895 : * </ul>
1896 : *
1897 : * @param pfnProgress the progress function to report completion.
1898 : *
1899 : * @param pProgressArg callback data for progress function.
1900 : *
1901 : *
1902 : * @return CE_None on success or CE_Failure on error.
1903 : */
1904 :
1905 0 : CPLErr GDALRasterizeLayersBuf(
1906 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
1907 : int nPixelSpace, int nLineSpace, int nLayerCount, OGRLayerH *pahLayers,
1908 : const char *pszDstProjection, double *padfDstGeoTransform,
1909 : GDALTransformerFunc pfnTransformer, void *pTransformArg, double dfBurnValue,
1910 : char **papszOptions, GDALProgressFunc pfnProgress, void *pProgressArg)
1911 :
1912 : {
1913 : /* -------------------------------------------------------------------- */
1914 : /* check eType, Avoid not supporting data types */
1915 : /* -------------------------------------------------------------------- */
1916 0 : if (GDALDataTypeIsComplex(eBufType) || eBufType <= GDT_Unknown ||
1917 0 : eBufType >= GDT_TypeCount)
1918 : {
1919 0 : CPLError(CE_Failure, CPLE_NotSupported,
1920 : "GDALRasterizeLayersBuf(): unsupported data type of eBufType");
1921 0 : return CE_Failure;
1922 : }
1923 :
1924 : /* -------------------------------------------------------------------- */
1925 : /* If pixel and line spaceing are defaulted assign reasonable */
1926 : /* value assuming a packed buffer. */
1927 : /* -------------------------------------------------------------------- */
1928 0 : int nTypeSizeBytes = GDALGetDataTypeSizeBytes(eBufType);
1929 0 : if (nPixelSpace == 0)
1930 : {
1931 0 : nPixelSpace = nTypeSizeBytes;
1932 : }
1933 0 : if (nPixelSpace < nTypeSizeBytes)
1934 : {
1935 0 : CPLError(CE_Failure, CPLE_NotSupported,
1936 : "GDALRasterizeLayersBuf(): unsupported value of nPixelSpace");
1937 0 : return CE_Failure;
1938 : }
1939 :
1940 0 : if (nLineSpace == 0)
1941 : {
1942 0 : nLineSpace = nPixelSpace * nBufXSize;
1943 : }
1944 0 : if (nLineSpace < nPixelSpace * nBufXSize)
1945 : {
1946 0 : CPLError(CE_Failure, CPLE_NotSupported,
1947 : "GDALRasterizeLayersBuf(): unsupported value of nLineSpace");
1948 0 : return CE_Failure;
1949 : }
1950 :
1951 0 : if (pfnProgress == nullptr)
1952 0 : pfnProgress = GDALDummyProgress;
1953 :
1954 : /* -------------------------------------------------------------------- */
1955 : /* Do some rudimentary arg checking. */
1956 : /* -------------------------------------------------------------------- */
1957 0 : if (nLayerCount == 0)
1958 0 : return CE_None;
1959 :
1960 : /* -------------------------------------------------------------------- */
1961 : /* Options */
1962 : /* -------------------------------------------------------------------- */
1963 0 : int bAllTouched = FALSE;
1964 0 : GDALBurnValueSrc eBurnValueSource = GBV_UserBurnValue;
1965 0 : GDALRasterMergeAlg eMergeAlg = GRMA_Replace;
1966 0 : if (GDALRasterizeOptions(papszOptions, &bAllTouched, &eBurnValueSource,
1967 0 : &eMergeAlg, nullptr) == CE_Failure)
1968 : {
1969 0 : return CE_Failure;
1970 : }
1971 :
1972 : /* ==================================================================== */
1973 : /* Read the specified layers transforming and rasterizing */
1974 : /* geometries. */
1975 : /* ==================================================================== */
1976 0 : CPLErr eErr = CE_None;
1977 0 : const char *pszBurnAttribute = CSLFetchNameValue(papszOptions, "ATTRIBUTE");
1978 :
1979 0 : pfnProgress(0.0, nullptr, pProgressArg);
1980 :
1981 0 : for (int iLayer = 0; iLayer < nLayerCount; iLayer++)
1982 : {
1983 0 : OGRLayer *poLayer = reinterpret_cast<OGRLayer *>(pahLayers[iLayer]);
1984 :
1985 0 : if (!poLayer)
1986 : {
1987 0 : CPLError(CE_Warning, CPLE_AppDefined,
1988 : "Layer element number %d is NULL, skipping.", iLayer);
1989 0 : continue;
1990 : }
1991 :
1992 : /* --------------------------------------------------------------------
1993 : */
1994 : /* If the layer does not contain any features just skip it. */
1995 : /* Do not force the feature count, so if driver doesn't know */
1996 : /* exact number of features, go down the normal way. */
1997 : /* --------------------------------------------------------------------
1998 : */
1999 0 : if (poLayer->GetFeatureCount(FALSE) == 0)
2000 0 : continue;
2001 :
2002 0 : int iBurnField = -1;
2003 0 : if (pszBurnAttribute)
2004 : {
2005 : iBurnField =
2006 0 : poLayer->GetLayerDefn()->GetFieldIndex(pszBurnAttribute);
2007 0 : if (iBurnField == -1)
2008 : {
2009 0 : CPLError(CE_Warning, CPLE_AppDefined,
2010 : "Failed to find field %s on layer %s, skipping.",
2011 0 : pszBurnAttribute, poLayer->GetLayerDefn()->GetName());
2012 0 : continue;
2013 : }
2014 : }
2015 :
2016 : /* --------------------------------------------------------------------
2017 : */
2018 : /* If we have no transformer, create the one from input file */
2019 : /* projection. Note that each layer can be georefernced */
2020 : /* separately. */
2021 : /* --------------------------------------------------------------------
2022 : */
2023 0 : bool bNeedToFreeTransformer = false;
2024 :
2025 0 : if (pfnTransformer == nullptr)
2026 : {
2027 0 : char *pszProjection = nullptr;
2028 0 : bNeedToFreeTransformer = true;
2029 :
2030 0 : OGRSpatialReference *poSRS = poLayer->GetSpatialRef();
2031 0 : if (!poSRS)
2032 : {
2033 0 : CPLError(CE_Warning, CPLE_AppDefined,
2034 : "Failed to fetch spatial reference on layer %s "
2035 : "to build transformer, assuming matching coordinate "
2036 : "systems.",
2037 0 : poLayer->GetLayerDefn()->GetName());
2038 : }
2039 : else
2040 : {
2041 0 : poSRS->exportToWkt(&pszProjection);
2042 : }
2043 :
2044 0 : pTransformArg = GDALCreateGenImgProjTransformer3(
2045 : pszProjection, nullptr, pszDstProjection, padfDstGeoTransform);
2046 0 : pfnTransformer = GDALGenImgProjTransform;
2047 :
2048 0 : CPLFree(pszProjection);
2049 : }
2050 :
2051 0 : for (auto &poFeat : poLayer)
2052 : {
2053 0 : OGRGeometry *poGeom = poFeat->GetGeometryRef();
2054 :
2055 0 : if (pszBurnAttribute)
2056 0 : dfBurnValue = poFeat->GetFieldAsDouble(iBurnField);
2057 :
2058 0 : gv_rasterize_one_shape(
2059 : static_cast<unsigned char *>(pData), 0, 0, nBufXSize, nBufYSize,
2060 : 1, eBufType, nPixelSpace, nLineSpace, 0, bAllTouched, poGeom,
2061 : GDT_Float64, &dfBurnValue, nullptr, eBurnValueSource, eMergeAlg,
2062 : pfnTransformer, pTransformArg);
2063 : }
2064 :
2065 0 : poLayer->ResetReading();
2066 :
2067 0 : if (!pfnProgress(1, "", pProgressArg))
2068 : {
2069 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
2070 0 : eErr = CE_Failure;
2071 : }
2072 :
2073 0 : if (bNeedToFreeTransformer)
2074 : {
2075 0 : GDALDestroyTransformer(pTransformArg);
2076 0 : pTransformArg = nullptr;
2077 0 : pfnTransformer = nullptr;
2078 : }
2079 : }
2080 :
2081 0 : return eErr;
2082 : }
|
