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 97 : GDALVectorGridAlgorithm::GDALVectorGridAlgorithm(bool standaloneStep)
37 : : GDALPipelineStepAlgorithm(
38 : NAME, DESCRIPTION, HELP_URL,
39 0 : ConstructorOptions()
40 97 : .SetStandaloneStep(standaloneStep)
41 194 : .SetOutputFormatCreateCapability(GDAL_DCAP_CREATE))
42 : {
43 97 : 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 6 : RegisterSubAlgorithm<GDALVectorGridAverageAlgorithm>();
62 6 : RegisterSubAlgorithm<GDALVectorGridInvdistAlgorithm>();
63 6 : RegisterSubAlgorithm<GDALVectorGridInvdistNNAlgorithm>();
64 6 : RegisterSubAlgorithm<GDALVectorGridLinearAlgorithm>();
65 6 : RegisterSubAlgorithm<GDALVectorGridNearestAlgorithm>();
66 6 : RegisterSubAlgorithm<GDALVectorGridMinimumAlgorithm>();
67 6 : RegisterSubAlgorithm<GDALVectorGridMaximumAlgorithm>();
68 6 : RegisterSubAlgorithm<GDALVectorGridRangeAlgorithm>();
69 6 : RegisterSubAlgorithm<GDALVectorGridCountAlgorithm>();
70 6 : RegisterSubAlgorithm<GDALVectorGridAverageDistanceAlgorithm>();
71 6 : RegisterSubAlgorithm<GDALVectorGridAverageDistancePointsAlgorithm>();
72 : }
73 97 : }
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).SetMutualExclusionGroup("layer-sql");
149 268 : AddArg("sql", 0, _("SQL statement"), &m_sql)
150 134 : .SetReadFromFileAtSyntaxAllowed()
151 268 : .SetMetaVar("<statement>|@<filename>")
152 134 : .SetRemoveSQLCommentsEnabled()
153 134 : .SetMutualExclusionGroup("layer-sql");
154 : AddBBOXArg(
155 : &m_bbox,
156 134 : _("Select only points contained within the specified bounding box"));
157 268 : AddArg("zfield", 0, _("Field name from which to get Z values."), &m_zField)
158 134 : .AddHiddenAlias("z-field");
159 : AddArg("zoffset", 0,
160 : _("Value to add to the Z field value (applied before zmultiply)"),
161 268 : &m_zOffset)
162 134 : .SetDefault(m_zOffset)
163 134 : .AddHiddenAlias("z-offset");
164 : AddArg("zmultiply", 0,
165 : _("Multiplication factor for the Z field value (applied after "
166 : "zoffset)"),
167 268 : &m_zMultiply)
168 134 : .SetDefault(m_zMultiply)
169 134 : .AddHiddenAlias("z-multiply");
170 :
171 134 : AddValidationAction(
172 93 : [this]()
173 : {
174 90 : bool ret = true;
175 90 : if (!m_targetResolution.empty() && m_targetExtent.empty())
176 : {
177 1 : ReportError(CE_Failure, CPLE_AppDefined,
178 : "'resolution' should be defined when 'extent' is.");
179 1 : ret = false;
180 : }
181 90 : return ret;
182 : });
183 134 : }
184 :
185 : /************************************************************************/
186 : /* GDALVectorGridAbstractAlgorithm::RunStep() */
187 : /************************************************************************/
188 :
189 77 : bool GDALVectorGridAbstractAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
190 : {
191 77 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
192 77 : CPLAssert(poSrcDS);
193 77 : CPLAssert(!m_outputDataset.GetDatasetRef());
194 :
195 154 : CPLStringList aosOptions;
196 :
197 154 : std::string outputFilename;
198 77 : if (m_standaloneStep)
199 : {
200 76 : outputFilename = m_outputDataset.GetName();
201 76 : if (!m_format.empty())
202 : {
203 72 : aosOptions.AddString("-of");
204 72 : aosOptions.AddString(m_format.c_str());
205 : }
206 :
207 77 : for (const std::string &co : m_creationOptions)
208 : {
209 1 : aosOptions.AddString("-co");
210 1 : aosOptions.AddString(co.c_str());
211 : }
212 : }
213 : else
214 : {
215 1 : outputFilename = CPLGenerateTempFilenameSafe("_grid.tif");
216 :
217 1 : aosOptions.AddString("-of");
218 1 : aosOptions.AddString("GTiff");
219 :
220 1 : aosOptions.AddString("-co");
221 1 : aosOptions.AddString("TILED=YES");
222 : }
223 :
224 77 : if (!m_targetExtent.empty())
225 : {
226 6 : aosOptions.AddString("-txe");
227 6 : aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[0]));
228 6 : aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[2]));
229 6 : aosOptions.AddString("-tye");
230 6 : aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[1]));
231 6 : aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[3]));
232 : }
233 :
234 77 : if (!m_bbox.empty())
235 : {
236 2 : aosOptions.AddString("-clipsrc");
237 10 : for (double v : m_bbox)
238 : {
239 8 : aosOptions.AddString(CPLSPrintf("%.17g", v));
240 : }
241 : }
242 :
243 77 : if (!m_targetResolution.empty())
244 : {
245 1 : aosOptions.AddString("-tr");
246 3 : for (double targetResolution : m_targetResolution)
247 : {
248 2 : aosOptions.AddString(CPLSPrintf("%.17g", targetResolution));
249 : }
250 : }
251 :
252 77 : if (!m_targetSize.empty())
253 : {
254 1 : aosOptions.AddString("-outsize");
255 3 : for (int targetSize : m_targetSize)
256 : {
257 2 : aosOptions.AddString(CPLSPrintf("%d", targetSize));
258 : }
259 : }
260 :
261 77 : if (!m_outputType.empty())
262 : {
263 77 : aosOptions.AddString("-ot");
264 77 : aosOptions.AddString(m_outputType.c_str());
265 : }
266 :
267 77 : if (!m_crs.empty())
268 : {
269 1 : aosOptions.AddString("-a_srs");
270 1 : aosOptions.AddString(m_crs.c_str());
271 : }
272 :
273 77 : if (m_sql.empty())
274 : {
275 78 : for (const std::string &layer : m_layers)
276 : {
277 3 : aosOptions.AddString("-l");
278 3 : aosOptions.AddString(layer.c_str());
279 : }
280 : }
281 : else
282 : {
283 2 : aosOptions.AddString("-sql");
284 2 : aosOptions.AddString(m_sql.c_str());
285 : }
286 :
287 77 : if (m_zOffset != 0)
288 : {
289 1 : aosOptions.AddString("-z_increase");
290 1 : aosOptions.AddString(CPLSPrintf("%.17g", m_zOffset));
291 : }
292 :
293 77 : if (m_zMultiply != 0)
294 : {
295 77 : aosOptions.AddString("-z_multiply");
296 77 : aosOptions.AddString(CPLSPrintf("%.17g", m_zMultiply));
297 : }
298 :
299 77 : if (!m_zField.empty())
300 : {
301 1 : aosOptions.AddString("-zfield");
302 1 : aosOptions.AddString(m_zField.c_str());
303 : }
304 76 : else if (m_sql.empty())
305 : {
306 75 : const auto CheckLayer = [this](OGRLayer *poLayer)
307 : {
308 : auto poFeat =
309 146 : std::unique_ptr<OGRFeature>(poLayer->GetNextFeature());
310 73 : poLayer->ResetReading();
311 73 : if (poFeat)
312 : {
313 73 : const auto poGeom = poFeat->GetGeometryRef();
314 73 : if (!poGeom || !poGeom->Is3D())
315 : {
316 2 : ReportError(
317 : CE_Warning, CPLE_AppDefined,
318 : "At least one geometry of layer '%s' lacks a Z "
319 : "component. You may need to set the 'zfield' argument",
320 2 : poLayer->GetName());
321 2 : return false;
322 : }
323 : }
324 71 : return true;
325 74 : };
326 :
327 74 : if (m_layers.empty())
328 : {
329 142 : for (auto poLayer : poSrcDS->GetLayers())
330 : {
331 71 : if (!CheckLayer(poLayer))
332 1 : break;
333 : }
334 : }
335 : else
336 : {
337 5 : for (const std::string &layerName : m_layers)
338 : {
339 3 : auto poLayer = poSrcDS->GetLayerByName(layerName.c_str());
340 3 : if (poLayer && !CheckLayer(poLayer))
341 1 : break;
342 : }
343 : }
344 : }
345 :
346 77 : aosOptions.AddString("-a");
347 77 : aosOptions.AddString(GetGridAlgorithm().c_str());
348 :
349 77 : bool bOK = false;
350 : std::unique_ptr<GDALGridOptions, decltype(&GDALGridOptionsFree)> psOptions{
351 77 : GDALGridOptionsNew(aosOptions.List(), nullptr), GDALGridOptionsFree};
352 :
353 77 : if (psOptions)
354 : {
355 77 : GDALGridOptionsSetProgress(psOptions.get(), ctxt.m_pfnProgress,
356 : ctxt.m_pProgressData);
357 :
358 : auto poRetDS = std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(
359 : GDALGrid(outputFilename.c_str(), GDALDataset::ToHandle(poSrcDS),
360 154 : psOptions.get(), nullptr)));
361 :
362 77 : if (poRetDS)
363 : {
364 75 : bOK = true;
365 :
366 75 : if (!m_standaloneStep)
367 : {
368 1 : VSIUnlink(outputFilename.c_str());
369 1 : poRetDS->MarkSuppressOnClose();
370 : }
371 :
372 75 : m_outputDataset.Set(std::move(poRetDS));
373 : }
374 : }
375 :
376 154 : return bOK;
377 : }
378 :
379 : /************************************************************************/
380 : /* GDALVectorGridAbstractAlgorithm::RunImpl() */
381 : /************************************************************************/
382 :
383 76 : bool GDALVectorGridAbstractAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
384 : void *pProgressData)
385 : {
386 76 : GDALPipelineStepRunContext stepCtxt;
387 76 : stepCtxt.m_pfnProgress = pfnProgress;
388 76 : stepCtxt.m_pProgressData = pProgressData;
389 76 : return RunPreStepPipelineValidations() && RunStep(stepCtxt);
390 : }
391 :
392 : /************************************************************************/
393 : /* GDALVectorGridAbstractAlgorithm::AddRadiusArg() */
394 : /************************************************************************/
395 :
396 134 : GDALInConstructionAlgorithmArg &GDALVectorGridAbstractAlgorithm::AddRadiusArg()
397 : {
398 268 : return AddArg("radius", 0, _("Radius of the search circle"), &m_radius)
399 268 : .SetMutualExclusionGroup("radius");
400 : }
401 :
402 : /************************************************************************/
403 : /* GDALVectorGridAbstractAlgorithm::AddRadius1AndRadius2Arg() */
404 : /************************************************************************/
405 :
406 112 : void GDALVectorGridAbstractAlgorithm::AddRadius1AndRadius2Arg()
407 : {
408 224 : AddArg("radius1", 0, _("First axis of the search ellipse"), &m_radius1)
409 112 : .SetMutualExclusionGroup("radius");
410 112 : AddArg("radius2", 0, _("Second axis of the search ellipse"), &m_radius2);
411 :
412 112 : AddValidationAction(
413 180 : [this]()
414 : {
415 76 : bool ret = true;
416 76 : if (m_radius1 > 0 && m_radius2 == 0)
417 : {
418 2 : ReportError(CE_Failure, CPLE_AppDefined,
419 : "'radius2' should be defined when 'radius1' is.");
420 2 : ret = false;
421 : }
422 74 : else if (m_radius2 > 0 && m_radius1 == 0)
423 : {
424 2 : ReportError(CE_Failure, CPLE_AppDefined,
425 : "'radius1' should be defined when 'radius2' is.");
426 2 : ret = false;
427 : }
428 76 : return ret;
429 : });
430 112 : }
431 :
432 : /************************************************************************/
433 : /* GDALVectorGridAbstractAlgorithm::AddAngleArg() */
434 : /************************************************************************/
435 :
436 112 : GDALInConstructionAlgorithmArg &GDALVectorGridAbstractAlgorithm::AddAngleArg()
437 : {
438 : return AddArg("angle", 0,
439 : _("Angle of search ellipse rotation in degrees (counter "
440 : "clockwise)"),
441 224 : &m_angle)
442 224 : .SetDefault(m_angle);
443 : }
444 :
445 : /************************************************************************/
446 : /* GDALVectorGridAbstractAlgorithm::AddMinPointsArg() */
447 : /************************************************************************/
448 :
449 : GDALInConstructionAlgorithmArg &
450 116 : GDALVectorGridAbstractAlgorithm::AddMinPointsArg()
451 : {
452 : return AddArg("min-points", 0, _("Minimum number of data points to use"),
453 232 : &m_minPoints)
454 232 : .SetDefault(m_minPoints);
455 : }
456 :
457 : /************************************************************************/
458 : /* GDALVectorGridAbstractAlgorithm::AddMaxPointsArg() */
459 : /************************************************************************/
460 :
461 : GDALInConstructionAlgorithmArg &
462 80 : GDALVectorGridAbstractAlgorithm::AddMaxPointsArg()
463 : {
464 : return AddArg("max-points", 0, _("Maximum number of data points to use"),
465 160 : &m_maxPoints)
466 160 : .SetDefault(m_maxPoints);
467 : }
468 :
469 : /************************************************************************/
470 : /* GDALVectorGridAbstractAlgorithm::AddMinMaxPointsPerQuadrantArg() */
471 : /************************************************************************/
472 :
473 116 : void GDALVectorGridAbstractAlgorithm::AddMinMaxPointsPerQuadrantArg()
474 : {
475 : AddArg("min-points-per-quadrant", 0,
476 : _("Minimum number of data points to use per quadrant"),
477 232 : &m_minPointsPerQuadrant)
478 116 : .SetDefault(m_minPointsPerQuadrant);
479 : AddArg("max-points-per-quadrant", 0,
480 : _("Maximum number of data points to use per quadrant"),
481 232 : &m_maxPointsPerQuadrant)
482 116 : .SetDefault(m_maxPointsPerQuadrant);
483 116 : }
484 :
485 : /************************************************************************/
486 : /* GDALVectorGridAbstractAlgorithm::AddNodataArg() */
487 : /************************************************************************/
488 :
489 134 : GDALInConstructionAlgorithmArg &GDALVectorGridAbstractAlgorithm::AddNodataArg()
490 : {
491 268 : return AddArg("nodata", 0, _("Target nodata value"), &m_nodata)
492 268 : .SetDefault(m_nodata);
493 : }
494 :
495 : GDALVectorGridAlgorithmStandalone::~GDALVectorGridAlgorithmStandalone() =
496 : default;
497 :
498 : //! @endcond
|
