Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: gdal "vector grid" subcommand
5 : * Author: Even Rouault <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalalg_vector_grid.h"
14 : #include "gdalalg_vector_grid_average.h"
15 : #include "gdalalg_vector_grid_data_metrics.h"
16 : #include "gdalalg_vector_grid_invdist.h"
17 : #include "gdalalg_vector_grid_invdistnn.h"
18 : #include "gdalalg_vector_grid_linear.h"
19 : #include "gdalalg_vector_grid_nearest.h"
20 :
21 : #include "cpl_conv.h"
22 : #include "gdal_priv.h"
23 : #include "gdal_utils.h"
24 : #include "ogrsf_frmts.h"
25 :
26 : //! @cond Doxygen_Suppress
27 :
28 : #ifndef _
29 : #define _(x) (x)
30 : #endif
31 :
32 : /************************************************************************/
33 : /* GDALVectorGridAlgorithm::GDALVectorGridAlgorithm() */
34 : /************************************************************************/
35 :
36 111 : GDALVectorGridAlgorithm::GDALVectorGridAlgorithm(bool standaloneStep)
37 : : GDALPipelineStepAlgorithm(
38 : NAME, DESCRIPTION, HELP_URL,
39 0 : ConstructorOptions()
40 111 : .SetStandaloneStep(standaloneStep)
41 222 : .SetOutputFormatCreateCapability(GDAL_DCAP_CREATE))
42 : {
43 111 : if (standaloneStep)
44 : {
45 91 : RegisterSubAlgorithm<GDALVectorGridAverageAlgorithmStandalone>();
46 91 : RegisterSubAlgorithm<GDALVectorGridInvdistAlgorithmStandalone>();
47 91 : RegisterSubAlgorithm<GDALVectorGridInvdistNNAlgorithmStandalone>();
48 91 : RegisterSubAlgorithm<GDALVectorGridLinearAlgorithmStandalone>();
49 91 : RegisterSubAlgorithm<GDALVectorGridNearestAlgorithmStandalone>();
50 91 : RegisterSubAlgorithm<GDALVectorGridMinimumAlgorithmStandalone>();
51 91 : RegisterSubAlgorithm<GDALVectorGridMaximumAlgorithmStandalone>();
52 91 : RegisterSubAlgorithm<GDALVectorGridRangeAlgorithmStandalone>();
53 91 : RegisterSubAlgorithm<GDALVectorGridCountAlgorithmStandalone>();
54 : RegisterSubAlgorithm<
55 91 : GDALVectorGridAverageDistanceAlgorithmStandalone>();
56 : RegisterSubAlgorithm<
57 91 : GDALVectorGridAverageDistancePointsAlgorithmStandalone>();
58 : }
59 : else
60 : {
61 20 : RegisterSubAlgorithm<GDALVectorGridAverageAlgorithm>();
62 20 : RegisterSubAlgorithm<GDALVectorGridInvdistAlgorithm>();
63 20 : RegisterSubAlgorithm<GDALVectorGridInvdistNNAlgorithm>();
64 20 : RegisterSubAlgorithm<GDALVectorGridLinearAlgorithm>();
65 20 : RegisterSubAlgorithm<GDALVectorGridNearestAlgorithm>();
66 20 : RegisterSubAlgorithm<GDALVectorGridMinimumAlgorithm>();
67 20 : RegisterSubAlgorithm<GDALVectorGridMaximumAlgorithm>();
68 20 : RegisterSubAlgorithm<GDALVectorGridRangeAlgorithm>();
69 20 : RegisterSubAlgorithm<GDALVectorGridCountAlgorithm>();
70 20 : RegisterSubAlgorithm<GDALVectorGridAverageDistanceAlgorithm>();
71 20 : RegisterSubAlgorithm<GDALVectorGridAverageDistancePointsAlgorithm>();
72 : }
73 111 : }
74 :
75 : /************************************************************************/
76 : /* GDALVectorGridAlgorithm::RunStep() */
77 : /************************************************************************/
78 :
79 0 : bool GDALVectorGridAlgorithm::RunStep(GDALPipelineStepRunContext &)
80 : {
81 0 : CPLError(CE_Failure, CPLE_AppDefined,
82 : "The Run() method should not be called directly on the \"gdal "
83 : "vector grid\" program.");
84 0 : return false;
85 : }
86 :
87 : /************************************************************************/
88 : /* GDALVectorGeomAlgorithm::RunImpl() */
89 : /************************************************************************/
90 :
91 1 : bool GDALVectorGridAlgorithm::RunImpl(GDALProgressFunc, void *)
92 : {
93 1 : CPLError(CE_Failure, CPLE_AppDefined,
94 : "The Run() method should not be called directly on the \"gdal "
95 : "vector grid\" program.");
96 1 : return false;
97 : }
98 :
99 : /************************************************************************/
100 : /* GDALVectorGridAbstractAlgorithm() */
101 : /************************************************************************/
102 :
103 134 : GDALVectorGridAbstractAlgorithm::GDALVectorGridAbstractAlgorithm(
104 : const std::string &name, const std::string &description,
105 134 : const std::string &helpURL, bool standaloneStep)
106 : : GDALPipelineStepAlgorithm(
107 : name, description, helpURL,
108 0 : ConstructorOptions()
109 134 : .SetStandaloneStep(standaloneStep)
110 268 : .SetOutputFormatCreateCapability(GDAL_DCAP_CREATE))
111 : {
112 134 : AddProgressArg();
113 134 : if (standaloneStep)
114 : {
115 111 : AddOutputFormatArg(&m_format).AddMetadataItem(
116 333 : GAAMDI_REQUIRED_CAPABILITIES, {GDAL_DCAP_RASTER, GDAL_DCAP_CREATE});
117 111 : AddOpenOptionsArg(&m_openOptions);
118 111 : AddInputFormatsArg(&m_inputFormats)
119 222 : .AddMetadataItem(GAAMDI_REQUIRED_CAPABILITIES, {GDAL_DCAP_VECTOR});
120 111 : AddInputDatasetArg(&m_inputDataset, GDAL_OF_VECTOR);
121 111 : AddOutputDatasetArg(&m_outputDataset, GDAL_OF_RASTER);
122 111 : AddCreationOptionsArg(&m_creationOptions);
123 : }
124 : AddArg("extent", 0, _("Set the target georeferenced extent"),
125 268 : &m_targetExtent)
126 134 : .SetMinCount(4)
127 134 : .SetMaxCount(4)
128 134 : .SetRepeatedArgAllowed(false)
129 134 : .SetMetaVar("<xmin>,<ymin>,<xmax>,<ymax>");
130 268 : AddArg("resolution", 0, _("Set the target resolution"), &m_targetResolution)
131 134 : .SetMinCount(2)
132 134 : .SetMaxCount(2)
133 134 : .SetRepeatedArgAllowed(false)
134 268 : .SetMetaVar("<xres>,<yres>")
135 134 : .SetMutualExclusionGroup("size-or-resolution");
136 : AddArg("size", 0, _("Set the target size in pixels and lines"),
137 268 : &m_targetSize)
138 134 : .SetMinCount(2)
139 134 : .SetMaxCount(2)
140 134 : .SetRepeatedArgAllowed(false)
141 268 : .SetMetaVar("<xsize>,<ysize>")
142 134 : .SetMutualExclusionGroup("size-or-resolution");
143 134 : AddOutputDataTypeArg(&m_outputType).SetDefault(m_outputType);
144 268 : AddArg("crs", 0, _("Override the projection for the output file"), &m_crs)
145 268 : .AddHiddenAlias("srs")
146 134 : .SetIsCRSArg(/*noneAllowed=*/false);
147 134 : AddOverwriteArg(&m_overwrite);
148 134 : AddLayerNameArg(&m_layers)
149 268 : .SetMutualExclusionGroup("layer-sql")
150 134 : .AddAlias("layer");
151 268 : AddArg("sql", 0, _("SQL statement"), &m_sql)
152 134 : .SetReadFromFileAtSyntaxAllowed()
153 268 : .SetMetaVar("<statement>|@<filename>")
154 134 : .SetRemoveSQLCommentsEnabled()
155 134 : .SetMutualExclusionGroup("layer-sql");
156 : AddBBOXArg(
157 : &m_bbox,
158 134 : _("Select only points contained within the specified bounding box"));
159 268 : AddArg("zfield", 0, _("Field name from which to get Z values."), &m_zField)
160 134 : .AddHiddenAlias("z-field");
161 : AddArg("zoffset", 0,
162 : _("Value to add to the Z field value (applied before zmultiply)"),
163 268 : &m_zOffset)
164 134 : .SetDefault(m_zOffset)
165 134 : .AddHiddenAlias("z-offset");
166 : AddArg("zmultiply", 0,
167 : _("Multiplication factor for the Z field value (applied after "
168 : "zoffset)"),
169 268 : &m_zMultiply)
170 134 : .SetDefault(m_zMultiply)
171 134 : .AddHiddenAlias("z-multiply");
172 :
173 134 : AddValidationAction(
174 93 : [this]()
175 : {
176 90 : bool ret = true;
177 90 : if (!m_targetResolution.empty() && m_targetExtent.empty())
178 : {
179 1 : ReportError(CE_Failure, CPLE_AppDefined,
180 : "'resolution' should be defined when 'extent' is.");
181 1 : ret = false;
182 : }
183 90 : return ret;
184 : });
185 134 : }
186 :
187 : /************************************************************************/
188 : /* GDALVectorGridAbstractAlgorithm::RunStep() */
189 : /************************************************************************/
190 :
191 77 : bool GDALVectorGridAbstractAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
192 : {
193 77 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
194 77 : CPLAssert(poSrcDS);
195 77 : CPLAssert(!m_outputDataset.GetDatasetRef());
196 :
197 154 : CPLStringList aosOptions;
198 :
199 154 : std::string outputFilename;
200 77 : if (m_standaloneStep)
201 : {
202 76 : outputFilename = m_outputDataset.GetName();
203 76 : if (!m_format.empty())
204 : {
205 72 : aosOptions.AddString("-of");
206 72 : aosOptions.AddString(m_format.c_str());
207 : }
208 :
209 77 : for (const std::string &co : m_creationOptions)
210 : {
211 1 : aosOptions.AddString("-co");
212 1 : aosOptions.AddString(co.c_str());
213 : }
214 : }
215 : else
216 : {
217 1 : outputFilename = CPLGenerateTempFilenameSafe("_grid.tif");
218 :
219 1 : aosOptions.AddString("-of");
220 1 : aosOptions.AddString("GTiff");
221 :
222 1 : aosOptions.AddString("-co");
223 1 : aosOptions.AddString("TILED=YES");
224 : }
225 :
226 77 : if (!m_targetExtent.empty())
227 : {
228 6 : aosOptions.AddString("-txe");
229 6 : aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[0]));
230 6 : aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[2]));
231 6 : aosOptions.AddString("-tye");
232 6 : aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[1]));
233 6 : aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[3]));
234 : }
235 :
236 77 : if (!m_bbox.empty())
237 : {
238 2 : aosOptions.AddString("-clipsrc");
239 10 : for (double v : m_bbox)
240 : {
241 8 : aosOptions.AddString(CPLSPrintf("%.17g", v));
242 : }
243 : }
244 :
245 77 : if (!m_targetResolution.empty())
246 : {
247 1 : aosOptions.AddString("-tr");
248 3 : for (double targetResolution : m_targetResolution)
249 : {
250 2 : aosOptions.AddString(CPLSPrintf("%.17g", targetResolution));
251 : }
252 : }
253 :
254 77 : if (!m_targetSize.empty())
255 : {
256 1 : aosOptions.AddString("-outsize");
257 3 : for (int targetSize : m_targetSize)
258 : {
259 2 : aosOptions.AddString(CPLSPrintf("%d", targetSize));
260 : }
261 : }
262 :
263 77 : if (!m_outputType.empty())
264 : {
265 77 : aosOptions.AddString("-ot");
266 77 : aosOptions.AddString(m_outputType.c_str());
267 : }
268 :
269 77 : if (!m_crs.empty())
270 : {
271 1 : aosOptions.AddString("-a_srs");
272 1 : aosOptions.AddString(m_crs.c_str());
273 : }
274 :
275 77 : if (m_sql.empty())
276 : {
277 78 : for (const std::string &layer : m_layers)
278 : {
279 3 : aosOptions.AddString("-l");
280 3 : aosOptions.AddString(layer.c_str());
281 : }
282 : }
283 : else
284 : {
285 2 : aosOptions.AddString("-sql");
286 2 : aosOptions.AddString(m_sql.c_str());
287 : }
288 :
289 77 : if (m_zOffset != 0)
290 : {
291 1 : aosOptions.AddString("-z_increase");
292 1 : aosOptions.AddString(CPLSPrintf("%.17g", m_zOffset));
293 : }
294 :
295 77 : if (m_zMultiply != 0)
296 : {
297 77 : aosOptions.AddString("-z_multiply");
298 77 : aosOptions.AddString(CPLSPrintf("%.17g", m_zMultiply));
299 : }
300 :
301 77 : if (!m_zField.empty())
302 : {
303 1 : aosOptions.AddString("-zfield");
304 1 : aosOptions.AddString(m_zField.c_str());
305 : }
306 76 : else if (m_sql.empty())
307 : {
308 75 : const auto CheckLayer = [this](OGRLayer *poLayer)
309 : {
310 : auto poFeat =
311 146 : std::unique_ptr<OGRFeature>(poLayer->GetNextFeature());
312 73 : poLayer->ResetReading();
313 73 : if (poFeat)
314 : {
315 73 : const auto poGeom = poFeat->GetGeometryRef();
316 73 : if (!poGeom || !poGeom->Is3D())
317 : {
318 2 : ReportError(
319 : CE_Warning, CPLE_AppDefined,
320 : "At least one geometry of layer '%s' lacks a Z "
321 : "component. You may need to set the 'zfield' argument",
322 2 : poLayer->GetName());
323 2 : return false;
324 : }
325 : }
326 71 : return true;
327 74 : };
328 :
329 74 : if (m_layers.empty())
330 : {
331 142 : for (auto poLayer : poSrcDS->GetLayers())
332 : {
333 71 : if (!CheckLayer(poLayer))
334 1 : break;
335 : }
336 : }
337 : else
338 : {
339 5 : for (const std::string &layerName : m_layers)
340 : {
341 3 : auto poLayer = poSrcDS->GetLayerByName(layerName.c_str());
342 3 : if (poLayer && !CheckLayer(poLayer))
343 1 : break;
344 : }
345 : }
346 : }
347 :
348 77 : aosOptions.AddString("-a");
349 77 : aosOptions.AddString(GetGridAlgorithm().c_str());
350 :
351 77 : bool bOK = false;
352 : std::unique_ptr<GDALGridOptions, decltype(&GDALGridOptionsFree)> psOptions{
353 77 : GDALGridOptionsNew(aosOptions.List(), nullptr), GDALGridOptionsFree};
354 :
355 77 : if (psOptions)
356 : {
357 77 : GDALGridOptionsSetProgress(psOptions.get(), ctxt.m_pfnProgress,
358 : ctxt.m_pProgressData);
359 :
360 : auto poRetDS = std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(
361 : GDALGrid(outputFilename.c_str(), GDALDataset::ToHandle(poSrcDS),
362 154 : psOptions.get(), nullptr)));
363 :
364 77 : if (poRetDS)
365 : {
366 75 : bOK = true;
367 :
368 75 : if (!m_standaloneStep)
369 : {
370 1 : VSIUnlink(outputFilename.c_str());
371 1 : poRetDS->MarkSuppressOnClose();
372 : }
373 :
374 75 : m_outputDataset.Set(std::move(poRetDS));
375 : }
376 : }
377 :
378 154 : return bOK;
379 : }
380 :
381 : /************************************************************************/
382 : /* GDALVectorGridAbstractAlgorithm::RunImpl() */
383 : /************************************************************************/
384 :
385 76 : bool GDALVectorGridAbstractAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
386 : void *pProgressData)
387 : {
388 76 : GDALPipelineStepRunContext stepCtxt;
389 76 : stepCtxt.m_pfnProgress = pfnProgress;
390 76 : stepCtxt.m_pProgressData = pProgressData;
391 76 : return RunPreStepPipelineValidations() && RunStep(stepCtxt);
392 : }
393 :
394 : /************************************************************************/
395 : /* GDALVectorGridAbstractAlgorithm::AddRadiusArg() */
396 : /************************************************************************/
397 :
398 134 : GDALInConstructionAlgorithmArg &GDALVectorGridAbstractAlgorithm::AddRadiusArg()
399 : {
400 268 : return AddArg("radius", 0, _("Radius of the search circle"), &m_radius)
401 268 : .SetMutualExclusionGroup("radius");
402 : }
403 :
404 : /************************************************************************/
405 : /* GDALVectorGridAbstractAlgorithm::AddRadius1AndRadius2Arg() */
406 : /************************************************************************/
407 :
408 112 : void GDALVectorGridAbstractAlgorithm::AddRadius1AndRadius2Arg()
409 : {
410 224 : AddArg("radius1", 0, _("First axis of the search ellipse"), &m_radius1)
411 112 : .SetMutualExclusionGroup("radius");
412 112 : AddArg("radius2", 0, _("Second axis of the search ellipse"), &m_radius2);
413 :
414 112 : AddValidationAction(
415 180 : [this]()
416 : {
417 76 : bool ret = true;
418 76 : if (m_radius1 > 0 && m_radius2 == 0)
419 : {
420 2 : ReportError(CE_Failure, CPLE_AppDefined,
421 : "'radius2' should be defined when 'radius1' is.");
422 2 : ret = false;
423 : }
424 74 : else if (m_radius2 > 0 && m_radius1 == 0)
425 : {
426 2 : ReportError(CE_Failure, CPLE_AppDefined,
427 : "'radius1' should be defined when 'radius2' is.");
428 2 : ret = false;
429 : }
430 76 : return ret;
431 : });
432 112 : }
433 :
434 : /************************************************************************/
435 : /* GDALVectorGridAbstractAlgorithm::AddAngleArg() */
436 : /************************************************************************/
437 :
438 112 : GDALInConstructionAlgorithmArg &GDALVectorGridAbstractAlgorithm::AddAngleArg()
439 : {
440 : return AddArg("angle", 0,
441 : _("Angle of search ellipse rotation in degrees (counter "
442 : "clockwise)"),
443 224 : &m_angle)
444 224 : .SetDefault(m_angle);
445 : }
446 :
447 : /************************************************************************/
448 : /* GDALVectorGridAbstractAlgorithm::AddMinPointsArg() */
449 : /************************************************************************/
450 :
451 : GDALInConstructionAlgorithmArg &
452 116 : GDALVectorGridAbstractAlgorithm::AddMinPointsArg()
453 : {
454 : return AddArg("min-points", 0, _("Minimum number of data points to use"),
455 232 : &m_minPoints)
456 232 : .SetDefault(m_minPoints);
457 : }
458 :
459 : /************************************************************************/
460 : /* GDALVectorGridAbstractAlgorithm::AddMaxPointsArg() */
461 : /************************************************************************/
462 :
463 : GDALInConstructionAlgorithmArg &
464 80 : GDALVectorGridAbstractAlgorithm::AddMaxPointsArg()
465 : {
466 : return AddArg("max-points", 0, _("Maximum number of data points to use"),
467 160 : &m_maxPoints)
468 160 : .SetDefault(m_maxPoints);
469 : }
470 :
471 : /************************************************************************/
472 : /* GDALVectorGridAbstractAlgorithm::AddMinMaxPointsPerQuadrantArg() */
473 : /************************************************************************/
474 :
475 116 : void GDALVectorGridAbstractAlgorithm::AddMinMaxPointsPerQuadrantArg()
476 : {
477 : AddArg("min-points-per-quadrant", 0,
478 : _("Minimum number of data points to use per quadrant"),
479 232 : &m_minPointsPerQuadrant)
480 116 : .SetDefault(m_minPointsPerQuadrant);
481 : AddArg("max-points-per-quadrant", 0,
482 : _("Maximum number of data points to use per quadrant"),
483 232 : &m_maxPointsPerQuadrant)
484 116 : .SetDefault(m_maxPointsPerQuadrant);
485 116 : }
486 :
487 : /************************************************************************/
488 : /* GDALVectorGridAbstractAlgorithm::AddNodataArg() */
489 : /************************************************************************/
490 :
491 134 : GDALInConstructionAlgorithmArg &GDALVectorGridAbstractAlgorithm::AddNodataArg()
492 : {
493 268 : return AddArg("nodata", 0, _("Target nodata value"), &m_nodata)
494 268 : .SetDefault(m_nodata);
495 : }
496 :
497 : GDALVectorGridAlgorithmStandalone::~GDALVectorGridAlgorithmStandalone() =
498 : default;
499 :
500 : //! @endcond
|
