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