Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: "gdal vector sort" 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_sort.h"
14 :
15 : #include "cpl_error.h"
16 : #include "gdal_alg.h"
17 : #include "gdal_priv.h"
18 : #include "gdalalg_vector_geom.h"
19 : #include "ogr_geometry.h"
20 :
21 : #include "ogr_geos.h"
22 :
23 : #include <cinttypes>
24 : #include <limits>
25 :
26 : #ifndef _
27 : #define _(x) (x)
28 : #endif
29 :
30 : //! @cond Doxygen_Suppress
31 :
32 49 : GDALVectorSortAlgorithm::GDALVectorSortAlgorithm(bool standaloneStep)
33 : : GDALVectorPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
34 49 : standaloneStep)
35 : {
36 : AddArg("geometry-field", 0, _("Name of geometry field to use in sort"),
37 49 : &m_geomField);
38 98 : AddArg("method", 0, _("Geometry sorting algorithm"), &m_sortMethod)
39 49 : .SetChoices("hilbert", "strtree")
40 49 : .SetDefault(m_sortMethod);
41 49 : }
42 :
43 : namespace
44 : {
45 :
46 8 : bool CreateDstFeatures(
47 : const std::vector<std::unique_ptr<OGRFeature>> &srcFeatures,
48 : const std::vector<size_t> &sortedIndices, OGRLayer &dstLayer)
49 : {
50 428 : for (size_t iSrcFeature : sortedIndices)
51 : {
52 420 : OGRFeature *poSrcFeature = srcFeatures[iSrcFeature].get();
53 420 : poSrcFeature->SetFDefnUnsafe(dstLayer.GetLayerDefn());
54 420 : poSrcFeature->SetFID(OGRNullFID);
55 :
56 420 : if (dstLayer.CreateFeature(poSrcFeature) != OGRERR_NONE)
57 : {
58 0 : return false;
59 : }
60 : }
61 :
62 8 : return true;
63 : }
64 :
65 : class GDALVectorHilbertSortDataset
66 : : public GDALVectorNonStreamingAlgorithmDataset
67 : {
68 : public:
69 : using GDALVectorNonStreamingAlgorithmDataset::
70 : GDALVectorNonStreamingAlgorithmDataset;
71 :
72 5 : bool Process(OGRLayer &srcLayer, OGRLayer &dstLayer,
73 : int geomFieldIndex) override
74 : {
75 10 : std::vector<std::unique_ptr<OGRFeature>> features;
76 :
77 225 : for (auto &feature : srcLayer)
78 : {
79 220 : features.emplace_back(feature.release());
80 : }
81 5 : OGREnvelope oLayerExtent;
82 :
83 10 : std::vector<OGREnvelope> envelopes(features.size());
84 225 : for (size_t i = 0; i < features.size(); i++)
85 : {
86 220 : const OGRFeature *poFeature = features[i].get();
87 : const OGRGeometry *poGeom =
88 220 : poFeature->GetGeomFieldRef(geomFieldIndex);
89 :
90 220 : if (poGeom != nullptr && !poGeom->IsEmpty())
91 : {
92 212 : poGeom->getEnvelope(&envelopes[i]);
93 212 : oLayerExtent.Merge(envelopes[i]);
94 : }
95 : }
96 :
97 : std::vector<std::pair<std::size_t, std::uint32_t>> hilbertCodes(
98 10 : features.size());
99 225 : for (std::size_t i = 0; i < features.size(); i++)
100 : {
101 220 : hilbertCodes[i].first = i;
102 :
103 220 : if (envelopes[i].IsInit())
104 : {
105 : double dfX, dfY;
106 212 : envelopes[i].Center(dfX, dfY);
107 212 : hilbertCodes[i].second =
108 212 : GDALHilbertCode(&oLayerExtent, dfX, dfY);
109 : }
110 : else
111 : {
112 8 : hilbertCodes[i].second =
113 8 : std::numeric_limits<std::uint32_t>::max();
114 : }
115 : }
116 :
117 5 : std::sort(hilbertCodes.begin(), hilbertCodes.end(),
118 1455 : [](const auto &a, const auto &b)
119 1455 : { return a.second < b.second; });
120 :
121 10 : std::vector<size_t> sortedIndices;
122 225 : for (const auto &sItem : hilbertCodes)
123 : {
124 220 : sortedIndices.push_back(sItem.first);
125 : }
126 :
127 10 : return CreateDstFeatures(features, sortedIndices, dstLayer);
128 : }
129 :
130 : private:
131 : CPL_DISALLOW_COPY_ASSIGN(GDALVectorHilbertSortDataset)
132 : };
133 :
134 : #ifdef HAVE_GEOS
135 : class GDALVectorSTRTreeSortDataset
136 : : public GDALVectorNonStreamingAlgorithmDataset
137 : {
138 : public:
139 : using GDALVectorNonStreamingAlgorithmDataset::
140 : GDALVectorNonStreamingAlgorithmDataset;
141 :
142 8 : ~GDALVectorSTRTreeSortDataset() override
143 4 : {
144 4 : if (m_geosContext)
145 : {
146 3 : finishGEOS_r(m_geosContext);
147 : }
148 8 : }
149 :
150 3 : bool Process(OGRLayer &srcLayer, OGRLayer &dstLayer,
151 : int geomFieldIndex) override
152 : {
153 6 : std::vector<std::unique_ptr<OGRFeature>> features;
154 6 : std::vector<size_t> sortedIndices;
155 :
156 203 : for (auto &feature : srcLayer)
157 : {
158 200 : features.emplace_back(feature.release());
159 : }
160 : // TODO: variant of this fn returning unique_ptr
161 3 : m_geosContext = OGRGeometry::createGEOSContext();
162 :
163 3 : auto TreeDeleter = [this](GEOSSTRtree *tree)
164 3 : { GEOSSTRtree_destroy_r(m_geosContext, tree); };
165 :
166 : std::unique_ptr<GEOSSTRtree, decltype(TreeDeleter)> poTree(
167 6 : GEOSSTRtree_create_r(m_geosContext, 10), TreeDeleter);
168 :
169 3 : OGREnvelope oGeomExtent;
170 6 : std::vector<size_t> nullIndices;
171 203 : for (size_t i = 0; i < features.size(); i++)
172 : {
173 200 : const OGRFeature *poFeature = features[i].get();
174 : const OGRGeometry *poGeom =
175 200 : poFeature->GetGeomFieldRef(geomFieldIndex);
176 :
177 200 : if (poGeom == nullptr || poGeom->IsEmpty())
178 : {
179 8 : nullIndices.push_back(i);
180 8 : continue;
181 : }
182 :
183 192 : poGeom->getEnvelope(&oGeomExtent);
184 192 : GEOSGeometry *poEnv = CreateGEOSEnvelope(oGeomExtent);
185 192 : if (poEnv == nullptr)
186 : {
187 0 : return false;
188 : }
189 :
190 192 : GEOSSTRtree_insert_r(m_geosContext, poTree.get(), poEnv,
191 : reinterpret_cast<void *>(i));
192 192 : GEOSGeom_destroy_r(m_geosContext, poEnv);
193 : }
194 :
195 : #if GEOS_VERSION_MAJOR > 3 || \
196 : (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 12)
197 3 : GEOSSTRtree_build_r(m_geosContext, poTree.get());
198 : #else
199 : if (!features.empty())
200 : {
201 : // Perform a dummy query to force tree construction.
202 : GEOSGeometry *poEnv = CreateGEOSEnvelope(oGeomExtent);
203 : GEOSSTRtree_query_r(
204 : m_geosContext, poTree.get(), poEnv, [](void *, void *) {},
205 : nullptr);
206 : GEOSGeom_destroy_r(m_geosContext, poEnv);
207 : }
208 : #endif
209 :
210 3 : GEOSSTRtree_iterate_r(
211 : m_geosContext, poTree.get(),
212 192 : [](void *item, void *userData)
213 : {
214 192 : static_cast<std::vector<size_t> *>(userData)->push_back(
215 192 : reinterpret_cast<size_t>(item));
216 192 : },
217 : &sortedIndices);
218 :
219 11 : for (size_t nullInd : nullIndices)
220 : {
221 8 : sortedIndices.push_back(nullInd);
222 : }
223 :
224 3 : return CreateDstFeatures(features, sortedIndices, dstLayer);
225 : }
226 :
227 : private:
228 : CPL_DISALLOW_COPY_ASSIGN(GDALVectorSTRTreeSortDataset)
229 :
230 : // FIXME: Duplicated from alg/zonal.cpp.
231 : // Put into OGRGeometryFactory?
232 192 : GEOSGeometry *CreateGEOSEnvelope(const OGREnvelope &oEnv) const
233 : {
234 192 : GEOSCoordSequence *seq = GEOSCoordSeq_create_r(m_geosContext, 2, 2);
235 192 : if (seq == nullptr)
236 : {
237 0 : return nullptr;
238 : }
239 192 : GEOSCoordSeq_setXY_r(m_geosContext, seq, 0, oEnv.MinX, oEnv.MinY);
240 192 : GEOSCoordSeq_setXY_r(m_geosContext, seq, 1, oEnv.MaxX, oEnv.MaxY);
241 192 : return GEOSGeom_createLineString_r(m_geosContext, seq);
242 : }
243 :
244 : GEOSContextHandle_t m_geosContext{nullptr};
245 : };
246 : #endif
247 :
248 : } // namespace
249 :
250 12 : bool GDALVectorSortAlgorithm::RunStep(GDALPipelineStepRunContext &)
251 : {
252 12 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
253 12 : std::unique_ptr<GDALVectorNonStreamingAlgorithmDataset> poDstDS;
254 :
255 12 : if (m_sortMethod == "hilbert")
256 : {
257 8 : poDstDS = std::make_unique<GDALVectorHilbertSortDataset>();
258 : }
259 : else
260 : {
261 : // Already checked for invalid method at arg parsing stage.
262 4 : CPLAssert(m_sortMethod == "strtree");
263 : #ifdef HAVE_GEOS
264 4 : poDstDS = std::make_unique<GDALVectorSTRTreeSortDataset>();
265 : #else
266 : CPLError(
267 : CE_Failure, CPLE_AppDefined,
268 : "--method strtree requires a GDAL build against the GEOS library.");
269 : return false;
270 : #endif
271 : }
272 :
273 22 : for (auto &&poSrcLayer : poSrcDS->GetLayers())
274 : {
275 12 : if (m_inputLayerNames.empty() ||
276 0 : std::find(m_inputLayerNames.begin(), m_inputLayerNames.end(),
277 12 : poSrcLayer->GetDescription()) != m_inputLayerNames.end())
278 : {
279 12 : const auto poSrcLayerDefn = poSrcLayer->GetLayerDefn();
280 12 : if (poSrcLayerDefn->GetGeomFieldCount() > 0)
281 : {
282 : const int geomFieldIndex =
283 10 : m_geomField.empty() ? 0
284 3 : : poSrcLayerDefn->GetGeomFieldIndex(
285 3 : m_geomField.c_str());
286 :
287 10 : if (geomFieldIndex == -1)
288 : {
289 2 : ReportError(
290 : CE_Failure, CPLE_AppDefined,
291 : "Specified geometry field '%s' does not exist in "
292 : "layer '%s'",
293 2 : m_geomField.c_str(), poSrcLayer->GetDescription());
294 2 : return false;
295 : }
296 :
297 16 : if (!poDstDS->AddProcessedLayer(*poSrcLayer,
298 8 : *poSrcLayer->GetLayerDefn(),
299 : geomFieldIndex))
300 : {
301 0 : return false;
302 : }
303 : }
304 2 : else if (m_inputLayerNames.empty())
305 : {
306 2 : poDstDS->AddPassThroughLayer(*poSrcLayer);
307 : }
308 : }
309 : }
310 :
311 10 : m_outputDataset.Set(std::move(poDstDS));
312 :
313 10 : return true;
314 : }
315 :
316 : GDALVectorSortAlgorithmStandalone::~GDALVectorSortAlgorithmStandalone() =
317 : default;
318 :
319 : //! @endcond
|
