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