Skip to content

Commit 0b2c609

Browse files
abellgithubabellgithub
committed
Merge remote-tracking branch 'upstream/master' into inputmasking
2 parents 856082c + 48ecd76 commit 0b2c609

27 files changed

+961
-140
lines changed

‎apps/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -68,6 +68,7 @@ add_library(
6868
gdalalg_vector.cpp
6969
gdalalg_vector_info.cpp
7070
gdalalg_vector_clip.cpp
71+
gdalalg_vector_clean_coverage.cpp
7172
gdalalg_vector_concat.cpp
7273
gdalalg_vector_convert.cpp
7374
gdalalg_vector_edit.cpp

‎apps/gdalalg_vector.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414

1515
#include "gdalalg_vector_info.h"
1616
#include "gdalalg_vector_buffer.h"
17+
#include "gdalalg_vector_clean_coverage.h"
1718
#include "gdalalg_vector_clip.h"
1819
#include "gdalalg_vector_concat.h"
1920
#include "gdalalg_vector_convert.h"
@@ -63,6 +64,7 @@ class GDALVectorAlgorithm final : public GDALAlgorithm
6364

6465
RegisterSubAlgorithm<GDALVectorInfoAlgorithm>();
6566
RegisterSubAlgorithm<GDALVectorBufferAlgorithmStandalone>();
67+
RegisterSubAlgorithm<GDALVectorCleanCoverageAlgorithmStandalone>();
6668
RegisterSubAlgorithm<GDALVectorClipAlgorithmStandalone>();
6769
RegisterSubAlgorithm<GDALVectorConcatAlgorithmStandalone>();
6870
RegisterSubAlgorithm<GDALVectorConvertAlgorithm>();

‎apps/gdalalg_vector_clean_coverage.cpp‎

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

0 commit comments

Comments
 (0)