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 131 : GDALVectorCheckGeometryAlgorithm::GDALVectorCheckGeometryAlgorithm(
30 131 : bool standaloneStep)
31 : : GDALVectorPipelineStepAlgorithm(
32 : NAME, DESCRIPTION, HELP_URL,
33 0 : ConstructorOptions()
34 131 : .SetStandaloneStep(standaloneStep)
35 262 : .SetNoCreateEmptyLayersArgument(standaloneStep))
36 : {
37 : AddArg("include-field", 0,
38 : _("Fields from input layer to include in output (special values: "
39 : "ALL and NONE)"),
40 262 : &m_includeFields)
41 131 : .SetDefault("NONE");
42 :
43 : AddArg("include-valid", 0,
44 : _("Include valid inputs in output, with empty geometry"),
45 131 : &m_includeValid);
46 :
47 : AddArg("geometry-field", 0, _("Name of geometry field to check"),
48 131 : &m_geomField);
49 131 : }
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 42 : GDALInvalidLocationLayer(OGRLayer &layer,
60 : const std::vector<int> &srcFieldIndices,
61 : bool bSingleLayerOutput, int srcGeomField,
62 : bool skipValid)
63 42 : : GDALVectorPipelineOutputLayer(layer),
64 : m_defn(OGRFeatureDefnRefCountedPtr::makeInstance(
65 : bSingleLayerOutput ? "error_location"
66 50 : : std::string("error_location_")
67 8 : .append(layer.GetDescription())
68 8 : .c_str())),
69 42 : m_geosContext(OGRGeometry::createGEOSContext()),
70 58 : m_srcGeomField(srcGeomField), m_skipValid(skipValid)
71 : {
72 42 : m_defn->SetGeomType(wkbMultiPoint);
73 :
74 42 : 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 84 : std::make_unique<OGRFieldDefn>(ERROR_DESCRIPTION_FIELD, OFTString);
88 42 : m_defn->AddFieldDefn(std::move(poDescriptionFieldDefn));
89 :
90 84 : m_defn->GetGeomFieldDefn(0)->SetSpatialRef(
91 42 : m_srcLayer.GetLayerDefn()
92 42 : ->GetGeomFieldDefn(m_srcGeomField)
93 42 : ->GetSpatialRef());
94 42 : }
95 :
96 : ~GDALInvalidLocationLayer() override;
97 :
98 13 : bool TestCapability(const char *) const override
99 : {
100 13 : return false;
101 : }
102 :
103 2 : const char *GetDescription() const override
104 : {
105 2 : return GetName();
106 : }
107 :
108 190 : const OGRFeatureDefn *GetLayerDefn() const override
109 : {
110 190 : 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 356 : bool TranslateFeature(
140 : std::unique_ptr<OGRFeature> poSrcFeature,
141 : std::vector<std::unique_ptr<OGRFeature>> &apoOutputFeatures) override
142 : {
143 : const OGRGeometry *poGeom =
144 356 : poSrcFeature->GetGeomFieldRef(m_srcGeomField);
145 0 : std::unique_ptr<OGRFeature> poErrorFeature;
146 :
147 356 : if (poGeom)
148 : {
149 356 : auto eType = wkbFlatten(poGeom->getGeometryType());
150 356 : GEOSGeometry *poGeosGeom = poGeom->exportToGEOS(m_geosContext);
151 :
152 356 : if (!poGeosGeom)
153 : {
154 : // Try to find a useful message / coordinate from
155 : // GEOS exception message.
156 6 : poErrorFeature = CreateFeatureFromLastError();
157 :
158 6 : if (eType == wkbPolygon)
159 : {
160 : const OGRLinearRing *poRing =
161 4 : poGeom->toPolygon()->getExteriorRing();
162 4 : if (poRing != nullptr && !poRing->IsEmpty())
163 : {
164 6 : auto poPoint = std::make_unique<OGRPoint>();
165 3 : poRing->StartPoint(poPoint.get());
166 3 : auto poMultiPoint = std::make_unique<OGRMultiPoint>();
167 3 : poMultiPoint->addGeometry(std::move(poPoint));
168 3 : poErrorFeature->SetGeometry(std::move(poMultiPoint));
169 : }
170 : else
171 : {
172 : // TODO get a point from somewhere else?
173 : }
174 : }
175 : }
176 : else
177 : {
178 350 : char *pszReason = nullptr;
179 350 : GEOSGeometry *location = nullptr;
180 350 : char ret = 1;
181 350 : bool warnAboutGeosVersion = false;
182 350 : bool checkedSimple = false;
183 :
184 : // check all geometry types for validity; IsSimple does not
185 : // detect NaN/Inf coordinates
186 350 : ret = GEOSisValidDetail_r(m_geosContext, poGeosGeom, 0,
187 : &pszReason, &location);
188 :
189 350 : if (ret == 1 &&
190 158 : (eType == wkbLineString || eType == wkbMultiLineString ||
191 158 : eType == wkbCircularString || eType == wkbCompoundCurve ||
192 : eType == wkbGeometryCollection))
193 : {
194 12 : checkedSimple = true;
195 : #if GEOS_VERSION_MAJOR > 3 || \
196 : (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 14)
197 12 : ret = GEOSisSimpleDetail_r(m_geosContext, poGeosGeom, 1,
198 : &location);
199 : #else
200 : ret = GEOSisSimple_r(m_geosContext, poGeosGeom);
201 : warnAboutGeosVersion = true;
202 : #endif
203 : }
204 :
205 350 : GEOSGeom_destroy_r(m_geosContext, poGeosGeom);
206 350 : if (ret == 0)
207 : {
208 192 : if (warnAboutGeosVersion)
209 : {
210 0 : CPLErrorOnce(
211 : CE_Warning, CPLE_AppDefined,
212 : "Detected a non-simple linear geometry, but "
213 : "cannot output self-intersection points "
214 : "because GEOS library version is < 3.14.");
215 : }
216 :
217 192 : poErrorFeature = std::make_unique<OGRFeature>(m_defn.get());
218 192 : if (pszReason == nullptr)
219 : {
220 8 : if (checkedSimple)
221 : {
222 8 : poErrorFeature->SetField(ERROR_DESCRIPTION_FIELD,
223 : "self-intersection");
224 : }
225 : }
226 : else
227 : {
228 184 : poErrorFeature->SetField(ERROR_DESCRIPTION_FIELD,
229 : pszReason);
230 184 : GEOSFree_r(m_geosContext, pszReason);
231 : }
232 :
233 192 : if (location != nullptr)
234 : {
235 : std::unique_ptr<OGRGeometry> poErrorGeom(
236 192 : OGRGeometryFactory::createFromGEOS(m_geosContext,
237 192 : location));
238 192 : GEOSGeom_destroy_r(m_geosContext, location);
239 :
240 192 : if (poErrorGeom->getGeometryType() == wkbPoint)
241 : {
242 : auto poMultiPoint =
243 376 : std::make_unique<OGRMultiPoint>();
244 188 : poMultiPoint->addGeometry(std::move(poErrorGeom));
245 188 : poErrorGeom = std::move(poMultiPoint);
246 : }
247 :
248 384 : poErrorGeom->assignSpatialReference(
249 384 : m_srcLayer.GetLayerDefn()
250 192 : ->GetGeomFieldDefn(m_srcGeomField)
251 384 : ->GetSpatialRef());
252 :
253 192 : poErrorFeature->SetGeometry(std::move(poErrorGeom));
254 : }
255 : }
256 158 : else if (ret == 2)
257 : {
258 0 : poErrorFeature = CreateFeatureFromLastError();
259 : }
260 : }
261 : }
262 :
263 356 : if (!poErrorFeature && !m_skipValid)
264 : {
265 6 : poErrorFeature = std::make_unique<OGRFeature>(m_defn.get());
266 : // TODO Set geometry to POINT EMPTY ?
267 : }
268 :
269 356 : if (poErrorFeature)
270 : {
271 204 : if (!m_srcFieldMap.empty())
272 : {
273 8 : poErrorFeature->SetFieldsFrom(
274 4 : poSrcFeature.get(), m_srcFieldMap.data(), false, false);
275 : }
276 204 : poErrorFeature->SetFID(poSrcFeature->GetFID());
277 204 : if (PassesFilters(poErrorFeature.get()))
278 203 : apoOutputFeatures.push_back(std::move(poErrorFeature));
279 : }
280 :
281 712 : return true;
282 : }
283 :
284 : CPL_DISALLOW_COPY_ASSIGN(GDALInvalidLocationLayer)
285 :
286 : private:
287 : std::vector<int> m_srcFieldMap{};
288 : const OGRFeatureDefnRefCountedPtr m_defn;
289 : const GEOSContextHandle_t m_geosContext;
290 : const int m_srcGeomField;
291 : const bool m_skipValid;
292 : };
293 :
294 84 : GDALInvalidLocationLayer::~GDALInvalidLocationLayer()
295 : {
296 42 : finishGEOS_r(m_geosContext);
297 84 : }
298 :
299 : #endif
300 :
301 42 : bool GDALVectorCheckGeometryAlgorithm::RunStep(GDALPipelineStepRunContext &)
302 : {
303 : #ifdef HAVE_GEOS
304 42 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
305 42 : CPLAssert(poSrcDS);
306 42 : CPLAssert(m_outputDataset.GetName().empty());
307 42 : CPLAssert(!m_outputDataset.GetDatasetRef());
308 :
309 42 : const bool bSingleLayerOutput = m_inputLayerNames.empty()
310 44 : ? poSrcDS->GetLayerCount() == 1
311 2 : : m_inputLayerNames.size() == 1;
312 :
313 84 : auto outDS = std::make_unique<GDALVectorPipelineOutputDataset>(*poSrcDS);
314 85 : for (auto &&poSrcLayer : poSrcDS->GetLayers())
315 : {
316 48 : if (m_inputLayerNames.empty() ||
317 0 : std::find(m_inputLayerNames.begin(), m_inputLayerNames.end(),
318 48 : poSrcLayer->GetDescription()) != m_inputLayerNames.end())
319 : {
320 46 : const auto poSrcLayerDefn = poSrcLayer->GetLayerDefn();
321 46 : if (poSrcLayerDefn->GetGeomFieldCount() == 0)
322 : {
323 2 : if (m_inputLayerNames.empty())
324 1 : continue;
325 1 : ReportError(CE_Failure, CPLE_AppDefined,
326 : "Specified layer '%s' has no geometry field",
327 1 : poSrcLayer->GetDescription());
328 3 : return false;
329 : }
330 :
331 : const int geomFieldIndex =
332 44 : m_geomField.empty()
333 44 : ? 0
334 2 : : poSrcLayerDefn->GetGeomFieldIndex(m_geomField.c_str());
335 :
336 44 : if (geomFieldIndex == -1)
337 : {
338 1 : ReportError(CE_Failure, CPLE_AppDefined,
339 : "Specified geometry field '%s' does not exist in "
340 : "layer '%s'",
341 1 : m_geomField.c_str(), poSrcLayer->GetDescription());
342 1 : return false;
343 : }
344 :
345 43 : std::vector<int> includeFieldIndices;
346 43 : if (!GetFieldIndices(m_includeFields,
347 : OGRLayer::ToHandle(poSrcLayer),
348 : includeFieldIndices))
349 : {
350 1 : return false;
351 : }
352 :
353 84 : outDS->AddLayer(*poSrcLayer,
354 42 : std::make_unique<GDALInvalidLocationLayer>(
355 : *poSrcLayer, includeFieldIndices,
356 : bSingleLayerOutput, geomFieldIndex,
357 84 : !m_includeValid));
358 : }
359 : }
360 :
361 39 : m_outputDataset.Set(std::move(outDS));
362 :
363 39 : return true;
364 : #else
365 : ReportError(CE_Failure, CPLE_AppDefined,
366 : "%s requires GDAL to be built against the GEOS library.", NAME);
367 : return false;
368 : #endif
369 : }
370 :
371 : GDALVectorCheckGeometryAlgorithmStandalone::
372 : ~GDALVectorCheckGeometryAlgorithmStandalone() = default;
373 :
374 : //! @endcond
|
