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