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