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 106 : GDALGridOptions()
97 106 : {
98 106 : void *l_pOptions = nullptr;
99 106 : GDALGridParseAlgorithmAndOptions(szAlgNameInvDist, &eAlgorithm,
100 : &l_pOptions);
101 106 : pOptions.reset(l_pOptions);
102 106 : }
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 64928 : void visit(const OGRPoint *p) override
278 : {
279 64928 : if (poClipSrc && !p->Within(poClipSrc))
280 20 : return;
281 :
282 64908 : if (iBurnField < 0 && std::isnan(p->getZ()))
283 1 : return;
284 :
285 64907 : adfX.push_back(p->getX());
286 64907 : adfY.push_back(p->getY());
287 64907 : if (iBurnField < 0)
288 129806 : adfZ.push_back((p->getZ() + dfIncreaseBurnValue) *
289 64903 : dfMultiplyBurnValue);
290 : else
291 4 : adfZ.push_back((dfBurnValue + dfIncreaseBurnValue) *
292 4 : 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 103 : 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 103 : int iBurnField = -1;
319 :
320 103 : if (!osBurnAttribute.empty())
321 : {
322 : iBurnField =
323 1 : poSrcLayer->GetLayerDefn()->GetFieldIndex(osBurnAttribute.c_str());
324 1 : 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 206 : GDALGridGeometryVisitor oVisitor;
337 103 : oVisitor.poClipSrc = poClipSrc;
338 103 : oVisitor.iBurnField = iBurnField;
339 103 : oVisitor.dfIncreaseBurnValue = dfIncreaseBurnValue;
340 103 : oVisitor.dfMultiplyBurnValue = dfMultiplyBurnValue;
341 :
342 64882 : for (auto &&poFeat : poSrcLayer)
343 : {
344 64779 : const OGRGeometry *poGeom = poFeat->GetGeometryRef();
345 64779 : if (poGeom)
346 : {
347 64779 : if (iBurnField >= 0)
348 : {
349 5 : if (!poFeat->IsFieldSetAndNotNull(iBurnField))
350 : {
351 1 : continue;
352 : }
353 4 : oVisitor.dfBurnValue = poFeat->GetFieldAsDouble(iBurnField);
354 : }
355 :
356 64778 : poGeom->accept(&oVisitor);
357 : }
358 : }
359 :
360 103 : 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 103 : if (!bIsXExtentSet || !bIsYExtentSet)
371 : {
372 0 : OGREnvelope sEnvelope;
373 0 : if (poSrcLayer->GetExtent(&sEnvelope, TRUE) == OGRERR_FAILURE)
374 : {
375 0 : return CE_Failure;
376 : }
377 :
378 0 : if (!bIsXExtentSet)
379 : {
380 0 : dfXMin = sEnvelope.MinX;
381 0 : dfXMax = sEnvelope.MaxX;
382 0 : bIsXExtentSet = true;
383 : }
384 :
385 0 : if (!bIsYExtentSet)
386 : {
387 0 : dfYMin = sEnvelope.MinY;
388 0 : dfYMax = sEnvelope.MaxY;
389 0 : bIsYExtentSet = true;
390 : }
391 : }
392 :
393 : // Produce north-up images
394 103 : if (dfYMin < dfYMax)
395 51 : std::swap(dfYMin, dfYMax);
396 :
397 : /* -------------------------------------------------------------------- */
398 : /* Perform gridding. */
399 : /* -------------------------------------------------------------------- */
400 :
401 103 : const double dfDeltaX = (dfXMax - dfXMin) / nXSize;
402 103 : const double dfDeltaY = (dfYMax - dfYMin) / nYSize;
403 :
404 103 : 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 103 : GDALRasterBand *poBand = poDstDS->GetRasterBand(nBand);
418 :
419 103 : int nBlockXSize = 0;
420 103 : int nBlockYSize = 0;
421 103 : const int nDataTypeSize = GDALGetDataTypeSizeBytes(eType);
422 :
423 : // Try to grow the work buffer up to 16 MB if it is smaller
424 103 : poBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
425 103 : if (nXSize == 0 || nYSize == 0 || nBlockXSize == 0 || nBlockYSize == 0)
426 0 : return CE_Failure;
427 :
428 103 : const int nDesiredBufferSize = 16 * 1024 * 1024;
429 103 : 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 103 : else if (nBlockXSize == nXSize && nBlockYSize < nYSize &&
439 7 : nBlockYSize < nDesiredBufferSize / (nXSize * nDataTypeSize))
440 : {
441 7 : const int nNewBlockYSize =
442 7 : nDesiredBufferSize / (nXSize * nDataTypeSize);
443 7 : nBlockYSize = (nNewBlockYSize / nBlockYSize) * nBlockYSize;
444 7 : if (nBlockYSize > nYSize)
445 7 : nBlockYSize = nYSize;
446 : }
447 103 : CPLDebug("GDAL_GRID", "Work buffer: %d * %d", nBlockXSize, nBlockYSize);
448 :
449 : std::unique_ptr<void, VSIFreeReleaser> pData(
450 206 : VSIMalloc3(nBlockXSize, nBlockYSize, nDataTypeSize));
451 103 : if (!pData)
452 : {
453 0 : CPLError(CE_Failure, CPLE_OutOfMemory, "Cannot allocate work buffer");
454 0 : return CE_Failure;
455 : }
456 :
457 103 : GIntBig nBlock = 0;
458 103 : const double dfBlockCount =
459 103 : static_cast<double>(DIV_ROUND_UP(nXSize, nBlockXSize)) *
460 103 : DIV_ROUND_UP(nYSize, nBlockYSize);
461 :
462 : struct GDALGridContextReleaser
463 : {
464 103 : void operator()(GDALGridContext *psContext)
465 : {
466 103 : GDALGridContextFree(psContext);
467 103 : }
468 : };
469 :
470 : std::unique_ptr<GDALGridContext, GDALGridContextReleaser> psContext(
471 : GDALGridContextCreate(eAlgorithm, pOptions,
472 103 : static_cast<int>(oVisitor.adfX.size()),
473 103 : &(oVisitor.adfX[0]), &(oVisitor.adfY[0]),
474 309 : &(oVisitor.adfZ[0]), TRUE));
475 103 : if (!psContext)
476 : {
477 0 : return CE_Failure;
478 : }
479 :
480 103 : CPLErr eErr = CE_None;
481 206 : for (int nYOffset = 0; nYOffset < nYSize && eErr == CE_None;
482 103 : nYOffset += nBlockYSize)
483 : {
484 206 : for (int nXOffset = 0; nXOffset < nXSize && eErr == CE_None;
485 103 : nXOffset += nBlockXSize)
486 : {
487 : std::unique_ptr<void, GDALScaledProgressReleaser> pScaledProgress(
488 : GDALCreateScaledProgress(
489 103 : static_cast<double>(nBlock) / dfBlockCount,
490 103 : static_cast<double>(nBlock + 1) / dfBlockCount, pfnProgress,
491 206 : pProgressData));
492 103 : nBlock++;
493 :
494 103 : int nXRequest = nBlockXSize;
495 103 : if (nXOffset > nXSize - nXRequest)
496 2 : nXRequest = nXSize - nXOffset;
497 :
498 103 : int nYRequest = nBlockYSize;
499 103 : if (nYOffset > nYSize - nYRequest)
500 2 : nYRequest = nYSize - nYOffset;
501 :
502 206 : eErr = GDALGridContextProcess(
503 103 : psContext.get(), dfXMin + dfDeltaX * nXOffset,
504 103 : dfXMin + dfDeltaX * (nXOffset + nXRequest),
505 103 : dfYMin + dfDeltaY * nYOffset,
506 103 : dfYMin + dfDeltaY * (nYOffset + nYRequest), nXRequest,
507 : nYRequest, eType, pData.get(), GDALScaledProgress,
508 : pScaledProgress.get());
509 :
510 103 : if (eErr == CE_None)
511 103 : eErr = poBand->RasterIO(GF_Write, nXOffset, nYOffset, nXRequest,
512 : nYRequest, pData.get(), nXRequest,
513 : nYRequest, eType, 0, 0, nullptr);
514 : }
515 : }
516 103 : if (eErr == CE_None && pfnProgress)
517 103 : pfnProgress(1.0, "", pProgressData);
518 :
519 103 : 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 106 : GDALDatasetH GDALGrid(const char *pszDest, GDALDatasetH hSrcDataset,
629 : const GDALGridOptions *psOptionsIn, int *pbUsageError)
630 :
631 : {
632 106 : 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 106 : 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 106 : std::unique_ptr<GDALGridOptions> psOptionsToFree;
650 106 : const GDALGridOptions *psOptions = psOptionsIn;
651 106 : if (psOptions == nullptr)
652 : {
653 0 : psOptionsToFree = std::make_unique<GDALGridOptions>();
654 0 : psOptions = psOptionsToFree.get();
655 : }
656 :
657 106 : GDALDataset *poSrcDS = GDALDataset::FromHandle(hSrcDataset);
658 :
659 153 : if (psOptions->osSQL.empty() && psOptions->aosLayers.empty() &&
660 47 : 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 106 : if ((psOptions->nXSize != 0 || psOptions->nYSize != 0) &&
671 105 : (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 212 : std::string osFormat;
682 106 : if (psOptions->osFormat.empty())
683 : {
684 54 : osFormat = GetOutputDriverForRaster(pszDest);
685 54 : if (osFormat.empty())
686 : {
687 0 : return nullptr;
688 : }
689 : }
690 : else
691 : {
692 52 : osFormat = psOptions->osFormat;
693 : }
694 :
695 106 : GDALDriverH hDriver = GDALGetDriverByName(osFormat.c_str());
696 106 : 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 configured and "
701 : "support output:\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 106 : int nLayerCount = psOptions->aosLayers.size();
725 106 : if (nLayerCount == 0 && psOptions->osSQL.empty())
726 47 : nLayerCount = 1; /* due to above check */
727 :
728 106 : int nBands = nLayerCount;
729 :
730 106 : if (!psOptions->osSQL.empty())
731 4 : nBands++;
732 :
733 : int nXSize;
734 : int nYSize;
735 106 : if (psOptions->dfXRes != 0 && psOptions->dfYRes != 0)
736 : {
737 1 : if ((psOptions->dfXMax == psOptions->dfXMin) ||
738 1 : (psOptions->dfYMax == psOptions->dfYMin))
739 : {
740 0 : CPLError(CE_Failure, CPLE_IllegalArg,
741 : "Invalid txe or tye parameters detected. Please check "
742 : "your -txe or -tye argument.");
743 :
744 0 : if (pbUsageError)
745 0 : *pbUsageError = TRUE;
746 0 : return nullptr;
747 : }
748 :
749 1 : double dfXSize = (std::fabs(psOptions->dfXMax - psOptions->dfXMin) +
750 1 : (psOptions->dfXRes / 2.0)) /
751 1 : psOptions->dfXRes;
752 1 : double dfYSize = (std::fabs(psOptions->dfYMax - psOptions->dfYMin) +
753 1 : (psOptions->dfYRes / 2.0)) /
754 1 : psOptions->dfYRes;
755 :
756 1 : if (dfXSize >= 1 && dfXSize <= INT_MAX && dfYSize >= 1 &&
757 : dfYSize <= INT_MAX)
758 : {
759 1 : nXSize = static_cast<int>(dfXSize);
760 1 : nYSize = static_cast<int>(dfYSize);
761 : }
762 : else
763 : {
764 0 : CPLError(
765 : CE_Failure, CPLE_IllegalArg,
766 : "Invalid output size detected. Please check your -tr argument");
767 :
768 0 : if (pbUsageError)
769 0 : *pbUsageError = TRUE;
770 0 : return nullptr;
771 1 : }
772 : }
773 : else
774 : {
775 : // FIXME
776 105 : nXSize = psOptions->nXSize;
777 105 : if (nXSize == 0)
778 0 : nXSize = 256;
779 105 : nYSize = psOptions->nYSize;
780 105 : if (nYSize == 0)
781 0 : nYSize = 256;
782 : }
783 :
784 : std::unique_ptr<GDALDataset> poDstDS(GDALDataset::FromHandle(GDALCreate(
785 106 : hDriver, pszDest, nXSize, nYSize, nBands, psOptions->eOutputType,
786 212 : psOptions->aosCreateOptions.List())));
787 106 : if (!poDstDS)
788 : {
789 0 : return nullptr;
790 : }
791 :
792 106 : if (psOptions->bNoDataSet)
793 : {
794 104 : for (int i = 1; i <= nBands; i++)
795 : {
796 52 : poDstDS->GetRasterBand(i)->SetNoDataValue(psOptions->dfNoDataValue);
797 : }
798 : }
799 :
800 106 : double dfXMin = psOptions->dfXMin;
801 106 : double dfYMin = psOptions->dfYMin;
802 106 : double dfXMax = psOptions->dfXMax;
803 106 : double dfYMax = psOptions->dfYMax;
804 106 : bool bIsXExtentSet = psOptions->bIsXExtentSet;
805 106 : bool bIsYExtentSet = psOptions->bIsYExtentSet;
806 106 : CPLErr eErr = CE_None;
807 :
808 : /* -------------------------------------------------------------------- */
809 : /* Process SQL request. */
810 : /* -------------------------------------------------------------------- */
811 :
812 106 : if (!psOptions->osSQL.empty())
813 : {
814 : OGRLayer *poLayer =
815 4 : poSrcDS->ExecuteSQL(psOptions->osSQL.c_str(),
816 4 : psOptions->poSpatialFilter.get(), nullptr);
817 4 : if (poLayer == nullptr)
818 : {
819 1 : return nullptr;
820 : }
821 :
822 : // Custom layer will be rasterized in the first band.
823 3 : eErr = ProcessLayer(
824 3 : poLayer, poDstDS.get(), psOptions->poSpatialFilter.get(), nXSize,
825 : nYSize, 1, bIsXExtentSet, bIsYExtentSet, dfXMin, dfXMax, dfYMin,
826 3 : dfYMax, psOptions->osBurnAttribute, psOptions->dfIncreaseBurnValue,
827 3 : psOptions->dfMultiplyBurnValue, psOptions->eOutputType,
828 3 : psOptions->eAlgorithm, psOptions->pOptions.get(), psOptions->bQuiet,
829 3 : psOptions->pfnProgress, psOptions->pProgressData);
830 :
831 3 : poSrcDS->ReleaseResultSet(poLayer);
832 : }
833 :
834 : /* -------------------------------------------------------------------- */
835 : /* Process each layer. */
836 : /* -------------------------------------------------------------------- */
837 210 : std::string osOutputSRS(psOptions->osOutputSRS);
838 205 : for (int i = 0; i < nLayerCount; i++)
839 : {
840 102 : auto poLayer = psOptions->aosLayers.empty()
841 102 : ? poSrcDS->GetLayer(0)
842 55 : : poSrcDS->GetLayerByName(psOptions->aosLayers[i]);
843 102 : if (!poLayer)
844 : {
845 1 : CPLError(CE_Failure, CPLE_AppDefined,
846 : "Unable to find layer \"%s\".",
847 1 : !psOptions->aosLayers.empty() && psOptions->aosLayers[i]
848 1 : ? psOptions->aosLayers[i]
849 : : "null");
850 1 : eErr = CE_Failure;
851 1 : break;
852 : }
853 :
854 101 : if (!psOptions->osWHERE.empty())
855 : {
856 2 : if (poLayer->SetAttributeFilter(psOptions->osWHERE.c_str()) !=
857 : OGRERR_NONE)
858 : {
859 1 : eErr = CE_Failure;
860 1 : break;
861 : }
862 : }
863 :
864 100 : if (psOptions->poSpatialFilter)
865 2 : poLayer->SetSpatialFilter(psOptions->poSpatialFilter.get());
866 :
867 : // Fetch the first meaningful SRS definition
868 100 : if (osOutputSRS.empty())
869 : {
870 100 : auto poSRS = poLayer->GetSpatialRef();
871 100 : if (poSRS)
872 91 : osOutputSRS = poSRS->exportToWkt();
873 : }
874 :
875 100 : eErr = ProcessLayer(
876 100 : poLayer, poDstDS.get(), psOptions->poSpatialFilter.get(), nXSize,
877 100 : nYSize, i + 1 + nBands - nLayerCount, bIsXExtentSet, bIsYExtentSet,
878 100 : dfXMin, dfXMax, dfYMin, dfYMax, psOptions->osBurnAttribute,
879 100 : psOptions->dfIncreaseBurnValue, psOptions->dfMultiplyBurnValue,
880 100 : psOptions->eOutputType, psOptions->eAlgorithm,
881 100 : psOptions->pOptions.get(), psOptions->bQuiet,
882 100 : psOptions->pfnProgress, psOptions->pProgressData);
883 100 : if (eErr != CE_None)
884 0 : break;
885 : }
886 :
887 : /* -------------------------------------------------------------------- */
888 : /* Apply geotransformation matrix. */
889 : /* -------------------------------------------------------------------- */
890 105 : double adfGeoTransform[6] = {dfXMin, (dfXMax - dfXMin) / nXSize,
891 : 0.0, dfYMin,
892 105 : 0.0, (dfYMax - dfYMin) / nYSize};
893 105 : poDstDS->SetGeoTransform(adfGeoTransform);
894 :
895 : /* -------------------------------------------------------------------- */
896 : /* Apply SRS definition if set. */
897 : /* -------------------------------------------------------------------- */
898 105 : if (!osOutputSRS.empty())
899 : {
900 91 : poDstDS->SetProjection(osOutputSRS.c_str());
901 : }
902 :
903 : /* -------------------------------------------------------------------- */
904 : /* End */
905 : /* -------------------------------------------------------------------- */
906 :
907 105 : if (eErr != CE_None)
908 : {
909 2 : return nullptr;
910 : }
911 :
912 103 : return GDALDataset::ToHandle(poDstDS.release());
913 : }
914 :
915 : /************************************************************************/
916 : /* GDALGridOptionsGetParser() */
917 : /************************************************************************/
918 :
919 : /*! @cond Doxygen_Suppress */
920 :
921 : static std::unique_ptr<GDALArgumentParser>
922 106 : GDALGridOptionsGetParser(GDALGridOptions *psOptions,
923 : GDALGridOptionsForBinary *psOptionsForBinary,
924 : int nCountClipSrc)
925 : {
926 : auto argParser = std::make_unique<GDALArgumentParser>(
927 106 : "gdal_grid", /* bForBinary=*/psOptionsForBinary != nullptr);
928 :
929 106 : argParser->add_description(
930 : _("Creates a regular grid (raster) from the scattered data read from a "
931 106 : "vector datasource."));
932 :
933 106 : argParser->add_epilog(_(
934 : "Available algorithms and parameters with their defaults:\n"
935 : " Inverse distance to a power (default)\n"
936 : " "
937 : "invdist:power=2.0:smoothing=0.0:radius1=0.0:radius2=0.0:angle=0.0:max_"
938 : "points=0:min_points=0:nodata=0.0\n"
939 : " Inverse distance to a power with nearest neighbor search\n"
940 : " "
941 : "invdistnn:power=2.0:radius=1.0:max_points=12:min_points=0:nodata=0\n"
942 : " Moving average\n"
943 : " "
944 : "average:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0\n"
945 : " Nearest neighbor\n"
946 : " nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0\n"
947 : " Various data metrics\n"
948 : " <metric "
949 : "name>:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0\n"
950 : " possible metrics are:\n"
951 : " minimum\n"
952 : " maximum\n"
953 : " range\n"
954 : " count\n"
955 : " average_distance\n"
956 : " average_distance_pts\n"
957 : " Linear\n"
958 : " linear:radius=-1.0:nodata=0.0\n"
959 : "\n"
960 106 : "For more details, consult https://gdal.org/programs/gdal_grid.html"));
961 :
962 : argParser->add_quiet_argument(
963 106 : psOptionsForBinary ? &psOptionsForBinary->bQuiet : nullptr);
964 :
965 106 : argParser->add_output_format_argument(psOptions->osFormat);
966 :
967 106 : argParser->add_output_type_argument(psOptions->eOutputType);
968 :
969 106 : argParser->add_argument("-txe")
970 212 : .metavar("<xmin> <xmax>")
971 106 : .nargs(2)
972 106 : .scan<'g', double>()
973 106 : .help(_("Set georeferenced X extents of output file to be created."));
974 :
975 106 : argParser->add_argument("-tye")
976 212 : .metavar("<ymin> <ymax>")
977 106 : .nargs(2)
978 106 : .scan<'g', double>()
979 106 : .help(_("Set georeferenced Y extents of output file to be created."));
980 :
981 106 : argParser->add_argument("-outsize")
982 212 : .metavar("<xsize> <ysize>")
983 106 : .nargs(2)
984 106 : .scan<'i', int>()
985 106 : .help(_("Set the size of the output file."));
986 :
987 106 : argParser->add_argument("-tr")
988 212 : .metavar("<xres> <yes>")
989 106 : .nargs(2)
990 106 : .scan<'g', double>()
991 106 : .help(_("Set target resolution."));
992 :
993 106 : argParser->add_creation_options_argument(psOptions->aosCreateOptions);
994 :
995 106 : argParser->add_argument("-zfield")
996 212 : .metavar("<field_name>")
997 106 : .store_into(psOptions->osBurnAttribute)
998 106 : .help(_("Field name from which to get Z values."));
999 :
1000 106 : argParser->add_argument("-z_increase")
1001 212 : .metavar("<increase_value>")
1002 106 : .store_into(psOptions->dfIncreaseBurnValue)
1003 : .help(_("Addition to the attribute field on the features to be used to "
1004 106 : "get a Z value from."));
1005 :
1006 106 : argParser->add_argument("-z_multiply")
1007 212 : .metavar("<multiply_value>")
1008 106 : .store_into(psOptions->dfMultiplyBurnValue)
1009 106 : .help(_("Multiplication ratio for the Z field.."));
1010 :
1011 106 : argParser->add_argument("-where")
1012 212 : .metavar("<expression>")
1013 106 : .store_into(psOptions->osWHERE)
1014 : .help(_("Query expression to be applied to select features to process "
1015 106 : "from the input layer(s)."));
1016 :
1017 106 : argParser->add_argument("-l")
1018 212 : .metavar("<layer_name>")
1019 106 : .append()
1020 55 : .action([psOptions](const std::string &s)
1021 161 : { psOptions->aosLayers.AddString(s.c_str()); })
1022 : .help(_("Layer(s) from the datasource that will be used for input "
1023 106 : "features."));
1024 :
1025 106 : argParser->add_argument("-sql")
1026 212 : .metavar("<select_statement>")
1027 106 : .store_into(psOptions->osSQL)
1028 : .help(_("SQL statement to be evaluated to produce a layer of features "
1029 106 : "to be processed."));
1030 :
1031 106 : argParser->add_argument("-spat")
1032 212 : .metavar("<xmin> <ymin> <xmax> <ymax>")
1033 106 : .nargs(4)
1034 106 : .scan<'g', double>()
1035 : .help(_("The area of interest. Only features within the rectangle will "
1036 106 : "be reported."));
1037 :
1038 106 : argParser->add_argument("-clipsrc")
1039 106 : .nargs(nCountClipSrc)
1040 212 : .metavar("[<xmin> <ymin> <xmax> <ymax>]|<WKT>|<datasource>|spat_extent")
1041 106 : .help(_("Clip geometries (in source SRS)."));
1042 :
1043 106 : argParser->add_argument("-clipsrcsql")
1044 212 : .metavar("<sql_statement>")
1045 106 : .store_into(psOptions->osClipSrcSQL)
1046 : .help(_("Select desired geometries from the source clip datasource "
1047 106 : "using an SQL query."));
1048 :
1049 106 : argParser->add_argument("-clipsrclayer")
1050 212 : .metavar("<layername>")
1051 106 : .store_into(psOptions->osClipSrcLayer)
1052 106 : .help(_("Select the named layer from the source clip datasource."));
1053 :
1054 106 : argParser->add_argument("-clipsrcwhere")
1055 212 : .metavar("<expression>")
1056 106 : .store_into(psOptions->osClipSrcWhere)
1057 : .help(_("Restrict desired geometries from the source clip layer based "
1058 106 : "on an attribute query."));
1059 :
1060 106 : argParser->add_argument("-a_srs")
1061 212 : .metavar("<srs_def>")
1062 : .action(
1063 0 : [psOptions](const std::string &osOutputSRSDef)
1064 : {
1065 0 : OGRSpatialReference oOutputSRS;
1066 :
1067 0 : if (oOutputSRS.SetFromUserInput(osOutputSRSDef.c_str()) !=
1068 : OGRERR_NONE)
1069 : {
1070 : throw std::invalid_argument(
1071 0 : std::string("Failed to process SRS definition: ")
1072 0 : .append(osOutputSRSDef));
1073 : }
1074 :
1075 0 : char *pszWKT = nullptr;
1076 0 : oOutputSRS.exportToWkt(&pszWKT);
1077 0 : if (pszWKT)
1078 0 : psOptions->osOutputSRS = pszWKT;
1079 0 : CPLFree(pszWKT);
1080 106 : })
1081 106 : .help(_("Assign an output SRS, but without reprojecting."));
1082 :
1083 106 : argParser->add_argument("-a")
1084 212 : .metavar("<algorithm>[[:<parameter1>=<value1>]...]")
1085 : .action(
1086 352 : [psOptions](const std::string &s)
1087 : {
1088 100 : const char *pszAlgorithm = s.c_str();
1089 100 : void *pOptions = nullptr;
1090 100 : if (GDALGridParseAlgorithmAndOptions(pszAlgorithm,
1091 : &psOptions->eAlgorithm,
1092 100 : &pOptions) != CE_None)
1093 : {
1094 : throw std::invalid_argument(
1095 0 : "Failed to process algorithm name and parameters");
1096 : }
1097 100 : psOptions->pOptions.reset(pOptions);
1098 :
1099 : const CPLStringList aosParams(
1100 200 : CSLTokenizeString2(pszAlgorithm, ":", FALSE));
1101 100 : const char *pszNoDataValue = aosParams.FetchNameValue("nodata");
1102 100 : if (pszNoDataValue != nullptr)
1103 : {
1104 52 : psOptions->bNoDataSet = true;
1105 52 : psOptions->dfNoDataValue = CPLAtofM(pszNoDataValue);
1106 : }
1107 206 : })
1108 : .help(_("Set the interpolation algorithm or data metric name and "
1109 106 : "(optionally) its parameters."));
1110 :
1111 106 : if (psOptionsForBinary)
1112 : {
1113 : argParser->add_open_options_argument(
1114 54 : &(psOptionsForBinary->aosOpenOptions));
1115 : }
1116 :
1117 106 : if (psOptionsForBinary)
1118 : {
1119 54 : argParser->add_argument("src_dataset_name")
1120 108 : .metavar("<src_dataset_name>")
1121 54 : .store_into(psOptionsForBinary->osSource)
1122 54 : .help(_("Input dataset."));
1123 :
1124 54 : argParser->add_argument("dst_dataset_name")
1125 108 : .metavar("<dst_dataset_name>")
1126 54 : .store_into(psOptionsForBinary->osDest)
1127 54 : .help(_("Output dataset."));
1128 : }
1129 :
1130 106 : return argParser;
1131 : }
1132 :
1133 : /*! @endcond */
1134 :
1135 : /************************************************************************/
1136 : /* GDALGridGetParserUsage() */
1137 : /************************************************************************/
1138 :
1139 0 : std::string GDALGridGetParserUsage()
1140 : {
1141 : try
1142 : {
1143 0 : GDALGridOptions sOptions;
1144 0 : GDALGridOptionsForBinary sOptionsForBinary;
1145 : auto argParser =
1146 0 : GDALGridOptionsGetParser(&sOptions, &sOptionsForBinary, 1);
1147 0 : return argParser->usage();
1148 : }
1149 0 : catch (const std::exception &err)
1150 : {
1151 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
1152 0 : err.what());
1153 0 : return std::string();
1154 : }
1155 : }
1156 :
1157 : /************************************************************************/
1158 : /* CHECK_HAS_ENOUGH_ADDITIONAL_ARGS() */
1159 : /************************************************************************/
1160 :
1161 : #ifndef CheckHasEnoughAdditionalArgs_defined
1162 : #define CheckHasEnoughAdditionalArgs_defined
1163 :
1164 1 : static bool CheckHasEnoughAdditionalArgs(CSLConstList papszArgv, int i,
1165 : int nExtraArg, int nArgc)
1166 : {
1167 1 : if (i + nExtraArg >= nArgc)
1168 : {
1169 0 : CPLError(CE_Failure, CPLE_IllegalArg,
1170 0 : "%s option requires %d argument%s", papszArgv[i], nExtraArg,
1171 : nExtraArg == 1 ? "" : "s");
1172 0 : return false;
1173 : }
1174 1 : return true;
1175 : }
1176 : #endif
1177 :
1178 : #define CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(nExtraArg) \
1179 : if (!CheckHasEnoughAdditionalArgs(papszArgv, i, nExtraArg, nArgc)) \
1180 : { \
1181 : return nullptr; \
1182 : }
1183 :
1184 : /************************************************************************/
1185 : /* GDALGridOptionsNew() */
1186 : /************************************************************************/
1187 :
1188 : /**
1189 : * Allocates a GDALGridOptions struct.
1190 : *
1191 : * @param papszArgv NULL terminated list of options (potentially including
1192 : * filename and open options too), or NULL. The accepted options are the ones of
1193 : * the <a href="/programs/gdal_translate.html">gdal_translate</a> utility.
1194 : * @param psOptionsForBinary (output) may be NULL (and should generally be
1195 : * NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with
1196 : * GDALGridOptionsForBinaryNew() prior to this
1197 : * function. Will be filled with potentially present filename, open options,...
1198 : * @return pointer to the allocated GDALGridOptions struct. Must be freed with
1199 : * GDALGridOptionsFree().
1200 : *
1201 : * @since GDAL 2.1
1202 : */
1203 :
1204 : GDALGridOptions *
1205 106 : GDALGridOptionsNew(char **papszArgv,
1206 : GDALGridOptionsForBinary *psOptionsForBinary)
1207 : {
1208 212 : auto psOptions = std::make_unique<GDALGridOptions>();
1209 :
1210 : /* -------------------------------------------------------------------- */
1211 : /* Pre-processing for custom syntax that ArgumentParser does not */
1212 : /* support. */
1213 : /* -------------------------------------------------------------------- */
1214 :
1215 212 : CPLStringList aosArgv;
1216 106 : const int nArgc = CSLCount(papszArgv);
1217 106 : int nCountClipSrc = 0;
1218 1742 : for (int i = 0;
1219 1742 : i < nArgc && papszArgv != nullptr && papszArgv[i] != nullptr; i++)
1220 : {
1221 1636 : if (EQUAL(papszArgv[i], "-clipsrc"))
1222 : {
1223 1 : if (nCountClipSrc)
1224 : {
1225 0 : CPLError(CE_Failure, CPLE_AppDefined, "Duplicate argument %s",
1226 0 : papszArgv[i]);
1227 0 : return nullptr;
1228 : }
1229 : // argparse doesn't handle well variable number of values
1230 : // just before the positional arguments, so we have to detect
1231 : // it manually and set the correct number.
1232 1 : nCountClipSrc = 1;
1233 1 : CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1234 1 : if (CPLGetValueType(papszArgv[i + 1]) != CPL_VALUE_STRING &&
1235 0 : i + 4 < nArgc)
1236 : {
1237 0 : nCountClipSrc = 4;
1238 : }
1239 :
1240 3 : for (int j = 0; j < 1 + nCountClipSrc; ++j)
1241 : {
1242 2 : aosArgv.AddString(papszArgv[i]);
1243 2 : ++i;
1244 : }
1245 1 : --i;
1246 : }
1247 :
1248 : else
1249 : {
1250 1635 : aosArgv.AddString(papszArgv[i]);
1251 : }
1252 : }
1253 :
1254 : try
1255 : {
1256 : auto argParser = GDALGridOptionsGetParser(
1257 212 : psOptions.get(), psOptionsForBinary, nCountClipSrc);
1258 :
1259 106 : argParser->parse_args_without_binary_name(aosArgv.List());
1260 :
1261 212 : if (auto oTXE = argParser->present<std::vector<double>>("-txe"))
1262 : {
1263 106 : psOptions->dfXMin = (*oTXE)[0];
1264 106 : psOptions->dfXMax = (*oTXE)[1];
1265 106 : psOptions->bIsXExtentSet = true;
1266 : }
1267 :
1268 212 : if (auto oTYE = argParser->present<std::vector<double>>("-tye"))
1269 : {
1270 106 : psOptions->dfYMin = (*oTYE)[0];
1271 106 : psOptions->dfYMax = (*oTYE)[1];
1272 106 : psOptions->bIsYExtentSet = true;
1273 : }
1274 :
1275 211 : if (auto oOutsize = argParser->present<std::vector<int>>("-outsize"))
1276 : {
1277 105 : psOptions->nXSize = (*oOutsize)[0];
1278 105 : psOptions->nYSize = (*oOutsize)[1];
1279 : }
1280 :
1281 106 : if (auto adfTargetRes = argParser->present<std::vector<double>>("-tr"))
1282 : {
1283 1 : psOptions->dfXRes = (*adfTargetRes)[0];
1284 1 : psOptions->dfYRes = (*adfTargetRes)[1];
1285 1 : if (psOptions->dfXRes <= 0 || psOptions->dfYRes <= 0)
1286 : {
1287 0 : CPLError(CE_Failure, CPLE_IllegalArg,
1288 : "Wrong value for -tr parameters.");
1289 0 : return nullptr;
1290 : }
1291 : }
1292 :
1293 107 : if (auto oSpat = argParser->present<std::vector<double>>("-spat"))
1294 : {
1295 2 : OGRLinearRing oRing;
1296 1 : const double dfMinX = (*oSpat)[0];
1297 1 : const double dfMinY = (*oSpat)[1];
1298 1 : const double dfMaxX = (*oSpat)[2];
1299 1 : const double dfMaxY = (*oSpat)[3];
1300 :
1301 1 : oRing.addPoint(dfMinX, dfMinY);
1302 1 : oRing.addPoint(dfMinX, dfMaxY);
1303 1 : oRing.addPoint(dfMaxX, dfMaxY);
1304 1 : oRing.addPoint(dfMaxX, dfMinY);
1305 1 : oRing.addPoint(dfMinX, dfMinY);
1306 :
1307 2 : auto poPolygon = std::make_unique<OGRPolygon>();
1308 1 : poPolygon->addRing(&oRing);
1309 1 : psOptions->poSpatialFilter = std::move(poPolygon);
1310 : }
1311 :
1312 106 : if (auto oClipSrc =
1313 106 : argParser->present<std::vector<std::string>>("-clipsrc"))
1314 : {
1315 1 : const std::string &osVal = (*oClipSrc)[0];
1316 :
1317 1 : psOptions->poClipSrc.reset();
1318 1 : psOptions->osClipSrcDS.clear();
1319 :
1320 : VSIStatBufL sStat;
1321 1 : psOptions->bClipSrc = true;
1322 1 : if (oClipSrc->size() == 4)
1323 : {
1324 0 : const double dfMinX = CPLAtofM((*oClipSrc)[0].c_str());
1325 0 : const double dfMinY = CPLAtofM((*oClipSrc)[1].c_str());
1326 0 : const double dfMaxX = CPLAtofM((*oClipSrc)[2].c_str());
1327 0 : const double dfMaxY = CPLAtofM((*oClipSrc)[3].c_str());
1328 :
1329 0 : OGRLinearRing oRing;
1330 :
1331 0 : oRing.addPoint(dfMinX, dfMinY);
1332 0 : oRing.addPoint(dfMinX, dfMaxY);
1333 0 : oRing.addPoint(dfMaxX, dfMaxY);
1334 0 : oRing.addPoint(dfMaxX, dfMinY);
1335 0 : oRing.addPoint(dfMinX, dfMinY);
1336 :
1337 0 : auto poPoly = std::make_unique<OGRPolygon>();
1338 0 : poPoly->addRing(&oRing);
1339 0 : psOptions->poClipSrc = std::move(poPoly);
1340 : }
1341 1 : else if ((STARTS_WITH_CI(osVal.c_str(), "POLYGON") ||
1342 1 : STARTS_WITH_CI(osVal.c_str(), "MULTIPOLYGON")) &&
1343 0 : VSIStatL(osVal.c_str(), &sStat) != 0)
1344 : {
1345 0 : OGRGeometry *poGeom = nullptr;
1346 0 : OGRGeometryFactory::createFromWkt(osVal.c_str(), nullptr,
1347 : &poGeom);
1348 0 : psOptions->poClipSrc.reset(poGeom);
1349 0 : if (psOptions->poClipSrc == nullptr)
1350 : {
1351 0 : CPLError(CE_Failure, CPLE_IllegalArg,
1352 : "Invalid geometry. Must be a valid POLYGON or "
1353 : "MULTIPOLYGON WKT");
1354 0 : return nullptr;
1355 : }
1356 : }
1357 1 : else if (EQUAL(osVal.c_str(), "spat_extent"))
1358 : {
1359 : // Nothing to do
1360 : }
1361 : else
1362 : {
1363 1 : psOptions->osClipSrcDS = osVal;
1364 : }
1365 : }
1366 :
1367 106 : if (psOptions->bClipSrc && !psOptions->osClipSrcDS.empty())
1368 : {
1369 2 : psOptions->poClipSrc = LoadGeometry(
1370 1 : psOptions->osClipSrcDS, psOptions->osClipSrcSQL,
1371 2 : psOptions->osClipSrcLayer, psOptions->osClipSrcWhere);
1372 1 : if (!psOptions->poClipSrc)
1373 : {
1374 0 : CPLError(CE_Failure, CPLE_AppDefined,
1375 : "Cannot load source clip geometry.");
1376 0 : return nullptr;
1377 : }
1378 : }
1379 105 : else if (psOptions->bClipSrc && !psOptions->poClipSrc &&
1380 0 : !psOptions->poSpatialFilter)
1381 : {
1382 0 : CPLError(CE_Failure, CPLE_AppDefined,
1383 : "-clipsrc must be used with -spat option or \n"
1384 : "a bounding box, WKT string or datasource must be "
1385 : "specified.");
1386 0 : return nullptr;
1387 : }
1388 :
1389 106 : if (psOptions->poSpatialFilter)
1390 : {
1391 1 : if (psOptions->poClipSrc)
1392 : {
1393 : auto poTemp = std::unique_ptr<OGRGeometry>(
1394 0 : psOptions->poSpatialFilter->Intersection(
1395 0 : psOptions->poClipSrc.get()));
1396 0 : if (poTemp)
1397 : {
1398 0 : psOptions->poSpatialFilter = std::move(poTemp);
1399 : }
1400 :
1401 0 : psOptions->poClipSrc.reset();
1402 : }
1403 : }
1404 : else
1405 : {
1406 105 : if (psOptions->poClipSrc)
1407 : {
1408 1 : psOptions->poSpatialFilter = std::move(psOptions->poClipSrc);
1409 : }
1410 : }
1411 :
1412 106 : return psOptions.release();
1413 : }
1414 0 : catch (const std::exception &err)
1415 : {
1416 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s", err.what());
1417 0 : return nullptr;
1418 : }
1419 : }
1420 :
1421 : /************************************************************************/
1422 : /* GDALGridOptionsFree() */
1423 : /************************************************************************/
1424 :
1425 : /**
1426 : * Frees the GDALGridOptions struct.
1427 : *
1428 : * @param psOptions the options struct for GDALGrid().
1429 : *
1430 : * @since GDAL 2.1
1431 : */
1432 :
1433 106 : void GDALGridOptionsFree(GDALGridOptions *psOptions)
1434 : {
1435 106 : delete psOptions;
1436 106 : }
1437 :
1438 : /************************************************************************/
1439 : /* GDALGridOptionsSetProgress() */
1440 : /************************************************************************/
1441 :
1442 : /**
1443 : * Set a progress function.
1444 : *
1445 : * @param psOptions the options struct for GDALGrid().
1446 : * @param pfnProgress the progress callback.
1447 : * @param pProgressData the user data for the progress callback.
1448 : *
1449 : * @since GDAL 2.1
1450 : */
1451 :
1452 54 : void GDALGridOptionsSetProgress(GDALGridOptions *psOptions,
1453 : GDALProgressFunc pfnProgress,
1454 : void *pProgressData)
1455 : {
1456 54 : psOptions->pfnProgress = pfnProgress;
1457 54 : psOptions->pProgressData = pProgressData;
1458 54 : if (pfnProgress == GDALTermProgress)
1459 54 : psOptions->bQuiet = false;
1460 54 : }
1461 :
1462 : #undef CHECK_HAS_ENOUGH_ADDITIONAL_ARGS
|
