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