1 /* ****************************************************************************
2 *
3 * Project: GDAL Utilities
4 * Purpose: GDAL scattered data gridding (interpolation) tool
5 * Author: Andrey Kiselev, dron@ak4719.spb.edu
6 *
7 * ****************************************************************************
8 * Copyright (c) 2007, Andrey Kiselev <dron@ak4719.spb.edu>
9 * Copyright (c) 2015, Even Rouault <even dot rouault at spatialys dot 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_utils.h"
32 #include "gdal_utils_priv.h"
33 #include "commonutils.h"
34
35 #include <cmath>
36 #include <cstdint>
37 #include <cstdio>
38 #include <cstdlib>
39 #include <algorithm>
40 #include <vector>
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_alg.h"
49 #include "gdal_priv.h"
50 #include "gdalgrid.h"
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 "ogr_srs_api.h"
57 #include "ogrsf_frmts.h"
58
59 CPL_CVSID("$Id: gdal_grid_lib.cpp 3facda23f7b853fb1ad72e4631ed4f70fd6c29c5 2021-04-22 19:00:35 +0200 Even Rouault $")
60
61 /************************************************************************/
62 /* GDALGridOptions */
63 /************************************************************************/
64
65 /** Options for use with GDALGrid(). GDALGridOptions* must be allocated
66 * and freed with GDALGridOptionsNew() and GDALGridOptionsFree() respectively.
67 */
68 struct GDALGridOptions
69 {
70 /*! output format. Use the short format name. */
71 char *pszFormat;
72
73 /*! allow or suppress progress monitor and other non-error output */
74 bool bQuiet;
75
76 /*! the progress function to use */
77 GDALProgressFunc pfnProgress;
78
79 /*! pointer to the progress data variable */
80 void *pProgressData;
81
82 char **papszLayers;
83 char *pszBurnAttribute;
84 double dfIncreaseBurnValue;
85 double dfMultiplyBurnValue;
86 char *pszWHERE;
87 char *pszSQL;
88 GDALDataType eOutputType;
89 char **papszCreateOptions;
90 int nXSize;
91 int nYSize;
92 double dfXRes;
93 double dfYRes;
94 double dfXMin;
95 double dfXMax;
96 double dfYMin;
97 double dfYMax;
98 bool bIsXExtentSet;
99 bool bIsYExtentSet;
100 GDALGridAlgorithm eAlgorithm;
101 void *pOptions;
102 char *pszOutputSRS;
103 OGRGeometry *poSpatialFilter;
104 bool bClipSrc;
105 OGRGeometry *poClipSrc;
106 char *pszClipSrcDS;
107 char *pszClipSrcSQL;
108 char *pszClipSrcLayer;
109 char *pszClipSrcWhere;
110 bool bNoDataSet;
111 double dfNoDataValue;
112 };
113
114 /************************************************************************/
115 /* GetAlgorithmName() */
116 /* */
117 /* Grids algorithm code into mnemonic name. */
118 /************************************************************************/
119
PrintAlgorithmAndOptions(GDALGridAlgorithm eAlgorithm,void * pOptions)120 static void PrintAlgorithmAndOptions( GDALGridAlgorithm eAlgorithm,
121 void *pOptions )
122 {
123 switch ( eAlgorithm )
124 {
125 case GGA_InverseDistanceToAPower:
126 {
127 printf("Algorithm name: \"%s\".\n", szAlgNameInvDist);
128 GDALGridInverseDistanceToAPowerOptions *pOptions2 =
129 static_cast<GDALGridInverseDistanceToAPowerOptions *>(pOptions);
130 CPLprintf("Options are "
131 "\"power=%f:smoothing=%f:radius1=%f:radius2=%f:angle=%f"
132 ":max_points=%lu:min_points=%lu:nodata=%f\"\n",
133 pOptions2->dfPower,
134 pOptions2->dfSmoothing,
135 pOptions2->dfRadius1,
136 pOptions2->dfRadius2,
137 pOptions2->dfAngle,
138 static_cast<unsigned long>(pOptions2->nMaxPoints),
139 static_cast<unsigned long>(pOptions2->nMinPoints),
140 pOptions2->dfNoDataValue);
141 break;
142 }
143 case GGA_InverseDistanceToAPowerNearestNeighbor:
144 {
145 printf("Algorithm name: \"%s\".\n", szAlgNameInvDistNearestNeighbor);
146 GDALGridInverseDistanceToAPowerNearestNeighborOptions *pOptions2 =
147 static_cast<GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(
148 pOptions);
149 CPLprintf("Options are "
150 "\"power=%f:smoothing=%f:radius=%f"
151 ":max_points=%lu:min_points=%lu:nodata=%f\"\n",
152 pOptions2->dfPower,
153 pOptions2->dfSmoothing,
154 pOptions2->dfRadius,
155 static_cast<unsigned long>(pOptions2->nMaxPoints),
156 static_cast<unsigned long>(pOptions2->nMinPoints),
157 pOptions2->dfNoDataValue);
158 break;
159 }
160 case GGA_MovingAverage:
161 {
162 printf("Algorithm name: \"%s\".\n", szAlgNameAverage);
163 GDALGridMovingAverageOptions *pOptions2 =
164 static_cast<GDALGridMovingAverageOptions *>(pOptions);
165 CPLprintf("Options are "
166 "\"radius1=%f:radius2=%f:angle=%f:min_points=%lu"
167 ":nodata=%f\"\n",
168 pOptions2->dfRadius1,
169 pOptions2->dfRadius2,
170 pOptions2->dfAngle,
171 static_cast<unsigned long>(pOptions2->nMinPoints),
172 pOptions2->dfNoDataValue);
173 break;
174 }
175 case GGA_NearestNeighbor:
176 {
177 printf("Algorithm name: \"%s\".\n", szAlgNameNearest);
178 GDALGridNearestNeighborOptions *pOptions2 =
179 static_cast<GDALGridNearestNeighborOptions *>(pOptions);
180 CPLprintf("Options are "
181 "\"radius1=%f:radius2=%f:angle=%f:nodata=%f\"\n",
182 pOptions2->dfRadius1,
183 pOptions2->dfRadius2,
184 pOptions2->dfAngle,
185 pOptions2->dfNoDataValue);
186 break;
187 }
188 case GGA_MetricMinimum:
189 {
190 printf("Algorithm name: \"%s\".\n", szAlgNameMinimum);
191 GDALGridDataMetricsOptions *pOptions2 =
192 static_cast<GDALGridDataMetricsOptions *>(pOptions);
193 CPLprintf("Options are "
194 "\"radius1=%f:radius2=%f:angle=%f:min_points=%lu"
195 ":nodata=%f\"\n",
196 pOptions2->dfRadius1,
197 pOptions2->dfRadius2,
198 pOptions2->dfAngle,
199 static_cast<unsigned long>(pOptions2->nMinPoints),
200 pOptions2->dfNoDataValue);
201 break;
202 }
203 case GGA_MetricMaximum:
204 {
205 printf("Algorithm name: \"%s\".\n", szAlgNameMaximum);
206 GDALGridDataMetricsOptions *pOptions2 =
207 static_cast<GDALGridDataMetricsOptions *>(pOptions);
208 CPLprintf("Options are "
209 "\"radius1=%f:radius2=%f:angle=%f:min_points=%lu"
210 ":nodata=%f\"\n",
211 pOptions2->dfRadius1,
212 pOptions2->dfRadius2,
213 pOptions2->dfAngle,
214 static_cast<unsigned long>(pOptions2->nMinPoints),
215 pOptions2->dfNoDataValue);
216 break;
217 }
218 case GGA_MetricRange:
219 {
220 printf("Algorithm name: \"%s\".\n", szAlgNameRange);
221 GDALGridDataMetricsOptions *pOptions2 =
222 static_cast<GDALGridDataMetricsOptions *>(pOptions);
223 CPLprintf("Options are "
224 "\"radius1=%f:radius2=%f:angle=%f:min_points=%lu"
225 ":nodata=%f\"\n",
226 pOptions2->dfRadius1,
227 pOptions2->dfRadius2,
228 pOptions2->dfAngle,
229 static_cast<unsigned long>(pOptions2->nMinPoints),
230 pOptions2->dfNoDataValue);
231 break;
232 }
233 case GGA_MetricCount:
234 {
235 printf("Algorithm name: \"%s\".\n", szAlgNameCount);
236 GDALGridDataMetricsOptions *pOptions2 =
237 static_cast<GDALGridDataMetricsOptions *>(pOptions);
238 CPLprintf("Options are "
239 "\"radius1=%f:radius2=%f:angle=%f:min_points=%lu"
240 ":nodata=%f\"\n",
241 pOptions2->dfRadius1,
242 pOptions2->dfRadius2,
243 pOptions2->dfAngle,
244 static_cast<unsigned long>(pOptions2->nMinPoints),
245 pOptions2->dfNoDataValue);
246 break;
247 }
248 case GGA_MetricAverageDistance:
249 {
250 printf("Algorithm name: \"%s\".\n", szAlgNameAverageDistance);
251 GDALGridDataMetricsOptions *pOptions2 =
252 static_cast<GDALGridDataMetricsOptions *>(pOptions);
253 CPLprintf("Options are "
254 "\"radius1=%f:radius2=%f:angle=%f:min_points=%lu"
255 ":nodata=%f\"\n",
256 pOptions2->dfRadius1,
257 pOptions2->dfRadius2,
258 pOptions2->dfAngle,
259 static_cast<unsigned long>(pOptions2->nMinPoints),
260 pOptions2->dfNoDataValue);
261 break;
262 }
263 case GGA_MetricAverageDistancePts:
264 {
265 printf("Algorithm name: \"%s\".\n", szAlgNameAverageDistancePts);
266 GDALGridDataMetricsOptions *pOptions2 =
267 static_cast<GDALGridDataMetricsOptions *>(pOptions);
268 CPLprintf("Options are "
269 "\"radius1=%f:radius2=%f:angle=%f:min_points=%lu"
270 ":nodata=%f\"\n",
271 pOptions2->dfRadius1,
272 pOptions2->dfRadius2,
273 pOptions2->dfAngle,
274 static_cast<unsigned long>(pOptions2->nMinPoints),
275 pOptions2->dfNoDataValue);
276 break;
277 }
278 case GGA_Linear:
279 {
280 printf("Algorithm name: \"%s\".\n", szAlgNameLinear );
281 GDALGridLinearOptions *pOptions2 =
282 static_cast<GDALGridLinearOptions *>(pOptions);
283 CPLprintf("Options are "
284 "\"radius=%f:nodata=%f\"\n",
285 pOptions2->dfRadius,
286 pOptions2->dfNoDataValue);
287 break;
288 }
289 default:
290 {
291 printf("Algorithm is unknown.\n");
292 break;
293 }
294 }
295 }
296
297 /************************************************************************/
298 /* ProcessGeometry() */
299 /* */
300 /* Extract point coordinates from the geometry reference and set the */
301 /* Z value as requested. Test whether we are in the clipped region */
302 /* before processing. */
303 /************************************************************************/
304
ProcessGeometry(OGRPoint * poGeom,OGRGeometry * poClipSrc,int iBurnField,double dfBurnValue,const double dfIncreaseBurnValue,const double dfMultiplyBurnValue,std::vector<double> & adfX,std::vector<double> & adfY,std::vector<double> & adfZ)305 static void ProcessGeometry( OGRPoint *poGeom, OGRGeometry *poClipSrc,
306 int iBurnField, double dfBurnValue,
307 const double dfIncreaseBurnValue,
308 const double dfMultiplyBurnValue,
309 std::vector<double> &adfX,
310 std::vector<double> &adfY,
311 std::vector<double> &adfZ )
312
313 {
314 if ( poClipSrc && !poGeom->Within(poClipSrc) )
315 return;
316
317 adfX.push_back( poGeom->getX() );
318 adfY.push_back( poGeom->getY() );
319 if ( iBurnField < 0 )
320 adfZ.push_back( (poGeom->getZ() + dfIncreaseBurnValue) * dfMultiplyBurnValue );
321 else
322 adfZ.push_back( (dfBurnValue + dfIncreaseBurnValue) * dfMultiplyBurnValue );
323 }
324
325 /************************************************************************/
326 /* ProcessCommonGeometry() */
327 /* */
328 /* Process recursively geometry and extract points. */
329 /************************************************************************/
330
ProcessCommonGeometry(OGRGeometry * poGeom,OGRGeometry * poClipSrc,int iBurnField,double dfBurnValue,const double dfIncreaseBurnValue,const double dfMultiplyBurnValue,std::vector<double> & adfX,std::vector<double> & adfY,std::vector<double> & adfZ)331 static void ProcessCommonGeometry(OGRGeometry* poGeom, OGRGeometry *poClipSrc,
332 int iBurnField, double dfBurnValue,
333 const double dfIncreaseBurnValue,
334 const double dfMultiplyBurnValue,
335 std::vector<double> &adfX,
336 std::vector<double> &adfY,
337 std::vector<double> &adfZ)
338 {
339 if (nullptr == poGeom)
340 return;
341
342 OGRwkbGeometryType eType = wkbFlatten(poGeom->getGeometryType());
343 switch (eType)
344 {
345 case wkbPoint:
346 return ProcessGeometry(poGeom->toPoint(), poClipSrc,
347 iBurnField, dfBurnValue, dfIncreaseBurnValue, dfMultiplyBurnValue, adfX, adfY, adfZ);
348 case wkbLinearRing:
349 case wkbLineString:
350 {
351 OGRLineString *poLS = poGeom->toLineString();
352 OGRPoint point;
353 for (int pointIndex = 0; pointIndex < poLS->getNumPoints(); pointIndex++)
354 {
355 poLS->getPoint(pointIndex, &point);
356 ProcessCommonGeometry(&point, poClipSrc,
357 iBurnField, dfBurnValue, dfIncreaseBurnValue, dfMultiplyBurnValue, adfX, adfY, adfZ);
358 }
359 }
360 break;
361 case wkbPolygon:
362 {
363 OGRPolygon* poPoly = poGeom->toPolygon();
364 OGRLinearRing* poRing = poPoly->getExteriorRing();
365 ProcessCommonGeometry(poRing, poClipSrc,
366 iBurnField, dfBurnValue, dfIncreaseBurnValue, dfMultiplyBurnValue, adfX, adfY, adfZ);
367
368 const int nRings = poPoly->getNumInteriorRings();
369 if (nRings > 0)
370 {
371 for (int ir = 0; ir < nRings; ++ir)
372 {
373 poRing = poPoly->getInteriorRing(ir);
374 ProcessCommonGeometry(poRing, poClipSrc,
375 iBurnField, dfBurnValue, dfIncreaseBurnValue, dfMultiplyBurnValue, adfX, adfY, adfZ);
376 }
377 }
378 }
379 break;
380 case wkbMultiPoint:
381 case wkbMultiPolygon:
382 case wkbMultiLineString:
383 case wkbGeometryCollection:
384 {
385 OGRGeometryCollection* pOGRGeometryCollection = poGeom->toGeometryCollection();
386 for (int i = 0; i < pOGRGeometryCollection->getNumGeometries(); ++i)
387 {
388 ProcessCommonGeometry(pOGRGeometryCollection->getGeometryRef(i), poClipSrc,
389 iBurnField, dfBurnValue, dfIncreaseBurnValue, dfMultiplyBurnValue, adfX, adfY, adfZ);
390 }
391 }
392 break;
393 case wkbUnknown:
394 case wkbNone:
395 default:
396 break;
397 }
398 }
399
400 /************************************************************************/
401 /* ProcessLayer() */
402 /* */
403 /* Process all the features in a layer selection, collecting */
404 /* geometries and burn values. */
405 /************************************************************************/
406
ProcessLayer(OGRLayerH hSrcLayer,GDALDatasetH hDstDS,OGRGeometry * poClipSrc,int nXSize,int nYSize,int nBand,bool & bIsXExtentSet,bool & bIsYExtentSet,double & dfXMin,double & dfXMax,double & dfYMin,double & dfYMax,const char * pszBurnAttribute,const double dfIncreaseBurnValue,const double dfMultiplyBurnValue,GDALDataType eType,GDALGridAlgorithm eAlgorithm,void * pOptions,bool bQuiet,GDALProgressFunc pfnProgress,void * pProgressData)407 static CPLErr ProcessLayer( OGRLayerH hSrcLayer, GDALDatasetH hDstDS,
408 OGRGeometry *poClipSrc,
409 int nXSize, int nYSize, int nBand,
410 bool& bIsXExtentSet, bool& bIsYExtentSet,
411 double& dfXMin, double& dfXMax,
412 double& dfYMin, double& dfYMax,
413 const char *pszBurnAttribute,
414 const double dfIncreaseBurnValue,
415 const double dfMultiplyBurnValue,
416 GDALDataType eType,
417 GDALGridAlgorithm eAlgorithm, void *pOptions,
418 bool bQuiet, GDALProgressFunc pfnProgress,
419 void* pProgressData )
420
421 {
422 /* -------------------------------------------------------------------- */
423 /* Get field index, and check. */
424 /* -------------------------------------------------------------------- */
425 int iBurnField = -1;
426
427 if ( pszBurnAttribute )
428 {
429 iBurnField = OGR_FD_GetFieldIndex( OGR_L_GetLayerDefn( hSrcLayer ),
430 pszBurnAttribute );
431 if( iBurnField == -1 )
432 {
433 printf( "Failed to find field %s on layer %s, skipping.\n",
434 pszBurnAttribute,
435 OGR_FD_GetName( OGR_L_GetLayerDefn( hSrcLayer ) ) );
436 return CE_Failure;
437 }
438 }
439
440 /* -------------------------------------------------------------------- */
441 /* Collect the geometries from this layer, and build list of */
442 /* values to be interpolated. */
443 /* -------------------------------------------------------------------- */
444 OGRFeature *poFeat;
445 std::vector<double> adfX, adfY, adfZ;
446
447 OGR_L_ResetReading( hSrcLayer );
448
449 while( (poFeat = reinterpret_cast<OGRFeature*>(OGR_L_GetNextFeature( hSrcLayer ))) != nullptr )
450 {
451 OGRGeometry *poGeom = poFeat->GetGeometryRef();
452 double dfBurnValue = 0.0;
453
454 if ( iBurnField >= 0 )
455 dfBurnValue = poFeat->GetFieldAsDouble( iBurnField );
456
457 ProcessCommonGeometry(poGeom, poClipSrc, iBurnField, dfBurnValue,
458 dfIncreaseBurnValue, dfMultiplyBurnValue, adfX, adfY, adfZ);
459
460 OGRFeature::DestroyFeature( poFeat );
461 }
462
463 if ( adfX.empty() )
464 {
465 printf( "No point geometry found on layer %s, skipping.\n",
466 OGR_FD_GetName( OGR_L_GetLayerDefn( hSrcLayer ) ) );
467 return CE_None;
468 }
469
470 /* -------------------------------------------------------------------- */
471 /* Compute grid geometry. */
472 /* -------------------------------------------------------------------- */
473 if ( !bIsXExtentSet || !bIsYExtentSet )
474 {
475 OGREnvelope sEnvelope;
476 OGR_L_GetExtent( hSrcLayer, &sEnvelope, TRUE );
477
478 if ( !bIsXExtentSet )
479 {
480 dfXMin = sEnvelope.MinX;
481 dfXMax = sEnvelope.MaxX;
482 bIsXExtentSet = true;
483 }
484
485 if ( !bIsYExtentSet )
486 {
487 dfYMin = sEnvelope.MinY;
488 dfYMax = sEnvelope.MaxY;
489 bIsYExtentSet = true;
490 }
491 }
492
493 /* -------------------------------------------------------------------- */
494 /* Perform gridding. */
495 /* -------------------------------------------------------------------- */
496
497 const double dfDeltaX = (dfXMax - dfXMin) / nXSize;
498 const double dfDeltaY = (dfYMax - dfYMin) / nYSize;
499
500 if ( !bQuiet )
501 {
502 printf( "Grid data type is \"%s\"\n", GDALGetDataTypeName(eType) );
503 printf("Grid size = (%lu %lu).\n",
504 static_cast<unsigned long>(nXSize),
505 static_cast<unsigned long>(nYSize));
506 CPLprintf( "Corner coordinates = (%f %f)-(%f %f).\n",
507 dfXMin - dfDeltaX / 2, dfYMax + dfDeltaY / 2,
508 dfXMax + dfDeltaX / 2, dfYMin - dfDeltaY / 2 );
509 CPLprintf( "Grid cell size = (%f %f).\n", dfDeltaX, dfDeltaY );
510 printf("Source point count = %lu.\n",
511 static_cast<unsigned long>(adfX.size()));
512 PrintAlgorithmAndOptions( eAlgorithm, pOptions );
513 printf("\n");
514 }
515
516 GDALRasterBandH hBand = GDALGetRasterBand( hDstDS, nBand );
517
518 int nBlockXSize = 0;
519 int nBlockYSize = 0;
520 const int nDataTypeSize = GDALGetDataTypeSizeBytes(eType);
521
522 // Try to grow the work buffer up to 16 MB if it is smaller
523 GDALGetBlockSize( hBand, &nBlockXSize, &nBlockYSize );
524 if( nXSize == 0 || nYSize == 0 || nBlockXSize == 0 || nBlockYSize == 0 )
525 return CE_Failure;
526
527 const int nDesiredBufferSize = 16*1024*1024;
528 if( nBlockXSize < nXSize && nBlockYSize < nYSize &&
529 nBlockXSize < nDesiredBufferSize / (nBlockYSize * nDataTypeSize) )
530 {
531 const int nNewBlockXSize =
532 nDesiredBufferSize / (nBlockYSize * nDataTypeSize);
533 nBlockXSize = (nNewBlockXSize / nBlockXSize) * nBlockXSize;
534 if( nBlockXSize > nXSize )
535 nBlockXSize = nXSize;
536 }
537 else if( nBlockXSize == nXSize && nBlockYSize < nYSize &&
538 nBlockYSize < nDesiredBufferSize / (nXSize * nDataTypeSize) )
539 {
540 const int nNewBlockYSize =
541 nDesiredBufferSize / (nXSize * nDataTypeSize);
542 nBlockYSize = (nNewBlockYSize / nBlockYSize) * nBlockYSize;
543 if( nBlockYSize > nYSize )
544 nBlockYSize = nYSize;
545 }
546 CPLDebug("GDAL_GRID", "Work buffer: %d * %d", nBlockXSize, nBlockYSize);
547
548 void *pData = VSIMalloc3(nBlockXSize, nBlockYSize, nDataTypeSize);
549 if( pData == nullptr )
550 {
551 CPLError(CE_Failure, CPLE_OutOfMemory, "Cannot allocate work buffer");
552 return CE_Failure;
553 }
554
555 int nBlock = 0;
556 const int nBlockCount =
557 ((nXSize + nBlockXSize - 1) / nBlockXSize) *
558 ((nYSize + nBlockYSize - 1) / nBlockYSize);
559
560 GDALGridContext* psContext = GDALGridContextCreate( eAlgorithm, pOptions,
561 static_cast<int>(adfX.size()),
562 &(adfX[0]), &(adfY[0]), &(adfZ[0]),
563 TRUE );
564 if( psContext == nullptr )
565 {
566 CPLFree( pData );
567 return CE_Failure;
568 }
569
570 CPLErr eErr = CE_None;
571 for ( int nYOffset = 0;
572 nYOffset < nYSize && eErr == CE_None;
573 nYOffset += nBlockYSize )
574 {
575 for ( int nXOffset = 0;
576 nXOffset < nXSize && eErr == CE_None;
577 nXOffset += nBlockXSize )
578 {
579 void *pScaledProgress = GDALCreateScaledProgress(
580 static_cast<double>(nBlock) / nBlockCount,
581 static_cast<double>(nBlock + 1) / nBlockCount,
582 pfnProgress, pProgressData);
583 nBlock ++;
584
585 int nXRequest = nBlockXSize;
586 if (nXOffset + nXRequest > nXSize)
587 nXRequest = nXSize - nXOffset;
588
589 int nYRequest = nBlockYSize;
590 if (nYOffset + nYRequest > nYSize)
591 nYRequest = nYSize - nYOffset;
592
593 eErr = GDALGridContextProcess( psContext,
594 dfXMin + dfDeltaX * nXOffset,
595 dfXMin + dfDeltaX * (nXOffset + nXRequest),
596 dfYMin + dfDeltaY * nYOffset,
597 dfYMin + dfDeltaY * (nYOffset + nYRequest),
598 nXRequest, nYRequest, eType, pData,
599 GDALScaledProgress, pScaledProgress );
600
601 if( eErr == CE_None )
602 eErr = GDALRasterIO( hBand, GF_Write, nXOffset, nYOffset,
603 nXRequest, nYRequest, pData,
604 nXRequest, nYRequest, eType, 0, 0 );
605
606 GDALDestroyScaledProgress( pScaledProgress );
607 }
608 }
609
610 GDALGridContextFree(psContext);
611
612 CPLFree( pData );
613 return eErr;
614 }
615
616 /************************************************************************/
617 /* LoadGeometry() */
618 /* */
619 /* Read geometries from the given dataset using specified filters and */
620 /* returns a collection of read geometries. */
621 /************************************************************************/
622
LoadGeometry(const char * pszDS,const char * pszSQL,const char * pszLyr,const char * pszWhere)623 static OGRGeometryCollection* LoadGeometry( const char* pszDS,
624 const char* pszSQL,
625 const char* pszLyr,
626 const char* pszWhere )
627 {
628 GDALDataset *poDS = static_cast<GDALDataset*>(
629 GDALOpenEx(pszDS, GDAL_OF_VECTOR, nullptr, nullptr, nullptr));
630 if ( poDS == nullptr )
631 return nullptr;
632
633 OGRLayer *poLyr = nullptr;
634 if ( pszSQL != nullptr )
635 poLyr = poDS->ExecuteSQL( pszSQL, nullptr, nullptr );
636 else if ( pszLyr != nullptr )
637 poLyr = poDS->GetLayerByName( pszLyr );
638 else
639 poLyr = poDS->GetLayer(0);
640
641 if ( poLyr == nullptr )
642 {
643 CPLError(CE_Failure, CPLE_AppDefined,
644 "Failed to identify source layer from datasource." );
645 GDALClose(poDS);
646 return nullptr;
647 }
648
649 if ( pszWhere )
650 poLyr->SetAttributeFilter( pszWhere );
651
652 OGRGeometryCollection *poGeom = nullptr;
653 OGRFeature *poFeat = nullptr;
654 while ( (poFeat = poLyr->GetNextFeature()) != nullptr )
655 {
656 OGRGeometry* poSrcGeom = poFeat->GetGeometryRef();
657 if ( poSrcGeom )
658 {
659 const OGRwkbGeometryType eType =
660 wkbFlatten(poSrcGeom->getGeometryType());
661
662 if ( poGeom == nullptr )
663 poGeom = new OGRMultiPolygon();
664
665 if ( eType == wkbPolygon )
666 {
667 poGeom->addGeometry( poSrcGeom );
668 }
669 else if ( eType == wkbMultiPolygon )
670 {
671 const int nGeomCount =
672 static_cast<OGRMultiPolygon *>(poSrcGeom)->
673 getNumGeometries();
674
675 for( int iGeom = 0; iGeom < nGeomCount; iGeom++ )
676 {
677 poGeom->addGeometry(
678 static_cast<OGRMultiPolygon *>(poSrcGeom)->
679 getGeometryRef(iGeom) );
680 }
681 }
682 else
683 {
684 CPLError(CE_Failure, CPLE_AppDefined,
685 "Geometry not of polygon type.");
686 OGRGeometryFactory::destroyGeometry( poGeom );
687 OGRFeature::DestroyFeature( poFeat );
688 if ( pszSQL != nullptr )
689 poDS->ReleaseResultSet( poLyr );
690 GDALClose(poDS);
691 return nullptr;
692 }
693 }
694
695 OGRFeature::DestroyFeature( poFeat );
696 }
697
698 if( pszSQL != nullptr )
699 poDS->ReleaseResultSet( poLyr );
700 GDALClose(poDS);
701
702 return poGeom;
703 }
704
705 /************************************************************************/
706 /* GDALGrid() */
707 /************************************************************************/
708
709 /**
710 * Create raster from the scattered data.
711 *
712 * This is the equivalent of the <a href="/programs/gdal_grid.html">gdal_grid</a> utility.
713 *
714 * GDALGridOptions* must be allocated and freed with GDALGridOptionsNew()
715 * and GDALGridOptionsFree() respectively.
716 *
717 * @param pszDest the destination dataset path.
718 * @param hSrcDataset the source dataset handle.
719 * @param psOptionsIn the options struct returned by GDALGridOptionsNew() or NULL.
720 * @param pbUsageError pointer to a integer output variable to store if any usage error has occurred or NULL.
721 * @return the output dataset (new dataset that must be closed using GDALClose()) or NULL in case of error.
722 *
723 * @since GDAL 2.1
724 */
725
GDALGrid(const char * pszDest,GDALDatasetH hSrcDataset,const GDALGridOptions * psOptionsIn,int * pbUsageError)726 GDALDatasetH GDALGrid( const char *pszDest, GDALDatasetH hSrcDataset,
727 const GDALGridOptions *psOptionsIn, int *pbUsageError )
728
729 {
730 if( hSrcDataset == nullptr )
731 {
732 CPLError( CE_Failure, CPLE_AppDefined, "No source dataset specified.");
733
734 if(pbUsageError)
735 *pbUsageError = TRUE;
736 return nullptr;
737 }
738 if( pszDest == nullptr )
739 {
740 CPLError( CE_Failure, CPLE_AppDefined, "No target dataset specified.");
741
742 if(pbUsageError)
743 *pbUsageError = TRUE;
744 return nullptr;
745 }
746
747 GDALGridOptions* psOptionsToFree = nullptr;
748 const GDALGridOptions* psOptions = psOptionsIn;
749 if( psOptions == nullptr )
750 {
751 psOptionsToFree = GDALGridOptionsNew(nullptr, nullptr);
752 psOptions = psOptionsToFree;
753 }
754
755 GDALDataset* poSrcDS = static_cast<GDALDataset *>(hSrcDataset);
756
757 if( psOptions->pszSQL == nullptr && psOptions->papszLayers == nullptr &&
758 poSrcDS->GetLayerCount() != 1 )
759 {
760 CPLError(CE_Failure, CPLE_NotSupported,
761 "Neither -sql nor -l are specified, but the source dataset has not one single layer.");
762 if( pbUsageError )
763 *pbUsageError = TRUE;
764 GDALGridOptionsFree(psOptionsToFree);
765 return nullptr;
766 }
767
768 if ( (psOptions->nXSize != 0 || psOptions->nYSize != 0) && (psOptions->dfXRes != 0 || psOptions->dfYRes != 0) ) {
769 CPLError(CE_Failure, CPLE_IllegalArg, "-outsize and -tr options cannot be used at the same time.");
770 GDALGridOptionsFree(psOptionsToFree);
771 return nullptr;
772 }
773
774 /* -------------------------------------------------------------------- */
775 /* Find the output driver. */
776 /* -------------------------------------------------------------------- */
777 CPLString osFormat;
778 if( psOptions->pszFormat == nullptr )
779 {
780 osFormat = GetOutputDriverForRaster(pszDest);
781 if( osFormat.empty() )
782 {
783 GDALGridOptionsFree(psOptionsToFree);
784 return nullptr;
785 }
786 }
787 else
788 {
789 osFormat = psOptions->pszFormat;
790 }
791
792 GDALDriverH hDriver = GDALGetDriverByName( osFormat );
793 if( hDriver == nullptr )
794 {
795 CPLError(CE_Failure, CPLE_AppDefined,
796 "Output driver `%s' not recognised.", osFormat.c_str() );
797 fprintf( stderr,
798 "The following format drivers are configured and support output:\n" );
799 for( int iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
800 {
801 hDriver = GDALGetDriver(iDr);
802
803 if( GDALGetMetadataItem( hDriver, GDAL_DCAP_RASTER, nullptr) != nullptr &&
804 ( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, nullptr ) != nullptr
805 || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY, nullptr ) != nullptr) )
806 {
807 fprintf( stderr, " %s: %s\n",
808 GDALGetDriverShortName( hDriver ),
809 GDALGetDriverLongName( hDriver ) );
810 }
811 }
812 printf( "\n" );
813 GDALGridOptionsFree(psOptionsToFree);
814 return nullptr;
815 }
816
817 /* -------------------------------------------------------------------- */
818 /* Create target raster file. */
819 /* -------------------------------------------------------------------- */
820 int nLayerCount = CSLCount(psOptions->papszLayers);
821 if( nLayerCount == 0 && psOptions->pszSQL == nullptr )
822 nLayerCount = 1; /* due to above check */
823
824 int nBands = nLayerCount;
825
826 if ( psOptions->pszSQL )
827 nBands++;
828
829 int nXSize;
830 int nYSize;
831 if ( psOptions->dfXRes != 0 && psOptions->dfYRes != 0 )
832 {
833 if ((psOptions->dfXMax == psOptions->dfXMin) || (psOptions->dfYMax == psOptions->dfYMin)) {
834 CPLError( CE_Failure, CPLE_IllegalArg,
835 "Invalid txe or tye parameters detected. Please check your -txe or -tye argument.");
836
837 if(pbUsageError)
838 *pbUsageError = TRUE;
839 GDALGridOptionsFree(psOptionsToFree);
840 return nullptr;
841 }
842
843 double dfXSize = (std::fabs(psOptions->dfXMax - psOptions->dfXMin) + (psOptions->dfXRes/2.0)) /
844 psOptions->dfXRes;
845 double dfYSize = (std::fabs(psOptions->dfYMax - psOptions->dfYMin) + (psOptions->dfYRes/2.0)) /
846 psOptions->dfYRes;
847
848 if (dfXSize >= 1 && dfXSize <= INT_MAX && dfYSize >= 1 && dfYSize <= INT_MAX) {
849 nXSize = static_cast<int>(dfXSize);
850 nYSize = static_cast<int>(dfYSize);
851 } else {
852 CPLError( CE_Failure, CPLE_IllegalArg, "Invalid output size detected. Please check your -tr argument");
853
854 if(pbUsageError)
855 *pbUsageError = TRUE;
856 GDALGridOptionsFree(psOptionsToFree);
857 return nullptr;
858 }
859 }
860 else
861 {
862 // FIXME
863 nXSize = psOptions->nXSize;
864 if ( nXSize == 0 )
865 nXSize = 256;
866 nYSize = psOptions->nYSize;
867 if ( nYSize == 0 )
868 nYSize = 256;
869 }
870
871 GDALDatasetH hDstDS =
872 GDALCreate(hDriver, pszDest, nXSize, nYSize, nBands,
873 psOptions->eOutputType, psOptions->papszCreateOptions);
874 if ( hDstDS == nullptr )
875 {
876 GDALGridOptionsFree(psOptionsToFree);
877 return nullptr;
878 }
879
880 if( psOptions->bNoDataSet )
881 {
882 for( int i = 1; i <= nBands; i++ )
883 {
884 GDALRasterBandH hBand = GDALGetRasterBand( hDstDS, i );
885 GDALSetRasterNoDataValue( hBand, psOptions->dfNoDataValue );
886 }
887 }
888
889 double dfXMin = psOptions->dfXMin;
890 double dfYMin = psOptions->dfYMin;
891 double dfXMax = psOptions->dfXMax;
892 double dfYMax = psOptions->dfYMax;
893 bool bIsXExtentSet = psOptions->bIsXExtentSet;
894 bool bIsYExtentSet = psOptions->bIsYExtentSet;
895 CPLErr eErr = CE_None;
896
897 /* -------------------------------------------------------------------- */
898 /* Process SQL request. */
899 /* -------------------------------------------------------------------- */
900
901 if( psOptions->pszSQL != nullptr )
902 {
903 OGRLayer* poLayer = poSrcDS->ExecuteSQL(psOptions->pszSQL,
904 psOptions->poSpatialFilter, nullptr );
905 if( poLayer != nullptr )
906 {
907 // Custom layer will be rasterized in the first band.
908 eErr = ProcessLayer( reinterpret_cast<OGRLayerH>(poLayer), hDstDS, psOptions->poSpatialFilter,
909 nXSize, nYSize, 1,
910 bIsXExtentSet, bIsYExtentSet,
911 dfXMin, dfXMax, dfYMin, dfYMax, psOptions->pszBurnAttribute,
912 psOptions->dfIncreaseBurnValue, psOptions->dfMultiplyBurnValue,
913 psOptions->eOutputType, psOptions->eAlgorithm, psOptions->pOptions,
914 psOptions->bQuiet, psOptions->pfnProgress, psOptions->pProgressData );
915
916 poSrcDS->ReleaseResultSet(poLayer);
917 }
918 }
919
920 /* -------------------------------------------------------------------- */
921 /* Process each layer. */
922 /* -------------------------------------------------------------------- */
923 char* pszOutputSRS =
924 psOptions->pszOutputSRS ? CPLStrdup(psOptions->pszOutputSRS) : nullptr;
925 for( int i = 0; i < nLayerCount; i++ )
926 {
927 OGRLayerH hLayer =
928 psOptions->papszLayers == nullptr
929 ? GDALDatasetGetLayer(hSrcDataset, 0)
930 : GDALDatasetGetLayerByName(hSrcDataset, psOptions->papszLayers[i]);
931 if( hLayer == nullptr )
932 {
933 CPLError(CE_Failure, CPLE_AppDefined, "Unable to find layer \"%s\", skipping.",
934 psOptions->papszLayers && psOptions->papszLayers[i] ?
935 psOptions->papszLayers[i] : "null" );
936 continue;
937 }
938
939 if( psOptions->pszWHERE )
940 {
941 if( OGR_L_SetAttributeFilter( hLayer, psOptions->pszWHERE ) != OGRERR_NONE )
942 break;
943 }
944
945 if ( psOptions->poSpatialFilter != nullptr )
946 OGR_L_SetSpatialFilter( hLayer, reinterpret_cast<OGRGeometryH>(psOptions->poSpatialFilter) );
947
948 // Fetch the first meaningful SRS definition
949 if ( !pszOutputSRS )
950 {
951 OGRSpatialReferenceH hSRS = OGR_L_GetSpatialRef( hLayer );
952 if ( hSRS )
953 OSRExportToWkt( hSRS, &pszOutputSRS );
954 }
955
956 eErr = ProcessLayer( hLayer, hDstDS, psOptions->poSpatialFilter, nXSize, nYSize,
957 i + 1 + nBands - nLayerCount,
958 bIsXExtentSet, bIsYExtentSet,
959 dfXMin, dfXMax, dfYMin, dfYMax, psOptions->pszBurnAttribute,
960 psOptions->dfIncreaseBurnValue, psOptions->dfMultiplyBurnValue,
961 psOptions->eOutputType, psOptions->eAlgorithm, psOptions->pOptions,
962 psOptions->bQuiet, psOptions->pfnProgress, psOptions->pProgressData );
963 if( eErr != CE_None )
964 break;
965 }
966
967 /* -------------------------------------------------------------------- */
968 /* Apply geotransformation matrix. */
969 /* -------------------------------------------------------------------- */
970 double adfGeoTransform[6] = {
971 dfXMin,
972 (dfXMax - dfXMin) / nXSize,
973 0.0,
974 dfYMin,
975 0.0,
976 (dfYMax - dfYMin) / nYSize };
977 GDALSetGeoTransform( hDstDS, adfGeoTransform );
978
979 /* -------------------------------------------------------------------- */
980 /* Apply SRS definition if set. */
981 /* -------------------------------------------------------------------- */
982 if ( pszOutputSRS )
983 {
984 GDALSetProjection( hDstDS, pszOutputSRS );
985 CPLFree(pszOutputSRS);
986 }
987
988 /* -------------------------------------------------------------------- */
989 /* End */
990 /* -------------------------------------------------------------------- */
991 GDALGridOptionsFree(psOptionsToFree);
992
993 if( eErr != CE_None )
994 {
995 GDALClose(hDstDS);
996 return nullptr;
997 }
998
999 return hDstDS;
1000 }
1001
1002 /************************************************************************/
1003 /* IsNumber() */
1004 /************************************************************************/
1005
IsNumber(const char * pszStr)1006 static bool IsNumber(const char* pszStr)
1007 {
1008 if (*pszStr == '-' || *pszStr == '+')
1009 pszStr ++;
1010 if (*pszStr == '.')
1011 pszStr ++;
1012 return *pszStr >= '0' && *pszStr <= '9';
1013 }
1014
1015 /************************************************************************/
1016 /* GDALGridOptionsNew() */
1017 /************************************************************************/
1018
1019 /**
1020 * Allocates a GDALGridOptions struct.
1021 *
1022 * @param papszArgv NULL terminated list of options (potentially including filename and open options too), or NULL.
1023 * The accepted options are the ones of the <a href="/programs/gdal_translate.html">gdal_translate</a> utility.
1024 * @param psOptionsForBinary (output) may be NULL (and should generally be NULL),
1025 * otherwise (gdal_translate_bin.cpp use case) must be allocated with
1026 * GDALGridOptionsForBinaryNew() prior to this function. Will be
1027 * filled with potentially present filename, open options,...
1028 * @return pointer to the allocated GDALGridOptions struct. Must be freed with GDALGridOptionsFree().
1029 *
1030 * @since GDAL 2.1
1031 */
1032
GDALGridOptionsNew(char ** papszArgv,GDALGridOptionsForBinary * psOptionsForBinary)1033 GDALGridOptions *GDALGridOptionsNew(char** papszArgv, GDALGridOptionsForBinary* psOptionsForBinary)
1034 {
1035 GDALGridOptions *psOptions =
1036 static_cast<GDALGridOptions *>(CPLCalloc(1, sizeof(GDALGridOptions)));
1037
1038 psOptions->pszFormat = nullptr;
1039 psOptions->bQuiet = true;
1040 psOptions->pfnProgress = GDALDummyProgress;
1041 psOptions->pProgressData = nullptr;
1042 psOptions->papszLayers = nullptr;
1043 psOptions->pszBurnAttribute = nullptr;
1044 psOptions->dfIncreaseBurnValue = 0.0;
1045 psOptions->dfMultiplyBurnValue = 1.0;
1046 psOptions->pszWHERE = nullptr;
1047 psOptions->pszSQL = nullptr;
1048 psOptions->eOutputType = GDT_Float64;
1049 psOptions->papszCreateOptions = nullptr;
1050 psOptions->nXSize = 0;
1051 psOptions->nYSize = 0;
1052 psOptions->dfXRes = 0;
1053 psOptions->dfYRes = 0;
1054 psOptions->dfXMin = 0.0;
1055 psOptions->dfXMax = 0.0;
1056 psOptions->dfYMin = 0.0;
1057 psOptions->dfYMax = 0.0;
1058 psOptions->bIsXExtentSet = false;
1059 psOptions->bIsYExtentSet = false;
1060 psOptions->eAlgorithm = GGA_InverseDistanceToAPower;
1061 psOptions->pOptions = nullptr;
1062 psOptions->pszOutputSRS = nullptr;
1063 psOptions->poSpatialFilter = nullptr;
1064 psOptions->poClipSrc = nullptr;
1065 psOptions->bClipSrc = false;
1066 psOptions->pszClipSrcDS = nullptr;
1067 psOptions->pszClipSrcSQL = nullptr;
1068 psOptions->pszClipSrcLayer = nullptr;
1069 psOptions->pszClipSrcWhere = nullptr;
1070 psOptions->bNoDataSet = false;
1071 psOptions->dfNoDataValue = 0;
1072
1073 ParseAlgorithmAndOptions( szAlgNameInvDist, &psOptions->eAlgorithm, &psOptions->pOptions );
1074
1075 bool bGotSourceFilename = false;
1076 bool bGotDestFilename = false;
1077 /* -------------------------------------------------------------------- */
1078 /* Handle command line arguments. */
1079 /* -------------------------------------------------------------------- */
1080 const int argc = CSLCount(papszArgv);
1081 for( int i = 0; i < argc && papszArgv != nullptr && papszArgv[i] != nullptr; i++ )
1082 {
1083 if( i < argc-1 && (EQUAL(papszArgv[i],"-of") || EQUAL(papszArgv[i],"-f")) )
1084 {
1085 ++i;
1086 CPLFree(psOptions->pszFormat);
1087 psOptions->pszFormat = CPLStrdup(papszArgv[i]);
1088 }
1089
1090 else if( EQUAL(papszArgv[i],"-q") || EQUAL(papszArgv[i],"-quiet") )
1091 {
1092 if( psOptionsForBinary )
1093 psOptionsForBinary->bQuiet = true;
1094 }
1095
1096 else if( EQUAL(papszArgv[i],"-ot") && papszArgv[i+1] )
1097 {
1098 for( int iType = 1; iType < GDT_TypeCount; iType++ )
1099 {
1100 if( GDALGetDataTypeName(static_cast<GDALDataType>(iType)) != nullptr
1101 && EQUAL(GDALGetDataTypeName(static_cast<GDALDataType>(iType)),
1102 papszArgv[i+1]) )
1103 {
1104 psOptions->eOutputType = static_cast<GDALDataType>(iType);
1105 }
1106 }
1107
1108 if( psOptions->eOutputType == GDT_Unknown )
1109 {
1110 CPLError(CE_Failure, CPLE_NotSupported,
1111 "Unknown output pixel type: %s.", papszArgv[i+1] );
1112 GDALGridOptionsFree(psOptions);
1113 return nullptr;
1114 }
1115 i++;
1116 }
1117
1118 else if( i+2 < argc && EQUAL(papszArgv[i],"-txe") )
1119 {
1120 psOptions->dfXMin = CPLAtof(papszArgv[++i]);
1121 psOptions->dfXMax = CPLAtof(papszArgv[++i]);
1122 psOptions->bIsXExtentSet = true;
1123 }
1124
1125 else if( i+2 < argc && EQUAL(papszArgv[i],"-tye") )
1126 {
1127 psOptions->dfYMin = CPLAtof(papszArgv[++i]);
1128 psOptions->dfYMax = CPLAtof(papszArgv[++i]);
1129 psOptions->bIsYExtentSet = true;
1130 }
1131
1132 else if( i+2 < argc && EQUAL(papszArgv[i],"-outsize") )
1133 {
1134 CPLAssert(papszArgv[i+1]);
1135 CPLAssert(papszArgv[i+2]);
1136 psOptions->nXSize = atoi(papszArgv[i+1]);
1137 psOptions->nYSize = atoi(papszArgv[i+2]);
1138 i += 2;
1139 }
1140
1141 else if( i+2 < argc && EQUAL(papszArgv[i],"-tr") )
1142 {
1143 psOptions->dfXRes = CPLAtofM(papszArgv[++i]);
1144 psOptions->dfYRes = CPLAtofM(papszArgv[++i]);
1145 if( psOptions->dfXRes <= 0 || psOptions->dfYRes <= 0 )
1146 {
1147 CPLError(CE_Failure, CPLE_IllegalArg, "Wrong value for -tr parameters.");
1148 GDALGridOptionsFree(psOptions);
1149 return nullptr;
1150 }
1151 }
1152
1153 else if( i+1 < argc && EQUAL(papszArgv[i],"-co") )
1154 {
1155 psOptions->papszCreateOptions = CSLAddString( psOptions->papszCreateOptions, papszArgv[++i] );
1156 }
1157
1158 else if( i+1 < argc && EQUAL(papszArgv[i],"-zfield") )
1159 {
1160 CPLFree(psOptions->pszBurnAttribute);
1161 psOptions->pszBurnAttribute = CPLStrdup(papszArgv[++i]);
1162 }
1163
1164 else if( i+1 < argc && EQUAL(papszArgv[i],"-z_increase") )
1165 {
1166 psOptions->dfIncreaseBurnValue = CPLAtof(papszArgv[++i]);
1167 }
1168
1169 else if( i+1 < argc && EQUAL(papszArgv[i],"-z_multiply") )
1170 {
1171 psOptions->dfMultiplyBurnValue = CPLAtof(papszArgv[++i]);
1172 }
1173
1174 else if( i+1 < argc && EQUAL(papszArgv[i],"-where") )
1175 {
1176 CPLFree(psOptions->pszWHERE);
1177 psOptions->pszWHERE = CPLStrdup(papszArgv[++i]);
1178 }
1179
1180 else if( i+1 < argc && EQUAL(papszArgv[i],"-l") )
1181 {
1182 psOptions->papszLayers = CSLAddString( psOptions->papszLayers, papszArgv[++i] );
1183 }
1184
1185 else if( i+1 < argc && EQUAL(papszArgv[i],"-sql") )
1186 {
1187 CPLFree(psOptions->pszSQL);
1188 psOptions->pszSQL = CPLStrdup(papszArgv[++i]);
1189 }
1190
1191 else if( i+4 < argc && EQUAL(papszArgv[i],"-spat") )
1192 {
1193 OGRLinearRing oRing;
1194
1195 oRing.addPoint( CPLAtof(papszArgv[i+1]), CPLAtof(papszArgv[i+2]) );
1196 oRing.addPoint( CPLAtof(papszArgv[i+1]), CPLAtof(papszArgv[i+4]) );
1197 oRing.addPoint( CPLAtof(papszArgv[i+3]), CPLAtof(papszArgv[i+4]) );
1198 oRing.addPoint( CPLAtof(papszArgv[i+3]), CPLAtof(papszArgv[i+2]) );
1199 oRing.addPoint( CPLAtof(papszArgv[i+1]), CPLAtof(papszArgv[i+2]) );
1200
1201 delete psOptions->poSpatialFilter;
1202 OGRPolygon *poPoly = new OGRPolygon();
1203 poPoly->addRing(&oRing);
1204 psOptions->poSpatialFilter = poPoly;
1205 i += 4;
1206 }
1207
1208 else if( EQUAL(papszArgv[i],"-clipsrc") )
1209 {
1210 if (i + 1 >= argc || papszArgv[i+1] == nullptr)
1211 {
1212 CPLError(CE_Failure, CPLE_IllegalArg, "%s option requires 1 or 4 arguments", papszArgv[i]);
1213 GDALGridOptionsFree(psOptions);
1214 return nullptr;
1215 }
1216
1217 VSIStatBufL sStat;
1218 psOptions->bClipSrc = true;
1219 if ( IsNumber(papszArgv[i+1])
1220 && papszArgv[i+2] != nullptr
1221 && papszArgv[i+3] != nullptr
1222 && papszArgv[i+4] != nullptr)
1223 {
1224 OGRLinearRing oRing;
1225
1226 oRing.addPoint( CPLAtof(papszArgv[i+1]), CPLAtof(papszArgv[i+2]) );
1227 oRing.addPoint( CPLAtof(papszArgv[i+1]), CPLAtof(papszArgv[i+4]) );
1228 oRing.addPoint( CPLAtof(papszArgv[i+3]), CPLAtof(papszArgv[i+4]) );
1229 oRing.addPoint( CPLAtof(papszArgv[i+3]), CPLAtof(papszArgv[i+2]) );
1230 oRing.addPoint( CPLAtof(papszArgv[i+1]), CPLAtof(papszArgv[i+2]) );
1231
1232 delete psOptions->poClipSrc;
1233 OGRPolygon *poPoly = static_cast<OGRPolygon *>(
1234 OGRGeometryFactory::createGeometry(wkbPolygon));
1235 poPoly->addRing(&oRing);
1236 psOptions->poClipSrc = poPoly;
1237 i += 4;
1238 }
1239 else if ((STARTS_WITH_CI(papszArgv[i+1], "POLYGON") ||
1240 STARTS_WITH_CI(papszArgv[i+1], "MULTIPOLYGON")) &&
1241 VSIStatL(papszArgv[i+1], &sStat) != 0)
1242 {
1243 delete psOptions->poClipSrc;
1244 OGRGeometryFactory::createFromWkt(papszArgv[i+1],
1245 nullptr, &psOptions->poClipSrc);
1246 if (psOptions->poClipSrc == nullptr)
1247 {
1248 CPLError(CE_Failure, CPLE_IllegalArg,
1249 "Invalid geometry. Must be a valid POLYGON or MULTIPOLYGON WKT");
1250 GDALGridOptionsFree(psOptions);
1251 return nullptr;
1252 }
1253 i ++;
1254 }
1255 else if (EQUAL(papszArgv[i+1], "spat_extent") )
1256 {
1257 i ++;
1258 }
1259 else
1260 {
1261 CPLFree(psOptions->pszClipSrcDS);
1262 psOptions->pszClipSrcDS = CPLStrdup(papszArgv[i+1]);
1263 i ++;
1264 }
1265 }
1266 else if( i+1 < argc && EQUAL(papszArgv[i],"-clipsrcsql") )
1267 {
1268 CPLFree(psOptions->pszClipSrcSQL);
1269 psOptions->pszClipSrcSQL = CPLStrdup(papszArgv[i+1]);
1270 i ++;
1271 }
1272 else if( i+1 < argc && EQUAL(papszArgv[i],"-clipsrclayer") )
1273 {
1274 CPLFree(psOptions->pszClipSrcLayer);
1275 psOptions->pszClipSrcLayer = CPLStrdup(papszArgv[i+1]);
1276 i ++;
1277 }
1278 else if( i+1 < argc && EQUAL(papszArgv[i],"-clipsrcwhere") )
1279 {
1280 CPLFree(psOptions->pszClipSrcWhere);
1281 psOptions->pszClipSrcWhere = CPLStrdup(papszArgv[i+1]);
1282 i ++;
1283 }
1284
1285 else if( i+1 < argc && EQUAL(papszArgv[i],"-a_srs") )
1286 {
1287 OGRSpatialReference oOutputSRS;
1288
1289 if( oOutputSRS.SetFromUserInput( papszArgv[i+1] ) != OGRERR_NONE )
1290 {
1291 CPLError(CE_Failure, CPLE_AppDefined,
1292 "Failed to process SRS definition: %s",
1293 papszArgv[i+1] );
1294 GDALGridOptionsFree(psOptions);
1295 return nullptr;
1296 }
1297
1298 CPLFree(psOptions->pszOutputSRS);
1299 oOutputSRS.exportToWkt( &(psOptions->pszOutputSRS) );
1300 i++;
1301 }
1302
1303 else if( i+1 < argc && EQUAL(papszArgv[i],"-a") )
1304 {
1305 const char* pszAlgorithm = papszArgv[++i];
1306 CPLFree(psOptions->pOptions);
1307 if ( ParseAlgorithmAndOptions( pszAlgorithm, &psOptions->eAlgorithm, &psOptions->pOptions )
1308 != CE_None )
1309 {
1310 CPLError(CE_Failure, CPLE_AppDefined,
1311 "Failed to process algorithm name and parameters" );
1312 GDALGridOptionsFree(psOptions);
1313 return nullptr;
1314 }
1315
1316 char **papszParams = CSLTokenizeString2( pszAlgorithm, ":", FALSE );
1317 const char* pszNoDataValue = CSLFetchNameValue( papszParams, "nodata" );
1318 if( pszNoDataValue != nullptr )
1319 {
1320 psOptions->bNoDataSet = true;
1321 psOptions->dfNoDataValue = CPLAtofM(pszNoDataValue);
1322 }
1323 CSLDestroy(papszParams);
1324 }
1325 else if( papszArgv[i][0] == '-' )
1326 {
1327 CPLError(CE_Failure, CPLE_NotSupported,
1328 "Unknown option name '%s'", papszArgv[i]);
1329 GDALGridOptionsFree(psOptions);
1330 return nullptr;
1331 }
1332 else if( !bGotSourceFilename )
1333 {
1334 bGotSourceFilename = true;
1335 if( psOptionsForBinary )
1336 psOptionsForBinary->pszSource = CPLStrdup(papszArgv[i]);
1337 }
1338 else if( !bGotDestFilename )
1339 {
1340 bGotDestFilename = true;
1341 if( psOptionsForBinary )
1342 psOptionsForBinary->pszDest = CPLStrdup(papszArgv[i]);
1343 }
1344 else
1345 {
1346 CPLError(CE_Failure, CPLE_NotSupported,
1347 "Too many command options '%s'", papszArgv[i]);
1348 GDALGridOptionsFree(psOptions);
1349 return nullptr;
1350 }
1351 }
1352
1353 if ( psOptions->bClipSrc && psOptions->pszClipSrcDS != nullptr )
1354 {
1355 psOptions->poClipSrc = LoadGeometry( psOptions->pszClipSrcDS, psOptions->pszClipSrcSQL,
1356 psOptions->pszClipSrcLayer, psOptions->pszClipSrcWhere );
1357 if ( psOptions->poClipSrc == nullptr )
1358 {
1359 CPLError(CE_Failure, CPLE_AppDefined, "Cannot load source clip geometry.");
1360 GDALGridOptionsFree(psOptions);
1361 return nullptr;
1362 }
1363 }
1364 else if ( psOptions->bClipSrc && psOptions->poClipSrc == nullptr && !psOptions->poSpatialFilter )
1365 {
1366 CPLError(CE_Failure, CPLE_AppDefined, "-clipsrc must be used with -spat option or \n"
1367 "a bounding box, WKT string or datasource must be "
1368 "specified.");
1369 GDALGridOptionsFree(psOptions);
1370 return nullptr;
1371 }
1372
1373 if ( psOptions->poSpatialFilter )
1374 {
1375 if ( psOptions->poClipSrc )
1376 {
1377 OGRGeometry *poTemp = psOptions->poSpatialFilter->Intersection( psOptions->poClipSrc );
1378 if ( poTemp )
1379 {
1380 delete psOptions->poSpatialFilter;
1381 psOptions->poSpatialFilter = poTemp;
1382 }
1383
1384 delete psOptions->poClipSrc;
1385 psOptions->poClipSrc = nullptr;
1386 }
1387 }
1388 else
1389 {
1390 if ( psOptions->poClipSrc )
1391 {
1392 psOptions->poSpatialFilter = psOptions->poClipSrc;
1393 psOptions->poClipSrc = nullptr;
1394 }
1395 }
1396
1397 if( psOptionsForBinary )
1398 {
1399 if( psOptions->pszFormat )
1400 {
1401 psOptionsForBinary->pszFormat = CPLStrdup(psOptions->pszFormat);
1402 }
1403 }
1404
1405 return psOptions;
1406 }
1407
1408 /************************************************************************/
1409 /* GDALGridOptionsFree() */
1410 /************************************************************************/
1411
1412 /**
1413 * Frees the GDALGridOptions struct.
1414 *
1415 * @param psOptions the options struct for GDALGrid().
1416 *
1417 * @since GDAL 2.1
1418 */
1419
GDALGridOptionsFree(GDALGridOptions * psOptions)1420 void GDALGridOptionsFree(GDALGridOptions *psOptions)
1421 {
1422 if( psOptions == nullptr )
1423 return;
1424
1425 CPLFree(psOptions->pszFormat);
1426 CSLDestroy(psOptions->papszLayers);
1427 CPLFree(psOptions->pszBurnAttribute);
1428 CPLFree(psOptions->pszWHERE);
1429 CPLFree(psOptions->pszSQL);
1430 CSLDestroy(psOptions->papszCreateOptions);
1431 CPLFree(psOptions->pOptions);
1432 CPLFree(psOptions->pszOutputSRS);
1433 delete psOptions->poSpatialFilter;
1434 delete psOptions->poClipSrc;
1435 CPLFree(psOptions->pszClipSrcDS);
1436 CPLFree(psOptions->pszClipSrcSQL);
1437 CPLFree(psOptions->pszClipSrcLayer);
1438 CPLFree(psOptions->pszClipSrcWhere);
1439 CPLFree(psOptions);
1440 }
1441
1442 /************************************************************************/
1443 /* GDALGridOptionsSetProgress() */
1444 /************************************************************************/
1445
1446 /**
1447 * Set a progress function.
1448 *
1449 * @param psOptions the options struct for GDALGrid().
1450 * @param pfnProgress the progress callback.
1451 * @param pProgressData the user data for the progress callback.
1452 *
1453 * @since GDAL 2.1
1454 */
1455
GDALGridOptionsSetProgress(GDALGridOptions * psOptions,GDALProgressFunc pfnProgress,void * pProgressData)1456 void GDALGridOptionsSetProgress( GDALGridOptions *psOptions,
1457 GDALProgressFunc pfnProgress,
1458 void *pProgressData )
1459 {
1460 psOptions->pfnProgress = pfnProgress;
1461 psOptions->pProgressData = pProgressData;
1462 if( pfnProgress == GDALTermProgress )
1463 psOptions->bQuiet = false;
1464 }
1465