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