Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: "clip" step of "vector pipeline", or "gdal vector clip" standalone
5 : * Author: Even Rouault <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalalg_vector_clip.h"
14 :
15 : #include "gdal_priv.h"
16 : #include "gdal_utils.h"
17 : #include "ogrsf_frmts.h"
18 :
19 : #include <algorithm>
20 :
21 : //! @cond Doxygen_Suppress
22 :
23 : #ifndef _
24 : #define _(x) (x)
25 : #endif
26 :
27 : /************************************************************************/
28 : /* GDALVectorClipAlgorithm::GDALVectorClipAlgorithm() */
29 : /************************************************************************/
30 :
31 74 : GDALVectorClipAlgorithm::GDALVectorClipAlgorithm(bool standaloneStep)
32 : : GDALVectorPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
33 74 : standaloneStep)
34 : {
35 74 : AddActiveLayerArg(&m_activeLayer);
36 74 : AddBBOXArg(&m_bbox, _("Clipping bounding box as xmin,ymin,xmax,ymax"))
37 74 : .SetMutualExclusionGroup("bbox-geometry-like");
38 148 : AddArg("bbox-crs", 0, _("CRS of clipping bounding box"), &m_bboxCrs)
39 148 : .SetIsCRSArg()
40 74 : .AddHiddenAlias("bbox_srs");
41 148 : AddArg("geometry", 0, _("Clipping geometry (WKT or GeoJSON)"), &m_geometry)
42 74 : .SetMutualExclusionGroup("bbox-geometry-like");
43 148 : AddArg("geometry-crs", 0, _("CRS of clipping geometry"), &m_geometryCrs)
44 148 : .SetIsCRSArg()
45 74 : .AddHiddenAlias("geometry_srs");
46 : AddArg("like", 0, _("Dataset to use as a template for bounds"),
47 148 : &m_likeDataset, GDAL_OF_RASTER | GDAL_OF_VECTOR)
48 148 : .SetMetaVar("DATASET")
49 74 : .SetMutualExclusionGroup("bbox-geometry-like");
50 : AddArg("like-sql", 0, ("SELECT statement to run on the 'like' dataset"),
51 148 : &m_likeSQL)
52 148 : .SetMetaVar("SELECT-STATEMENT")
53 74 : .SetMutualExclusionGroup("sql-where");
54 : AddArg("like-layer", 0, ("Name of the layer of the 'like' dataset"),
55 148 : &m_likeLayer)
56 74 : .SetMetaVar("LAYER-NAME");
57 : AddArg("like-where", 0, ("WHERE SQL clause to run on the 'like' dataset"),
58 148 : &m_likeWhere)
59 148 : .SetMetaVar("WHERE-EXPRESSION")
60 74 : .SetMutualExclusionGroup("sql-where");
61 74 : }
62 :
63 : /************************************************************************/
64 : /* GDALVectorClipAlgorithmLayer */
65 : /************************************************************************/
66 :
67 : namespace
68 : {
69 : class GDALVectorClipAlgorithmLayer final : public GDALVectorPipelineOutputLayer
70 : {
71 : public:
72 32 : GDALVectorClipAlgorithmLayer(OGRLayer &oSrcLayer,
73 : std::unique_ptr<OGRGeometry> poClipGeom)
74 32 : : GDALVectorPipelineOutputLayer(oSrcLayer),
75 32 : m_poClipGeom(std::move(poClipGeom)),
76 32 : m_eSrcLayerGeomType(oSrcLayer.GetGeomType()),
77 64 : m_eFlattenSrcLayerGeomType(wkbFlatten(m_eSrcLayerGeomType)),
78 64 : m_bSrcLayerGeomTypeIsCollection(OGR_GT_IsSubClassOf(
79 32 : m_eFlattenSrcLayerGeomType, wkbGeometryCollection)),
80 64 : m_poFeatureDefn(oSrcLayer.GetLayerDefn()->Clone())
81 : {
82 32 : SetDescription(oSrcLayer.GetDescription());
83 32 : SetMetadata(oSrcLayer.GetMetadata());
84 32 : oSrcLayer.SetSpatialFilter(m_poClipGeom.get());
85 32 : m_poFeatureDefn->Reference();
86 32 : }
87 :
88 64 : ~GDALVectorClipAlgorithmLayer() override
89 32 : {
90 32 : m_poFeatureDefn->Release();
91 64 : }
92 :
93 426 : const OGRFeatureDefn *GetLayerDefn() const override
94 : {
95 426 : return m_poFeatureDefn;
96 : }
97 :
98 500 : void TranslateFeature(
99 : std::unique_ptr<OGRFeature> poSrcFeature,
100 : std::vector<std::unique_ptr<OGRFeature>> &apoOutFeatures) override
101 : {
102 0 : std::unique_ptr<OGRGeometry> poIntersection;
103 500 : auto poGeom = poSrcFeature->GetGeometryRef();
104 :
105 500 : if (poGeom == nullptr)
106 0 : return;
107 :
108 500 : poIntersection.reset(poGeom->Intersection(m_poClipGeom.get()));
109 500 : if (!poIntersection)
110 0 : return;
111 1000 : poIntersection->assignSpatialReference(
112 500 : m_poFeatureDefn->GetGeomFieldDefn(0)->GetSpatialRef());
113 :
114 500 : poSrcFeature->SetFDefnUnsafe(m_poFeatureDefn);
115 :
116 500 : const auto eSrcGeomType = wkbFlatten(poGeom->getGeometryType());
117 : const auto eFeatGeomType =
118 500 : wkbFlatten(poIntersection->getGeometryType());
119 500 : if (eFeatGeomType != eSrcGeomType &&
120 8 : m_eFlattenSrcLayerGeomType != wkbUnknown &&
121 8 : m_eFlattenSrcLayerGeomType != eFeatGeomType)
122 : {
123 : // If the intersection is a collection of geometry and the
124 : // layer geometry type is of non-collection type, create
125 : // one feature per element of the collection.
126 10 : if (!m_bSrcLayerGeomTypeIsCollection &&
127 2 : OGR_GT_IsSubClassOf(eFeatGeomType, wkbGeometryCollection))
128 : {
129 : auto poGeomColl = std::unique_ptr<OGRGeometryCollection>(
130 2 : poIntersection.release()->toGeometryCollection());
131 3 : for (const auto *poSubGeom : poGeomColl.get())
132 : {
133 : auto poDstFeature =
134 4 : std::unique_ptr<OGRFeature>(poSrcFeature->Clone());
135 2 : poDstFeature->SetGeometry(poSubGeom);
136 2 : apoOutFeatures.push_back(std::move(poDstFeature));
137 : }
138 : }
139 7 : else if (OGR_GT_GetCollection(eFeatGeomType) ==
140 7 : m_eFlattenSrcLayerGeomType)
141 : {
142 5 : poIntersection.reset(OGRGeometryFactory::forceTo(
143 5 : poIntersection.release(), m_eSrcLayerGeomType));
144 5 : poSrcFeature->SetGeometryDirectly(poIntersection.release());
145 5 : apoOutFeatures.push_back(std::move(poSrcFeature));
146 : }
147 2 : else if (m_eFlattenSrcLayerGeomType == wkbGeometryCollection)
148 : {
149 2 : auto poGeomColl = std::make_unique<OGRGeometryCollection>();
150 1 : poGeomColl->addGeometry(std::move(poIntersection));
151 1 : poSrcFeature->SetGeometryDirectly(poGeomColl.release());
152 1 : apoOutFeatures.push_back(std::move(poSrcFeature));
153 8 : }
154 : // else discard geometries of incompatible type with the
155 : // layer geometry type
156 : }
157 : else
158 : {
159 492 : poSrcFeature->SetGeometryDirectly(poIntersection.release());
160 492 : apoOutFeatures.push_back(std::move(poSrcFeature));
161 : }
162 : }
163 :
164 62 : int TestCapability(const char *pszCap) const override
165 : {
166 62 : if (EQUAL(pszCap, OLCStringsAsUTF8) ||
167 49 : EQUAL(pszCap, OLCCurveGeometries) || EQUAL(pszCap, OLCZGeometries))
168 20 : return m_srcLayer.TestCapability(pszCap);
169 42 : return false;
170 : }
171 :
172 : private:
173 : std::unique_ptr<OGRGeometry> const m_poClipGeom{};
174 : const OGRwkbGeometryType m_eSrcLayerGeomType;
175 : const OGRwkbGeometryType m_eFlattenSrcLayerGeomType;
176 : const bool m_bSrcLayerGeomTypeIsCollection;
177 : OGRFeatureDefn *const m_poFeatureDefn;
178 :
179 : CPL_DISALLOW_COPY_ASSIGN(GDALVectorClipAlgorithmLayer)
180 : };
181 :
182 : } // namespace
183 :
184 : /************************************************************************/
185 : /* GDALVectorClipAlgorithm::RunStep() */
186 : /************************************************************************/
187 :
188 36 : bool GDALVectorClipAlgorithm::RunStep(GDALPipelineStepRunContext &)
189 : {
190 36 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
191 36 : CPLAssert(poSrcDS);
192 :
193 36 : CPLAssert(m_outputDataset.GetName().empty());
194 36 : CPLAssert(!m_outputDataset.GetDatasetRef());
195 :
196 36 : const int nLayerCount = poSrcDS->GetLayerCount();
197 36 : bool bSrcLayerHasSRS = false;
198 56 : for (int i = 0; i < nLayerCount; ++i)
199 : {
200 35 : auto poSrcLayer = poSrcDS->GetLayer(i);
201 35 : if (poSrcLayer &&
202 37 : (m_activeLayer.empty() ||
203 72 : m_activeLayer == poSrcLayer->GetDescription()) &&
204 34 : poSrcLayer->GetSpatialRef())
205 : {
206 15 : bSrcLayerHasSRS = true;
207 15 : break;
208 : }
209 : }
210 :
211 72 : auto [poClipGeom, errMsg] = GetClipGeometry();
212 36 : if (!poClipGeom)
213 : {
214 8 : ReportError(CE_Failure, CPLE_AppDefined, "%s", errMsg.c_str());
215 8 : return false;
216 : }
217 :
218 28 : auto poLikeDS = m_likeDataset.GetDatasetRef();
219 29 : if (bSrcLayerHasSRS && !poClipGeom->getSpatialReference() && poLikeDS &&
220 1 : poLikeDS->GetLayerCount() == 0)
221 : {
222 1 : ReportError(CE_Warning, CPLE_AppDefined,
223 : "Dataset '%s' has no CRS. Assuming its CRS is the "
224 : "same as the input vector.",
225 1 : poLikeDS->GetDescription());
226 : }
227 :
228 28 : auto outDS = std::make_unique<GDALVectorPipelineOutputDataset>(*poSrcDS);
229 :
230 28 : bool ret = true;
231 61 : for (int i = 0; ret && i < nLayerCount; ++i)
232 : {
233 33 : auto poSrcLayer = poSrcDS->GetLayer(i);
234 33 : ret = (poSrcLayer != nullptr);
235 33 : if (ret)
236 : {
237 35 : if (m_activeLayer.empty() ||
238 2 : m_activeLayer == poSrcLayer->GetDescription())
239 : {
240 : auto poClipGeomForLayer =
241 64 : std::unique_ptr<OGRGeometry>(poClipGeom->clone());
242 40 : if (poClipGeomForLayer->getSpatialReference() &&
243 8 : poSrcLayer->GetSpatialRef())
244 : {
245 14 : ret = poClipGeomForLayer->transformTo(
246 7 : poSrcLayer->GetSpatialRef()) == OGRERR_NONE;
247 : }
248 32 : if (ret)
249 : {
250 64 : outDS->AddLayer(
251 : *poSrcLayer,
252 64 : std::make_unique<GDALVectorClipAlgorithmLayer>(
253 32 : *poSrcLayer, std::move(poClipGeomForLayer)));
254 : }
255 : }
256 : else
257 : {
258 2 : outDS->AddLayer(
259 : *poSrcLayer,
260 2 : std::make_unique<GDALVectorPipelinePassthroughLayer>(
261 : *poSrcLayer));
262 : }
263 : }
264 : }
265 :
266 28 : if (ret)
267 28 : m_outputDataset.Set(std::move(outDS));
268 :
269 28 : return ret;
270 : }
271 :
272 : GDALVectorClipAlgorithmStandalone::~GDALVectorClipAlgorithmStandalone() =
273 : default;
274 :
275 : //! @endcond
|