Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: "gdal vector simplify-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_simplify_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 29 : GDALVectorSimplifyCoverageAlgorithm::GDALVectorSimplifyCoverageAlgorithm(
31 29 : bool standaloneStep)
32 : : GDALVectorPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
33 29 : standaloneStep)
34 : {
35 29 : AddActiveLayerArg(&m_activeLayer);
36 : AddArg("tolerance", 0, _("Distance tolerance for simplification."),
37 58 : &m_opts.tolerance)
38 29 : .SetRequired()
39 29 : .SetMinValueIncluded(0);
40 : AddArg("preserve-boundary", 0,
41 : _("Whether the exterior boundary should be preserved."),
42 29 : &m_opts.preserveBoundary);
43 29 : }
44 :
45 : #if defined HAVE_GEOS && \
46 : (GEOS_VERSION_MAJOR > 3 || \
47 : (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 12))
48 :
49 : class GDALVectorSimplifyCoverageOutputDataset
50 : : public GDALVectorNonStreamingAlgorithmDataset
51 : {
52 : public:
53 8 : GDALVectorSimplifyCoverageOutputDataset(
54 : const GDALVectorSimplifyCoverageAlgorithm::Options &opts)
55 8 : : m_opts(opts)
56 : {
57 8 : }
58 :
59 : ~GDALVectorSimplifyCoverageOutputDataset() override;
60 :
61 7 : bool Process(OGRLayer &srcLayer, OGRLayer &dstLayer) override
62 : {
63 14 : std::vector<OGRFeatureUniquePtr> features;
64 14 : std::vector<GEOSGeometry *> geoms;
65 7 : m_poGeosContext = OGRGeometry::createGEOSContext();
66 :
67 : // Copy features from srcLayer into dstLayer, converting
68 : // their geometries to GEOS
69 47 : for (auto &feature : srcLayer)
70 : {
71 43 : const OGRGeometry *fgeom = feature->GetGeometryRef();
72 :
73 : // Avoid segfault with non-polygonal inputs on GEOS < 3.12.2
74 : // Later versions produce an error instead
75 : const auto eFGType =
76 43 : fgeom ? wkbFlatten(fgeom->getGeometryType()) : wkbUnknown;
77 43 : if (eFGType != wkbPolygon && eFGType != wkbMultiPolygon &&
78 3 : eFGType != wkbCurvePolygon && eFGType != wkbMultiSurface)
79 : {
80 3 : for (auto &geom : geoms)
81 : {
82 0 : GEOSGeom_destroy_r(m_poGeosContext, geom);
83 : }
84 3 : CPLError(CE_Failure, CPLE_AppDefined,
85 : "Coverage simplification can only be performed on "
86 : "polygonal geometries. Feature %" PRId64
87 : " does not have one",
88 3 : static_cast<int64_t>(feature->GetFID()));
89 3 : return false;
90 : }
91 :
92 : GEOSGeometry *geosGeom =
93 40 : fgeom->exportToGEOS(m_poGeosContext, false);
94 40 : if (!geosGeom)
95 : {
96 : // should not happen normally
97 0 : for (auto &geom : geoms)
98 : {
99 0 : GEOSGeom_destroy_r(m_poGeosContext, geom);
100 : }
101 0 : CPLError(CE_Failure, CPLE_AppDefined,
102 : "Geometry of feature %" PRId64
103 : " failed to convert to GEOS",
104 0 : static_cast<int64_t>(feature->GetFID()));
105 0 : return false;
106 : }
107 40 : geoms.push_back(geosGeom);
108 :
109 40 : feature->SetGeometry(nullptr); // free some memory
110 40 : feature->SetFDefnUnsafe(dstLayer.GetLayerDefn());
111 :
112 40 : features.push_back(std::move(feature));
113 : }
114 :
115 : // Perform coverage simplifciation
116 4 : GEOSGeometry *coll = GEOSGeom_createCollection_r(
117 : m_poGeosContext, GEOS_GEOMETRYCOLLECTION, geoms.data(),
118 4 : static_cast<unsigned int>(geoms.size()));
119 :
120 4 : if (coll == nullptr)
121 : {
122 0 : for (auto &geom : geoms)
123 : {
124 0 : GEOSGeom_destroy_r(m_poGeosContext, geom);
125 : }
126 0 : return false;
127 : }
128 :
129 8 : GEOSGeometry *geos_result = GEOSCoverageSimplifyVW_r(
130 4 : m_poGeosContext, coll, m_opts.tolerance, m_opts.preserveBoundary);
131 4 : GEOSGeom_destroy_r(m_poGeosContext, coll);
132 :
133 4 : if (geos_result == nullptr)
134 : {
135 0 : return false;
136 : }
137 :
138 4 : m_papoGeosResults = GEOSGeom_releaseCollection_r(
139 : m_poGeosContext, geos_result, &m_nGeosResultSize);
140 4 : GEOSGeom_destroy_r(m_poGeosContext, geos_result);
141 4 : CPLAssert(features.size() == m_nGeosResultSize);
142 :
143 : // Create features with the modified geometries
144 44 : for (size_t i = 0; i < features.size(); i++)
145 : {
146 40 : GEOSGeometry *dstGeom = m_papoGeosResults[i];
147 :
148 : std::unique_ptr<OGRGeometry> poSimplified(
149 40 : OGRGeometryFactory::createFromGEOS(m_poGeosContext, dstGeom));
150 40 : GEOSGeom_destroy_r(m_poGeosContext, m_papoGeosResults[i]);
151 40 : m_papoGeosResults[i] = nullptr;
152 :
153 40 : if (poSimplified == nullptr)
154 : {
155 0 : CPLError(CE_Failure, CPLE_AppDefined,
156 : "Failed to convert result from GEOS");
157 0 : return false;
158 : }
159 80 : poSimplified->assignSpatialReference(
160 40 : dstLayer.GetLayerDefn()->GetGeomFieldDefn(0)->GetSpatialRef());
161 40 : features[i]->SetGeometry(std::move(poSimplified));
162 :
163 40 : if (dstLayer.CreateFeature(features[i].get()) != CE_None)
164 : {
165 0 : return false;
166 : }
167 40 : features[i].reset();
168 : }
169 :
170 4 : return true;
171 : }
172 :
173 : private:
174 : CPL_DISALLOW_COPY_ASSIGN(GDALVectorSimplifyCoverageOutputDataset)
175 :
176 : const GDALVectorSimplifyCoverageAlgorithm::Options &m_opts;
177 : GEOSContextHandle_t m_poGeosContext{nullptr};
178 : GEOSGeometry **m_papoGeosResults{nullptr};
179 : unsigned int m_nGeosResultSize{0};
180 : };
181 :
182 16 : GDALVectorSimplifyCoverageOutputDataset::
183 8 : ~GDALVectorSimplifyCoverageOutputDataset()
184 : {
185 8 : if (m_poGeosContext != nullptr)
186 : {
187 47 : for (size_t i = 0; i < m_nGeosResultSize; i++)
188 : {
189 40 : GEOSGeom_destroy_r(m_poGeosContext, m_papoGeosResults[i]);
190 : }
191 7 : GEOSFree_r(m_poGeosContext, m_papoGeosResults);
192 7 : finishGEOS_r(m_poGeosContext);
193 : }
194 16 : }
195 :
196 8 : bool GDALVectorSimplifyCoverageAlgorithm::RunStep(GDALPipelineStepRunContext &)
197 : {
198 8 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
199 : auto poDstDS =
200 16 : std::make_unique<GDALVectorSimplifyCoverageOutputDataset>(m_opts);
201 :
202 8 : bool bFoundActiveLayer = false;
203 :
204 15 : for (auto &&poSrcLayer : poSrcDS->GetLayers())
205 : {
206 14 : if (m_activeLayer.empty() ||
207 4 : m_activeLayer == poSrcLayer->GetDescription())
208 : {
209 7 : if (!poDstDS->AddProcessedLayer(*poSrcLayer))
210 : {
211 3 : return false;
212 : }
213 4 : bFoundActiveLayer = true;
214 : }
215 : else
216 : {
217 3 : poDstDS->AddPassThroughLayer(*poSrcLayer);
218 : }
219 : }
220 :
221 5 : if (!bFoundActiveLayer)
222 : {
223 1 : ReportError(CE_Failure, CPLE_AppDefined,
224 : "Specified layer '%s' was not found",
225 : m_activeLayer.c_str());
226 1 : return false;
227 : }
228 :
229 4 : m_outputDataset.Set(std::move(poDstDS));
230 :
231 4 : return true;
232 : }
233 :
234 : #else
235 :
236 : bool GDALVectorSimplifyCoverageAlgorithm::RunStep(GDALPipelineStepRunContext &)
237 : {
238 : ReportError(CE_Failure, CPLE_AppDefined,
239 : "%s requires GDAL to be built against version 3.12 or later of "
240 : "the GEOS library.",
241 : NAME);
242 : return false;
243 : }
244 : #endif // HAVE_GEOS
245 :
246 : GDALVectorSimplifyCoverageAlgorithmStandalone::
247 : ~GDALVectorSimplifyCoverageAlgorithmStandalone() = default;
248 :
249 : //! @endcond
|
