Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: "gdal vector check-geometry" subcommand
5 : * Author: Daniel Baston
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2025, ISciences LLC
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalalg_vector_check_geometry.h"
14 :
15 : #include "cpl_error.h"
16 : #include "gdal_priv.h"
17 : #include "gdalalg_vector_geom.h"
18 : #include "ogr_geometry.h"
19 : #include "ogr_geos.h"
20 :
21 : #include <cinttypes>
22 :
23 : #ifndef _
24 : #define _(x) (x)
25 : #endif
26 :
27 : //! @cond Doxygen_Suppress
28 :
29 107 : GDALVectorCheckGeometryAlgorithm::GDALVectorCheckGeometryAlgorithm(
30 107 : bool standaloneStep)
31 : : GDALVectorPipelineStepAlgorithm(
32 : NAME, DESCRIPTION, HELP_URL,
33 0 : ConstructorOptions()
34 107 : .SetStandaloneStep(standaloneStep)
35 214 : .SetNoCreateEmptyLayersArgument(standaloneStep))
36 : {
37 : AddArg("include-field", 0,
38 : _("Fields from input layer to include in output (special values: "
39 : "ALL and NONE)"),
40 214 : &m_includeFields)
41 107 : .SetDefault("NONE");
42 :
43 : AddArg("include-valid", 0,
44 : _("Include valid inputs in output, with empty geometry"),
45 107 : &m_includeValid);
46 :
47 : AddArg("geometry-field", 0, _("Name of geometry field to check"),
48 107 : &m_geomField);
49 107 : }
50 :
51 : #ifdef HAVE_GEOS
52 :
53 : class GDALInvalidLocationLayer final : public GDALVectorPipelineOutputLayer
54 : {
55 : private:
56 : static constexpr const char *ERROR_DESCRIPTION_FIELD = "error";
57 :
58 : public:
59 29 : GDALInvalidLocationLayer(OGRLayer &layer,
60 : const std::vector<int> &srcFieldIndices,
61 : bool bSingleLayerOutput, int srcGeomField,
62 : bool skipValid)
63 29 : : GDALVectorPipelineOutputLayer(layer),
64 : m_defn(OGRFeatureDefnRefCountedPtr::makeInstance(
65 : bSingleLayerOutput ? "error_location"
66 37 : : std::string("error_location_")
67 8 : .append(layer.GetDescription())
68 8 : .c_str())),
69 29 : m_geosContext(OGRGeometry::createGEOSContext()),
70 45 : m_srcGeomField(srcGeomField), m_skipValid(skipValid)
71 : {
72 29 : m_defn->SetGeomType(wkbMultiPoint);
73 :
74 29 : if (!srcFieldIndices.empty())
75 : {
76 4 : const OGRFeatureDefn &srcDefn = *layer.GetLayerDefn();
77 4 : m_srcFieldMap.resize(srcDefn.GetFieldCount(), -1);
78 4 : int iDstField = 0;
79 13 : for (int iSrcField : srcFieldIndices)
80 : {
81 9 : m_defn->AddFieldDefn(srcDefn.GetFieldDefn(iSrcField));
82 9 : m_srcFieldMap[iSrcField] = iDstField++;
83 : }
84 : }
85 :
86 : auto poDescriptionFieldDefn =
87 58 : std::make_unique<OGRFieldDefn>(ERROR_DESCRIPTION_FIELD, OFTString);
88 29 : m_defn->AddFieldDefn(std::move(poDescriptionFieldDefn));
89 :
90 58 : m_defn->GetGeomFieldDefn(0)->SetSpatialRef(
91 29 : m_srcLayer.GetLayerDefn()
92 29 : ->GetGeomFieldDefn(m_srcGeomField)
93 29 : ->GetSpatialRef());
94 29 : }
95 :
96 : ~GDALInvalidLocationLayer() override;
97 :
98 13 : int TestCapability(const char *) const override
99 : {
100 13 : return false;
101 : }
102 :
103 1 : const char *GetDescription() const override
104 : {
105 1 : return GetName();
106 : }
107 :
108 187 : const OGRFeatureDefn *GetLayerDefn() const override
109 : {
110 187 : return m_defn.get();
111 : }
112 :
113 6 : std::unique_ptr<OGRFeature> CreateFeatureFromLastError() const
114 : {
115 6 : auto poErrorFeature = std::make_unique<OGRFeature>(m_defn.get());
116 :
117 12 : std::string msg = CPLGetLastErrorMsg();
118 :
119 : // Trim GEOS exception name
120 6 : const auto subMsgPos = msg.find(": ");
121 6 : if (subMsgPos != std::string::npos)
122 : {
123 6 : msg = msg.substr(subMsgPos + strlen(": "));
124 : }
125 :
126 : // Trim newline from end of GEOS exception message
127 6 : if (!msg.empty() && msg.back() == '\n')
128 : {
129 3 : msg.pop_back();
130 : }
131 :
132 6 : poErrorFeature->SetField(ERROR_DESCRIPTION_FIELD, msg.c_str());
133 :
134 6 : CPLErrorReset();
135 :
136 12 : return poErrorFeature;
137 : }
138 :
139 338 : void TranslateFeature(
140 : std::unique_ptr<OGRFeature> poSrcFeature,
141 : std::vector<std::unique_ptr<OGRFeature>> &apoOutputFeatures) override
142 : {
143 : const OGRGeometry *poGeom =
144 338 : poSrcFeature->GetGeomFieldRef(m_srcGeomField);
145 338 : std::unique_ptr<OGRFeature> poErrorFeature;
146 :
147 338 : if (poGeom)
148 : {
149 338 : if (poGeom->getDimension() < 1)
150 : {
151 1 : CPLErrorOnce(CE_Warning, CPLE_AppDefined,
152 : "Point geometry passed to 'gdal vector "
153 : "check-geometry'. Point geometries are "
154 : "always valid/simple.");
155 : }
156 : else
157 : {
158 337 : auto eType = wkbFlatten(poGeom->getGeometryType());
159 337 : GEOSGeometry *poGeosGeom = poGeom->exportToGEOS(m_geosContext);
160 :
161 337 : if (!poGeosGeom)
162 : {
163 : // Try to find a useful message / coordinate from
164 : // GEOS exception message.
165 6 : poErrorFeature = CreateFeatureFromLastError();
166 :
167 6 : if (eType == wkbPolygon)
168 : {
169 : const OGRLinearRing *poRing =
170 4 : poGeom->toPolygon()->getExteriorRing();
171 4 : if (poRing != nullptr && !poRing->IsEmpty())
172 : {
173 6 : auto poPoint = std::make_unique<OGRPoint>();
174 3 : poRing->StartPoint(poPoint.get());
175 : auto poMultiPoint =
176 3 : std::make_unique<OGRMultiPoint>();
177 3 : poMultiPoint->addGeometry(std::move(poPoint));
178 6 : poErrorFeature->SetGeometry(
179 3 : std::move(poMultiPoint));
180 : }
181 : else
182 : {
183 : // TODO get a point from somewhere else?
184 : }
185 : }
186 : }
187 : else
188 : {
189 331 : char *pszReason = nullptr;
190 331 : GEOSGeometry *location = nullptr;
191 331 : char ret = 1;
192 331 : bool warnAboutGeosVersion = false;
193 331 : bool checkedSimple = false;
194 :
195 331 : if (eType == wkbPolygon || eType == wkbMultiPolygon ||
196 10 : eType == wkbCurvePolygon || eType == wkbMultiSurface ||
197 : eType == wkbGeometryCollection)
198 : {
199 323 : ret = GEOSisValidDetail_r(m_geosContext, poGeosGeom, 0,
200 : &pszReason, &location);
201 : }
202 :
203 331 : if (eType == wkbLineString || eType == wkbMultiLineString ||
204 325 : eType == wkbCircularString ||
205 323 : eType == wkbCompoundCurve ||
206 149 : (ret == 1 && eType == wkbGeometryCollection))
207 : {
208 10 : checkedSimple = true;
209 : #if GEOS_VERSION_MAJOR > 3 || \
210 : (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 14)
211 10 : ret = GEOSisSimpleDetail_r(m_geosContext, poGeosGeom, 1,
212 : &location);
213 : #else
214 : ret = GEOSisSimple_r(m_geosContext, poGeosGeom);
215 : warnAboutGeosVersion = true;
216 : #endif
217 : }
218 :
219 331 : GEOSGeom_destroy_r(m_geosContext, poGeosGeom);
220 331 : if (ret == 0)
221 : {
222 182 : if (warnAboutGeosVersion)
223 : {
224 0 : CPLErrorOnce(
225 : CE_Warning, CPLE_AppDefined,
226 : "Detected a non-simple linear geometry, but "
227 : "cannot output self-intersection points "
228 : "because GEOS library version is < 3.14.");
229 : }
230 :
231 : poErrorFeature =
232 182 : std::make_unique<OGRFeature>(m_defn.get());
233 182 : if (pszReason == nullptr)
234 : {
235 8 : if (checkedSimple)
236 : {
237 8 : poErrorFeature->SetField(
238 : ERROR_DESCRIPTION_FIELD,
239 : "self-intersection");
240 : }
241 : }
242 : else
243 : {
244 174 : poErrorFeature->SetField(ERROR_DESCRIPTION_FIELD,
245 : pszReason);
246 174 : GEOSFree_r(m_geosContext, pszReason);
247 : }
248 :
249 182 : if (location != nullptr)
250 : {
251 : std::unique_ptr<OGRGeometry> poErrorGeom(
252 : OGRGeometryFactory::createFromGEOS(
253 182 : m_geosContext, location));
254 182 : GEOSGeom_destroy_r(m_geosContext, location);
255 :
256 182 : if (poErrorGeom->getGeometryType() == wkbPoint)
257 : {
258 : auto poMultiPoint =
259 356 : std::make_unique<OGRMultiPoint>();
260 356 : poMultiPoint->addGeometry(
261 178 : std::move(poErrorGeom));
262 178 : poErrorGeom = std::move(poMultiPoint);
263 : }
264 :
265 364 : poErrorGeom->assignSpatialReference(
266 364 : m_srcLayer.GetLayerDefn()
267 182 : ->GetGeomFieldDefn(m_srcGeomField)
268 364 : ->GetSpatialRef());
269 :
270 182 : poErrorFeature->SetGeometry(std::move(poErrorGeom));
271 : }
272 : }
273 149 : else if (ret == 2)
274 : {
275 0 : poErrorFeature = CreateFeatureFromLastError();
276 : }
277 : }
278 : }
279 : }
280 :
281 338 : if (!poErrorFeature && !m_skipValid)
282 : {
283 6 : poErrorFeature = std::make_unique<OGRFeature>(m_defn.get());
284 : // TODO Set geometry to POINT EMPTY ?
285 : }
286 :
287 338 : if (poErrorFeature)
288 : {
289 194 : if (!m_srcFieldMap.empty())
290 : {
291 8 : poErrorFeature->SetFieldsFrom(
292 4 : poSrcFeature.get(), m_srcFieldMap.data(), false, false);
293 : }
294 194 : poErrorFeature->SetFID(poSrcFeature->GetFID());
295 194 : apoOutputFeatures.push_back(std::move(poErrorFeature));
296 : }
297 338 : }
298 :
299 : CPL_DISALLOW_COPY_ASSIGN(GDALInvalidLocationLayer)
300 :
301 : private:
302 : std::vector<int> m_srcFieldMap{};
303 : const OGRFeatureDefnRefCountedPtr m_defn;
304 : const GEOSContextHandle_t m_geosContext;
305 : const int m_srcGeomField;
306 : const bool m_skipValid;
307 : };
308 :
309 58 : GDALInvalidLocationLayer::~GDALInvalidLocationLayer()
310 : {
311 29 : finishGEOS_r(m_geosContext);
312 58 : }
313 :
314 : #endif
315 :
316 29 : bool GDALVectorCheckGeometryAlgorithm::RunStep(GDALPipelineStepRunContext &)
317 : {
318 : #ifdef HAVE_GEOS
319 29 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
320 29 : CPLAssert(poSrcDS);
321 29 : CPLAssert(m_outputDataset.GetName().empty());
322 29 : CPLAssert(!m_outputDataset.GetDatasetRef());
323 :
324 29 : const bool bSingleLayerOutput = m_inputLayerNames.empty()
325 31 : ? poSrcDS->GetLayerCount() == 1
326 2 : : m_inputLayerNames.size() == 1;
327 :
328 58 : auto outDS = std::make_unique<GDALVectorPipelineOutputDataset>(*poSrcDS);
329 59 : for (auto &&poSrcLayer : poSrcDS->GetLayers())
330 : {
331 35 : if (m_inputLayerNames.empty() ||
332 0 : std::find(m_inputLayerNames.begin(), m_inputLayerNames.end(),
333 35 : poSrcLayer->GetDescription()) != m_inputLayerNames.end())
334 : {
335 33 : const auto poSrcLayerDefn = poSrcLayer->GetLayerDefn();
336 33 : if (poSrcLayerDefn->GetGeomFieldCount() == 0)
337 : {
338 2 : if (m_inputLayerNames.empty())
339 1 : continue;
340 1 : ReportError(CE_Failure, CPLE_AppDefined,
341 : "Specified layer '%s' has no geometry field",
342 1 : poSrcLayer->GetDescription());
343 3 : return false;
344 : }
345 :
346 : const int geomFieldIndex =
347 31 : m_geomField.empty()
348 31 : ? 0
349 2 : : poSrcLayerDefn->GetGeomFieldIndex(m_geomField.c_str());
350 :
351 31 : if (geomFieldIndex == -1)
352 : {
353 1 : ReportError(CE_Failure, CPLE_AppDefined,
354 : "Specified geometry field '%s' does not exist in "
355 : "layer '%s'",
356 1 : m_geomField.c_str(), poSrcLayer->GetDescription());
357 1 : return false;
358 : }
359 :
360 30 : std::vector<int> includeFieldIndices;
361 30 : if (!GetFieldIndices(m_includeFields,
362 : OGRLayer::ToHandle(poSrcLayer),
363 : includeFieldIndices))
364 : {
365 1 : return false;
366 : }
367 :
368 58 : outDS->AddLayer(*poSrcLayer,
369 29 : std::make_unique<GDALInvalidLocationLayer>(
370 : *poSrcLayer, includeFieldIndices,
371 : bSingleLayerOutput, geomFieldIndex,
372 58 : !m_includeValid));
373 : }
374 : }
375 :
376 26 : m_outputDataset.Set(std::move(outDS));
377 :
378 26 : return true;
379 : #else
380 : ReportError(CE_Failure, CPLE_AppDefined,
381 : "%s requires GDAL to be built against the GEOS library.", NAME);
382 : return false;
383 : #endif
384 : }
385 :
386 : GDALVectorCheckGeometryAlgorithmStandalone::
387 : ~GDALVectorCheckGeometryAlgorithmStandalone() = default;
388 :
389 : //! @endcond
|
