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