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 109 : GDALVectorCheckGeometryAlgorithm::GDALVectorCheckGeometryAlgorithm(
30 109 : bool standaloneStep)
31 : : GDALVectorPipelineStepAlgorithm(
32 : NAME, DESCRIPTION, HELP_URL,
33 0 : ConstructorOptions()
34 109 : .SetStandaloneStep(standaloneStep)
35 218 : .SetNoCreateEmptyLayersArgument(standaloneStep))
36 : {
37 : AddArg("include-field", 0,
38 : _("Fields from input layer to include in output (special values: "
39 : "ALL and NONE)"),
40 218 : &m_includeFields)
41 109 : .SetDefault("NONE");
42 :
43 : AddArg("include-valid", 0,
44 : _("Include valid inputs in output, with empty geometry"),
45 109 : &m_includeValid);
46 :
47 : AddArg("geometry-field", 0, _("Name of geometry field to check"),
48 109 : &m_geomField);
49 109 : }
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 : bool 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 0 : 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 :
298 676 : return true;
299 : }
300 :
301 : CPL_DISALLOW_COPY_ASSIGN(GDALInvalidLocationLayer)
302 :
303 : private:
304 : std::vector<int> m_srcFieldMap{};
305 : const OGRFeatureDefnRefCountedPtr m_defn;
306 : const GEOSContextHandle_t m_geosContext;
307 : const int m_srcGeomField;
308 : const bool m_skipValid;
309 : };
310 :
311 58 : GDALInvalidLocationLayer::~GDALInvalidLocationLayer()
312 : {
313 29 : finishGEOS_r(m_geosContext);
314 58 : }
315 :
316 : #endif
317 :
318 29 : bool GDALVectorCheckGeometryAlgorithm::RunStep(GDALPipelineStepRunContext &)
319 : {
320 : #ifdef HAVE_GEOS
321 29 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
322 29 : CPLAssert(poSrcDS);
323 29 : CPLAssert(m_outputDataset.GetName().empty());
324 29 : CPLAssert(!m_outputDataset.GetDatasetRef());
325 :
326 29 : const bool bSingleLayerOutput = m_inputLayerNames.empty()
327 31 : ? poSrcDS->GetLayerCount() == 1
328 2 : : m_inputLayerNames.size() == 1;
329 :
330 58 : auto outDS = std::make_unique<GDALVectorPipelineOutputDataset>(*poSrcDS);
331 59 : for (auto &&poSrcLayer : poSrcDS->GetLayers())
332 : {
333 35 : if (m_inputLayerNames.empty() ||
334 0 : std::find(m_inputLayerNames.begin(), m_inputLayerNames.end(),
335 35 : poSrcLayer->GetDescription()) != m_inputLayerNames.end())
336 : {
337 33 : const auto poSrcLayerDefn = poSrcLayer->GetLayerDefn();
338 33 : if (poSrcLayerDefn->GetGeomFieldCount() == 0)
339 : {
340 2 : if (m_inputLayerNames.empty())
341 1 : continue;
342 1 : ReportError(CE_Failure, CPLE_AppDefined,
343 : "Specified layer '%s' has no geometry field",
344 1 : poSrcLayer->GetDescription());
345 3 : return false;
346 : }
347 :
348 : const int geomFieldIndex =
349 31 : m_geomField.empty()
350 31 : ? 0
351 2 : : poSrcLayerDefn->GetGeomFieldIndex(m_geomField.c_str());
352 :
353 31 : if (geomFieldIndex == -1)
354 : {
355 1 : ReportError(CE_Failure, CPLE_AppDefined,
356 : "Specified geometry field '%s' does not exist in "
357 : "layer '%s'",
358 1 : m_geomField.c_str(), poSrcLayer->GetDescription());
359 1 : return false;
360 : }
361 :
362 30 : std::vector<int> includeFieldIndices;
363 30 : if (!GetFieldIndices(m_includeFields,
364 : OGRLayer::ToHandle(poSrcLayer),
365 : includeFieldIndices))
366 : {
367 1 : return false;
368 : }
369 :
370 58 : outDS->AddLayer(*poSrcLayer,
371 29 : std::make_unique<GDALInvalidLocationLayer>(
372 : *poSrcLayer, includeFieldIndices,
373 : bSingleLayerOutput, geomFieldIndex,
374 58 : !m_includeValid));
375 : }
376 : }
377 :
378 26 : m_outputDataset.Set(std::move(outDS));
379 :
380 26 : return true;
381 : #else
382 : ReportError(CE_Failure, CPLE_AppDefined,
383 : "%s requires GDAL to be built against the GEOS library.", NAME);
384 : return false;
385 : #endif
386 : }
387 :
388 : GDALVectorCheckGeometryAlgorithmStandalone::
389 : ~GDALVectorCheckGeometryAlgorithmStandalone() = default;
390 :
391 : //! @endcond
|
