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