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 71 : GDALVectorCheckGeometryAlgorithm::GDALVectorCheckGeometryAlgorithm(
30 71 : bool standaloneStep)
31 : : GDALVectorPipelineStepAlgorithm(
32 : NAME, DESCRIPTION, HELP_URL,
33 0 : ConstructorOptions()
34 71 : .SetStandaloneStep(standaloneStep)
35 142 : .SetNoCreateEmptyLayersArgument(standaloneStep))
36 : {
37 : AddArg("include-field", 0,
38 : _("Fields from input layer to include in output (special values: "
39 : "ALL and NONE)"),
40 142 : &m_includeFields)
41 71 : .SetDefault("NONE");
42 :
43 : AddArg("include-valid", 0,
44 : _("Include valid inputs in output, with empty geometry"),
45 71 : &m_includeValid);
46 :
47 : AddArg("geometry-field", 0, _("Name of geometry field to check"),
48 71 : &m_geomField);
49 71 : }
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 37 : m_defn(OGRFeatureDefn::CreateFeatureDefn(
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 58 : m_srcGeomField(srcGeomField), m_skipValid(skipValid)
71 : {
72 29 : m_defn->Reference();
73 29 : m_defn->SetGeomType(wkbMultiPoint);
74 :
75 29 : if (!srcFieldIndices.empty())
76 : {
77 4 : const OGRFeatureDefn &srcDefn = *layer.GetLayerDefn();
78 4 : m_srcFieldMap.resize(srcDefn.GetFieldCount(), -1);
79 4 : int iDstField = 0;
80 13 : for (int iSrcField : srcFieldIndices)
81 : {
82 9 : m_defn->AddFieldDefn(srcDefn.GetFieldDefn(iSrcField));
83 9 : m_srcFieldMap[iSrcField] = iDstField++;
84 : }
85 : }
86 :
87 : auto poDescriptionFieldDefn =
88 58 : std::make_unique<OGRFieldDefn>(ERROR_DESCRIPTION_FIELD, OFTString);
89 29 : m_defn->AddFieldDefn(std::move(poDescriptionFieldDefn));
90 :
91 58 : m_defn->GetGeomFieldDefn(0)->SetSpatialRef(
92 29 : m_srcLayer.GetLayerDefn()
93 29 : ->GetGeomFieldDefn(m_srcGeomField)
94 29 : ->GetSpatialRef());
95 29 : }
96 :
97 : ~GDALInvalidLocationLayer() override;
98 :
99 13 : int TestCapability(const char *) const override
100 : {
101 13 : return false;
102 : }
103 :
104 1 : const char *GetDescription() const override
105 : {
106 1 : return GetName();
107 : }
108 :
109 187 : const OGRFeatureDefn *GetLayerDefn() const override
110 : {
111 187 : return m_defn;
112 : }
113 :
114 6 : std::unique_ptr<OGRFeature> CreateFeatureFromLastError() const
115 : {
116 6 : auto poErrorFeature = std::make_unique<OGRFeature>(m_defn);
117 :
118 12 : std::string msg = CPLGetLastErrorMsg();
119 :
120 : // Trim GEOS exception name
121 6 : const auto subMsgPos = msg.find(": ");
122 6 : if (subMsgPos != std::string::npos)
123 : {
124 6 : msg = msg.substr(subMsgPos + strlen(": "));
125 : }
126 :
127 : // Trim newline from end of GEOS exception message
128 6 : if (!msg.empty() && msg.back() == '\n')
129 : {
130 3 : msg.pop_back();
131 : }
132 :
133 6 : poErrorFeature->SetField(ERROR_DESCRIPTION_FIELD, msg.c_str());
134 :
135 6 : CPLErrorReset();
136 :
137 12 : return poErrorFeature;
138 : }
139 :
140 338 : void TranslateFeature(
141 : std::unique_ptr<OGRFeature> poSrcFeature,
142 : std::vector<std::unique_ptr<OGRFeature>> &apoOutputFeatures) override
143 : {
144 : const OGRGeometry *poGeom =
145 338 : poSrcFeature->GetGeomFieldRef(m_srcGeomField);
146 338 : std::unique_ptr<OGRFeature> poErrorFeature;
147 :
148 338 : if (poGeom)
149 : {
150 338 : if (poGeom->getDimension() < 1)
151 : {
152 1 : CPLErrorOnce(CE_Warning, CPLE_AppDefined,
153 : "Point geometry passed to 'gdal vector "
154 : "check-geometry'. Point geometries are "
155 : "always valid/simple.");
156 : }
157 : else
158 : {
159 337 : auto eType = wkbFlatten(poGeom->getGeometryType());
160 337 : GEOSGeometry *poGeosGeom = poGeom->exportToGEOS(m_geosContext);
161 :
162 337 : if (!poGeosGeom)
163 : {
164 : // Try to find a useful message / coordinate from
165 : // GEOS exception message.
166 6 : poErrorFeature = CreateFeatureFromLastError();
167 :
168 6 : if (eType == wkbPolygon)
169 : {
170 : const OGRLinearRing *poRing =
171 4 : poGeom->toPolygon()->getExteriorRing();
172 4 : if (poRing != nullptr && !poRing->IsEmpty())
173 : {
174 6 : auto poPoint = std::make_unique<OGRPoint>();
175 3 : poRing->StartPoint(poPoint.get());
176 : auto poMultiPoint =
177 3 : std::make_unique<OGRMultiPoint>();
178 3 : poMultiPoint->addGeometry(std::move(poPoint));
179 6 : poErrorFeature->SetGeometry(
180 3 : std::move(poMultiPoint));
181 : }
182 : else
183 : {
184 : // TODO get a point from somewhere else?
185 : }
186 : }
187 : }
188 : else
189 : {
190 331 : char *pszReason = nullptr;
191 331 : GEOSGeometry *location = nullptr;
192 331 : char ret = 1;
193 331 : bool warnAboutGeosVersion = false;
194 331 : bool checkedSimple = false;
195 :
196 331 : if (eType == wkbPolygon || eType == wkbMultiPolygon ||
197 10 : eType == wkbCurvePolygon || eType == wkbMultiSurface ||
198 : eType == wkbGeometryCollection)
199 : {
200 323 : ret = GEOSisValidDetail_r(m_geosContext, poGeosGeom, 0,
201 : &pszReason, &location);
202 : }
203 :
204 331 : if (eType == wkbLineString || eType == wkbMultiLineString ||
205 325 : eType == wkbCircularString ||
206 323 : eType == wkbCompoundCurve ||
207 149 : (ret == 1 && eType == wkbGeometryCollection))
208 : {
209 10 : checkedSimple = true;
210 : #if GEOS_VERSION_MAJOR > 3 || \
211 : (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 14)
212 10 : ret = GEOSisSimpleDetail_r(m_geosContext, poGeosGeom, 1,
213 : &location);
214 : #else
215 : ret = GEOSisSimple_r(m_geosContext, poGeosGeom);
216 : warnAboutGeosVersion = true;
217 : #endif
218 : }
219 :
220 331 : GEOSGeom_destroy_r(m_geosContext, poGeosGeom);
221 331 : if (ret == 0)
222 : {
223 182 : if (warnAboutGeosVersion)
224 : {
225 0 : CPLErrorOnce(
226 : CE_Warning, CPLE_AppDefined,
227 : "Detected a non-simple linear geometry, but "
228 : "cannot output self-intersection points "
229 : "because GEOS library version is < 3.14.");
230 : }
231 :
232 182 : poErrorFeature = std::make_unique<OGRFeature>(m_defn);
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);
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 : OGRFeatureDefn *const 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 : m_defn->Release();
312 29 : finishGEOS_r(m_geosContext);
313 58 : }
314 :
315 : #endif
316 :
317 29 : bool GDALVectorCheckGeometryAlgorithm::RunStep(GDALPipelineStepRunContext &)
318 : {
319 : #ifdef HAVE_GEOS
320 29 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
321 29 : CPLAssert(poSrcDS);
322 29 : CPLAssert(m_outputDataset.GetName().empty());
323 29 : CPLAssert(!m_outputDataset.GetDatasetRef());
324 :
325 29 : const bool bSingleLayerOutput = m_inputLayerNames.empty()
326 31 : ? poSrcDS->GetLayerCount() == 1
327 2 : : m_inputLayerNames.size() == 1;
328 :
329 58 : auto outDS = std::make_unique<GDALVectorPipelineOutputDataset>(*poSrcDS);
330 59 : for (auto &&poSrcLayer : poSrcDS->GetLayers())
331 : {
332 35 : if (m_inputLayerNames.empty() ||
333 0 : std::find(m_inputLayerNames.begin(), m_inputLayerNames.end(),
334 35 : poSrcLayer->GetDescription()) != m_inputLayerNames.end())
335 : {
336 33 : const auto poSrcLayerDefn = poSrcLayer->GetLayerDefn();
337 33 : if (poSrcLayerDefn->GetGeomFieldCount() == 0)
338 : {
339 2 : if (m_inputLayerNames.empty())
340 1 : continue;
341 1 : ReportError(CE_Failure, CPLE_AppDefined,
342 : "Specified layer '%s' has no geometry field",
343 1 : poSrcLayer->GetDescription());
344 3 : return false;
345 : }
346 :
347 : const int geomFieldIndex =
348 31 : m_geomField.empty()
349 31 : ? 0
350 2 : : poSrcLayerDefn->GetGeomFieldIndex(m_geomField.c_str());
351 :
352 31 : if (geomFieldIndex == -1)
353 : {
354 1 : ReportError(CE_Failure, CPLE_AppDefined,
355 : "Specified geometry field '%s' does not exist in "
356 : "layer '%s'",
357 1 : m_geomField.c_str(), poSrcLayer->GetDescription());
358 1 : return false;
359 : }
360 :
361 30 : std::vector<int> includeFieldIndices;
362 30 : if (!GetFieldIndices(m_includeFields,
363 : OGRLayer::ToHandle(poSrcLayer),
364 : includeFieldIndices))
365 : {
366 1 : return false;
367 : }
368 :
369 58 : outDS->AddLayer(*poSrcLayer,
370 29 : std::make_unique<GDALInvalidLocationLayer>(
371 : *poSrcLayer, includeFieldIndices,
372 : bSingleLayerOutput, geomFieldIndex,
373 58 : !m_includeValid));
374 : }
375 : }
376 :
377 26 : m_outputDataset.Set(std::move(outDS));
378 :
379 26 : return true;
380 : #else
381 : ReportError(CE_Failure, CPLE_AppDefined,
382 : "%s requires GDAL to be built against the GEOS library.", NAME);
383 : return false;
384 : #endif
385 : }
386 :
387 : GDALVectorCheckGeometryAlgorithmStandalone::
388 : ~GDALVectorCheckGeometryAlgorithmStandalone() = default;
389 :
390 : //! @endcond
|
