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