Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: "gdal vector clean-coverage" 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_clean_coverage.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 : #include "ogrsf_frmts.h"
21 :
22 : #include <cinttypes>
23 :
24 : #ifndef _
25 : #define _(x) (x)
26 : #endif
27 :
28 : //! @cond Doxygen_Suppress
29 :
30 36 : GDALVectorCleanCoverageAlgorithm::GDALVectorCleanCoverageAlgorithm(
31 36 : bool standaloneStep)
32 : : GDALVectorPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
33 36 : standaloneStep)
34 : {
35 36 : AddActiveLayerArg(&m_activeLayer);
36 : AddArg("snapping-distance", 0, _("Distance tolerance for snapping nodes"),
37 72 : &m_opts.snappingTolerance)
38 36 : .SetMinValueIncluded(0);
39 :
40 : AddArg("merge-strategy", 0,
41 : _("Algorithm to assign overlaps to neighboring polygons"),
42 72 : &m_opts.mergeStrategy)
43 180 : .SetChoices({"longest-border", "max-area", "min-area", "min-index"});
44 :
45 : AddArg("maximum-gap-width", 0, _("Maximum width of a gap to be closed"),
46 72 : &m_opts.maximumGapWidth)
47 36 : .SetMinValueIncluded(0);
48 36 : }
49 :
50 : #if defined HAVE_GEOS && \
51 : (GEOS_VERSION_MAJOR > 3 || \
52 : (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 14))
53 :
54 : class GDALVectorCleanCoverageOutputDataset
55 : : public GDALVectorNonStreamingAlgorithmDataset
56 : {
57 : public:
58 14 : GDALVectorCleanCoverageOutputDataset(
59 : const GDALVectorCleanCoverageAlgorithm::Options &opts)
60 14 : : m_opts(opts)
61 : {
62 14 : m_poGeosContext = OGRGeometry::createGEOSContext();
63 14 : }
64 :
65 : ~GDALVectorCleanCoverageOutputDataset() override;
66 :
67 13 : GEOSCoverageCleanParams *GetCoverageCleanParams() const
68 : {
69 : GEOSCoverageCleanParams *params =
70 13 : GEOSCoverageCleanParams_create_r(m_poGeosContext);
71 :
72 13 : if (!params)
73 : {
74 0 : CPLError(CE_Failure, CPLE_AppDefined,
75 : "Failed to create coverage clean parameters");
76 0 : return nullptr;
77 : }
78 :
79 13 : if (!GEOSCoverageCleanParams_setSnappingDistance_r(
80 13 : m_poGeosContext, params, m_opts.snappingTolerance))
81 : {
82 0 : CPLError(CE_Failure, CPLE_AppDefined,
83 : "Failed to set snapping tolerance");
84 0 : GEOSCoverageCleanParams_destroy_r(m_poGeosContext, params);
85 0 : return nullptr;
86 : }
87 :
88 13 : if (!GEOSCoverageCleanParams_setGapMaximumWidth_r(
89 13 : m_poGeosContext, params, m_opts.maximumGapWidth))
90 : {
91 0 : CPLError(CE_Failure, CPLE_AppDefined,
92 : "Failed to set maximum gap width");
93 0 : GEOSCoverageCleanParams_destroy_r(m_poGeosContext, params);
94 0 : return nullptr;
95 : }
96 :
97 : int mergeStrategy;
98 13 : if (m_opts.mergeStrategy == "longest-border")
99 : {
100 10 : mergeStrategy = GEOS_MERGE_LONGEST_BORDER;
101 : }
102 3 : else if (m_opts.mergeStrategy == "max-area")
103 : {
104 1 : mergeStrategy = GEOS_MERGE_MAX_AREA;
105 : }
106 2 : else if (m_opts.mergeStrategy == "min-area")
107 : {
108 1 : mergeStrategy = GEOS_MERGE_MIN_AREA;
109 : }
110 1 : else if (m_opts.mergeStrategy == "min-index")
111 : {
112 1 : mergeStrategy = GEOS_MERGE_MIN_INDEX;
113 : }
114 : else
115 : {
116 0 : CPLError(CE_Failure, CPLE_AppDefined,
117 : "Unknown overlap merge strategy: %s",
118 0 : m_opts.mergeStrategy.c_str());
119 0 : GEOSCoverageCleanParams_destroy_r(m_poGeosContext, params);
120 0 : return nullptr;
121 : }
122 :
123 13 : if (!GEOSCoverageCleanParams_setOverlapMergeStrategy_r(
124 13 : m_poGeosContext, params, mergeStrategy))
125 : {
126 0 : CPLError(CE_Failure, CPLE_AppDefined,
127 : "Failed to set overlap merge strategy");
128 0 : GEOSCoverageCleanParams_destroy_r(m_poGeosContext, params);
129 0 : return nullptr;
130 : }
131 :
132 13 : return params;
133 : }
134 :
135 13 : bool Process(OGRLayer &srcLayer, OGRLayer &dstLayer) override
136 : {
137 26 : std::vector<OGRFeatureUniquePtr> features;
138 26 : std::vector<GEOSGeometry *> geoms;
139 :
140 13 : GEOSCoverageCleanParams *params = GetCoverageCleanParams();
141 :
142 : // Copy features from srcLayer into dstLayer, converting
143 : // their geometries to GEOS
144 57 : for (auto &feature : srcLayer)
145 : {
146 47 : const OGRGeometry *fgeom = feature->GetGeometryRef();
147 :
148 : const auto eFGType =
149 47 : fgeom ? wkbFlatten(fgeom->getGeometryType()) : wkbUnknown;
150 47 : if (eFGType != wkbPolygon && eFGType != wkbMultiPolygon &&
151 3 : eFGType != wkbCurvePolygon && eFGType != wkbMultiSurface)
152 : {
153 3 : for (auto &geom : geoms)
154 : {
155 0 : GEOSGeom_destroy_r(m_poGeosContext, geom);
156 : }
157 3 : GEOSCoverageCleanParams_destroy_r(m_poGeosContext, params);
158 3 : CPLError(CE_Failure, CPLE_AppDefined,
159 : "Coverage cleaning can only be performed on "
160 : "polygonal geometries. Feature %" PRId64
161 : " does not have one",
162 3 : static_cast<int64_t>(feature->GetFID()));
163 3 : return false;
164 : }
165 :
166 : GEOSGeometry *geosGeom =
167 44 : fgeom->exportToGEOS(m_poGeosContext, false);
168 44 : if (!geosGeom)
169 : {
170 : // should not happen normally
171 0 : for (auto &geom : geoms)
172 : {
173 0 : GEOSGeom_destroy_r(m_poGeosContext, geom);
174 : }
175 0 : GEOSCoverageCleanParams_destroy_r(m_poGeosContext, params);
176 0 : CPLError(CE_Failure, CPLE_AppDefined,
177 : "Geometry of feature %" PRId64
178 : " failed to convert to GEOS",
179 0 : static_cast<int64_t>(feature->GetFID()));
180 0 : return false;
181 : }
182 44 : geoms.push_back(geosGeom);
183 :
184 44 : feature->SetGeometry(nullptr); // free some memory
185 44 : feature->SetFDefnUnsafe(dstLayer.GetLayerDefn());
186 :
187 44 : features.push_back(std::move(feature));
188 : }
189 :
190 : // Perform coverage cleaning
191 10 : GEOSGeometry *coll = GEOSGeom_createCollection_r(
192 : m_poGeosContext, GEOS_GEOMETRYCOLLECTION, geoms.data(),
193 10 : static_cast<unsigned int>(geoms.size()));
194 :
195 10 : if (coll == nullptr)
196 : {
197 0 : for (auto &geom : geoms)
198 : {
199 0 : GEOSGeom_destroy_r(m_poGeosContext, geom);
200 : }
201 0 : GEOSCoverageCleanParams_destroy_r(m_poGeosContext, params);
202 0 : return false;
203 : }
204 :
205 : GEOSGeometry *geos_result =
206 10 : GEOSCoverageCleanWithParams_r(m_poGeosContext, coll, params);
207 10 : GEOSGeom_destroy_r(m_poGeosContext, coll);
208 10 : GEOSCoverageCleanParams_destroy_r(m_poGeosContext, params);
209 :
210 10 : if (geos_result == nullptr)
211 : {
212 0 : return false;
213 : }
214 :
215 10 : m_papoGeosResults = GEOSGeom_releaseCollection_r(
216 : m_poGeosContext, geos_result, &m_nGeosResultSize);
217 10 : GEOSGeom_destroy_r(m_poGeosContext, geos_result);
218 10 : CPLAssert(features.size() == m_nGeosResultSize);
219 :
220 : // Create features with the modified geometries
221 54 : for (size_t i = 0; i < features.size(); i++)
222 : {
223 44 : GEOSGeometry *dstGeom = m_papoGeosResults[i];
224 :
225 : std::unique_ptr<OGRGeometry> poSimplified(
226 44 : OGRGeometryFactory::createFromGEOS(m_poGeosContext, dstGeom));
227 44 : GEOSGeom_destroy_r(m_poGeosContext, m_papoGeosResults[i]);
228 44 : m_papoGeosResults[i] = nullptr;
229 :
230 44 : if (poSimplified == nullptr)
231 : {
232 0 : CPLError(CE_Failure, CPLE_AppDefined,
233 : "Failed to convert result from GEOS");
234 0 : return false;
235 : }
236 88 : poSimplified->assignSpatialReference(
237 44 : dstLayer.GetLayerDefn()->GetGeomFieldDefn(0)->GetSpatialRef());
238 44 : features[i]->SetGeometry(std::move(poSimplified));
239 :
240 44 : if (dstLayer.CreateFeature(features[i].get()) != CE_None)
241 : {
242 0 : return false;
243 : }
244 44 : features[i].reset();
245 : }
246 :
247 10 : return true;
248 : }
249 :
250 : private:
251 : CPL_DISALLOW_COPY_ASSIGN(GDALVectorCleanCoverageOutputDataset)
252 :
253 : const GDALVectorCleanCoverageAlgorithm::Options &m_opts;
254 : GEOSContextHandle_t m_poGeosContext{nullptr};
255 : GEOSGeometry **m_papoGeosResults{nullptr};
256 : unsigned int m_nGeosResultSize{0};
257 : };
258 :
259 28 : GDALVectorCleanCoverageOutputDataset::~GDALVectorCleanCoverageOutputDataset()
260 : {
261 14 : if (m_poGeosContext != nullptr)
262 : {
263 58 : for (size_t i = 0; i < m_nGeosResultSize; i++)
264 : {
265 44 : GEOSGeom_destroy_r(m_poGeosContext, m_papoGeosResults[i]);
266 : }
267 14 : GEOSFree_r(m_poGeosContext, m_papoGeosResults);
268 14 : finishGEOS_r(m_poGeosContext);
269 : }
270 28 : }
271 :
272 14 : bool GDALVectorCleanCoverageAlgorithm::RunStep(GDALPipelineStepRunContext &)
273 : {
274 14 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
275 : auto poDstDS =
276 28 : std::make_unique<GDALVectorCleanCoverageOutputDataset>(m_opts);
277 :
278 14 : bool bFoundActiveLayer = false;
279 :
280 27 : for (auto &&poSrcLayer : poSrcDS->GetLayers())
281 : {
282 20 : if (m_activeLayer.empty() ||
283 4 : m_activeLayer == poSrcLayer->GetDescription())
284 : {
285 13 : if (!poDstDS->AddProcessedLayer(*poSrcLayer))
286 : {
287 3 : return false;
288 : }
289 10 : bFoundActiveLayer = true;
290 : }
291 : else
292 : {
293 3 : poDstDS->AddPassThroughLayer(*poSrcLayer);
294 : }
295 : }
296 :
297 11 : if (!bFoundActiveLayer)
298 : {
299 1 : ReportError(CE_Failure, CPLE_AppDefined,
300 : "Specified layer '%s' was not found",
301 : m_activeLayer.c_str());
302 1 : return false;
303 : }
304 :
305 10 : m_outputDataset.Set(std::move(poDstDS));
306 :
307 10 : return true;
308 : }
309 :
310 : #else
311 :
312 : bool GDALVectorCleanCoverageAlgorithm::RunStep(GDALPipelineStepRunContext &)
313 : {
314 : ReportError(CE_Failure, CPLE_AppDefined,
315 : "%s requires GDAL to be built against version 3.14 or later of "
316 : "the GEOS library.",
317 : NAME);
318 : return false;
319 : }
320 : #endif // HAVE_GEOS
321 :
322 : GDALVectorCleanCoverageAlgorithmStandalone::
323 : ~GDALVectorCleanCoverageAlgorithmStandalone() = default;
324 :
325 : //! @endcond
|
