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