Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: "as-features" step of "gdal pipeline"
5 : * Author: Daniel Baston
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2025, ISciences, LLC
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalalg_raster_as_features.h"
14 : #include "gdalalg_vector_pipeline.h"
15 :
16 : #include "cpl_conv.h"
17 : #include "gdal_priv.h"
18 : #include "gdal_alg.h"
19 : #include "ogrsf_frmts.h"
20 :
21 : #include <cmath>
22 : #include <limits>
23 : #include <optional>
24 :
25 : //! @cond Doxygen_Suppress
26 :
27 : #ifndef _
28 : #define _(x) (x)
29 : #endif
30 :
31 18 : GDALRasterAsFeaturesAlgorithm::GDALRasterAsFeaturesAlgorithm(
32 18 : bool standaloneStep)
33 : : GDALPipelineStepAlgorithm(
34 : NAME, DESCRIPTION, HELP_URL,
35 0 : ConstructorOptions()
36 18 : .SetStandaloneStep(standaloneStep)
37 18 : .SetAddUpsertArgument(false)
38 18 : .SetAddSkipErrorsArgument(false)
39 36 : .SetOutputFormatCreateCapability(GDAL_DCAP_CREATE))
40 : {
41 18 : m_outputLayerName = "pixels";
42 :
43 18 : if (standaloneStep)
44 : {
45 18 : AddRasterInputArgs(false, false);
46 18 : AddVectorOutputArgs(false, false);
47 : }
48 :
49 18 : AddBandArg(&m_bands);
50 36 : AddArg("geometry-type", 0, _("Geometry type"), &m_geomTypeName)
51 18 : .SetChoices("none", "point", "polygon")
52 18 : .SetDefault(m_geomTypeName);
53 : AddArg("skip-nodata", 0, _("Omit NoData pixels from the result"),
54 18 : &m_skipNoData);
55 : AddArg("include-xy", 0, _("Include fields for cell center coordinates"),
56 18 : &m_includeXY);
57 : AddArg("include-row-col", 0, _("Include columns for row and column"),
58 18 : &m_includeRowCol);
59 18 : }
60 :
61 : GDALRasterAsFeaturesAlgorithm::~GDALRasterAsFeaturesAlgorithm() = default;
62 :
63 : GDALRasterAsFeaturesAlgorithmStandalone::
64 : ~GDALRasterAsFeaturesAlgorithmStandalone() = default;
65 :
66 : struct RasterAsFeaturesOptions
67 : {
68 : OGRwkbGeometryType geomType{wkbNone};
69 : bool includeXY{false};
70 : bool includeRowCol{false};
71 : bool skipNoData{false};
72 : std::vector<int> bands{};
73 : };
74 :
75 : class GDALRasterAsFeaturesLayer final
76 : : public OGRLayer,
77 : public OGRGetNextFeatureThroughRaw<GDALRasterAsFeaturesLayer>
78 : {
79 : public:
80 : static constexpr const char *ROW_FIELD = "ROW";
81 : static constexpr const char *COL_FIELD = "COL";
82 : static constexpr const char *X_FIELD = "CENTER_X";
83 : static constexpr const char *Y_FIELD = "CENTER_Y";
84 :
85 8598 : DEFINE_GET_NEXT_FEATURE_THROUGH_RAW(GDALRasterAsFeaturesLayer)
86 :
87 14 : GDALRasterAsFeaturesLayer(GDALDataset &ds, RasterAsFeaturesOptions options)
88 14 : : m_ds(ds), m_it(GDALRasterBand::WindowIterator(
89 28 : m_ds.GetRasterXSize(), m_ds.GetRasterYSize(),
90 28 : m_ds.GetRasterXSize(), m_ds.GetRasterYSize(), 0, 0)),
91 : m_end(GDALRasterBand::WindowIterator(
92 28 : m_ds.GetRasterXSize(), m_ds.GetRasterYSize(),
93 28 : m_ds.GetRasterXSize(), m_ds.GetRasterYSize(), 1, 0)),
94 14 : m_includeXY(options.includeXY),
95 14 : m_includeRowCol(options.includeRowCol),
96 14 : m_excludeNoDataPixels(options.skipNoData)
97 : {
98 : // TODO: Handle Int64, UInt64
99 14 : m_ds.GetGeoTransform(m_gt);
100 :
101 14 : int nBands = m_ds.GetRasterCount();
102 14 : m_bands.resize(nBands);
103 31 : for (int i = 1; i <= nBands; i++)
104 : {
105 17 : m_bands[i - 1] = i;
106 : }
107 :
108 : // TODO: Handle per-band NoData values
109 14 : if (nBands > 0)
110 : {
111 : int hasNoData;
112 13 : double noData = m_ds.GetRasterBand(1)->GetNoDataValue(&hasNoData);
113 13 : if (hasNoData)
114 : {
115 1 : m_noData = noData;
116 : }
117 : }
118 :
119 14 : m_defn = new OGRFeatureDefn();
120 14 : if (options.geomType == wkbNone)
121 : {
122 11 : m_defn->SetGeomType(wkbNone);
123 : }
124 : else
125 : {
126 3 : m_defn->GetGeomFieldDefn(0)->SetType(options.geomType);
127 3 : m_defn->GetGeomFieldDefn(0)->SetSpatialRef(ds.GetSpatialRef());
128 : }
129 14 : m_defn->Reference();
130 :
131 14 : if (m_includeXY)
132 : {
133 2 : auto xField = std::make_unique<OGRFieldDefn>(X_FIELD, OFTReal);
134 2 : auto yField = std::make_unique<OGRFieldDefn>(Y_FIELD, OFTReal);
135 1 : m_defn->AddFieldDefn(std::move(xField));
136 1 : m_defn->AddFieldDefn(std::move(yField));
137 : }
138 14 : if (m_includeRowCol)
139 : {
140 : auto rowField =
141 8 : std::make_unique<OGRFieldDefn>(ROW_FIELD, OFTInteger);
142 : auto colField =
143 8 : std::make_unique<OGRFieldDefn>(COL_FIELD, OFTInteger);
144 4 : m_defn->AddFieldDefn(std::move(rowField));
145 4 : m_defn->AddFieldDefn(std::move(colField));
146 : }
147 31 : for (int band : m_bands)
148 : {
149 34 : CPLString fieldName = CPLSPrintf("BAND_%d", band);
150 : auto bandField =
151 17 : std::make_unique<OGRFieldDefn>(fieldName.c_str(), OFTReal);
152 17 : m_defn->AddFieldDefn(std::move(bandField));
153 17 : m_bandFields.push_back(m_defn->GetFieldIndex(fieldName));
154 : }
155 :
156 14 : GDALRasterAsFeaturesLayer::ResetReading();
157 14 : }
158 :
159 : ~GDALRasterAsFeaturesLayer() override;
160 :
161 28 : void ResetReading() override
162 : {
163 28 : if (m_ds.GetRasterCount() > 0)
164 : {
165 26 : GDALRasterBand *poFirstBand = m_ds.GetRasterBand(1);
166 26 : CPLAssert(poFirstBand); // appease clang scan-build
167 26 : m_it = poFirstBand->IterateWindows().begin();
168 26 : m_end = poFirstBand->IterateWindows().end();
169 : }
170 28 : }
171 :
172 15 : int TestCapability(const char *pszCap) const override
173 : {
174 16 : return EQUAL(pszCap, OLCFastFeatureCount) &&
175 16 : m_poFilterGeom == nullptr && m_poAttrQuery == nullptr &&
176 16 : !m_excludeNoDataPixels;
177 : }
178 :
179 1 : GIntBig GetFeatureCount(int bForce) override
180 : {
181 1 : if (m_poFilterGeom == nullptr && m_poAttrQuery == nullptr &&
182 1 : !m_excludeNoDataPixels)
183 : {
184 1 : return static_cast<GIntBig>(m_ds.GetRasterXSize()) *
185 1 : m_ds.GetRasterYSize();
186 : }
187 0 : return OGRLayer::GetFeatureCount(bForce);
188 : }
189 :
190 74 : OGRFeatureDefn *GetLayerDefn() const override
191 : {
192 74 : return m_defn;
193 : }
194 :
195 8598 : OGRFeature *GetNextRawFeature()
196 : {
197 8598 : if (m_row >= m_window.nYSize && !NextWindow())
198 : {
199 14 : return nullptr;
200 : }
201 :
202 17168 : std::unique_ptr<OGRFeature> feature;
203 :
204 8604 : while (m_row < m_window.nYSize)
205 : {
206 8604 : const double *pSrcVal = reinterpret_cast<double *>(m_buf.data()) +
207 8604 : (m_bands.size() * m_row * m_window.nXSize +
208 8604 : m_col * m_bands.size());
209 :
210 : const bool emitFeature =
211 8604 : !m_excludeNoDataPixels || !IsNoData(*pSrcVal);
212 :
213 8604 : if (emitFeature)
214 : {
215 8584 : feature.reset(OGRFeature::CreateFeature(m_defn));
216 :
217 27164 : for (int fieldPos : m_bandFields)
218 : {
219 18580 : feature->SetField(fieldPos, *pSrcVal);
220 18580 : pSrcVal++;
221 : }
222 :
223 8584 : const double line = m_window.nYOff + m_row;
224 8584 : const double pixel = m_window.nXOff + m_col;
225 :
226 8584 : if (m_includeRowCol)
227 : {
228 404 : feature->SetField(ROW_FIELD, static_cast<GIntBig>(line));
229 404 : feature->SetField(COL_FIELD, static_cast<GIntBig>(pixel));
230 : }
231 8584 : if (m_includeXY)
232 : {
233 : double x, y;
234 400 : m_gt.Apply(pixel + 0.5, line + 0.5, &x, &y);
235 400 : feature->SetField(X_FIELD, x);
236 400 : feature->SetField(Y_FIELD, y);
237 : }
238 :
239 0 : std::unique_ptr<OGRGeometry> geom;
240 8584 : const auto geomType = m_defn->GetGeomType();
241 8584 : if (geomType == wkbPoint)
242 : {
243 : double x, y;
244 2900 : m_gt.Apply(pixel + 0.5, line + 0.5, &x, &y);
245 :
246 2900 : geom = std::make_unique<OGRPoint>(x, y);
247 5800 : geom->assignSpatialReference(
248 2900 : m_defn->GetGeomFieldDefn(0)->GetSpatialRef());
249 : }
250 5684 : else if (geomType == wkbPolygon)
251 : {
252 : double x, y;
253 :
254 800 : auto lr = std::make_unique<OGRLinearRing>();
255 :
256 400 : m_gt.Apply(pixel, line, &x, &y);
257 400 : lr->addPoint(x, y);
258 400 : m_gt.Apply(pixel, line + 1, &x, &y);
259 400 : lr->addPoint(x, y);
260 400 : m_gt.Apply(pixel + 1, line + 1, &x, &y);
261 400 : lr->addPoint(x, y);
262 400 : m_gt.Apply(pixel + 1, line, &x, &y);
263 400 : lr->addPoint(x, y);
264 400 : m_gt.Apply(pixel, line, &x, &y);
265 400 : lr->addPoint(x, y);
266 :
267 800 : auto poly = std::make_unique<OGRPolygon>();
268 400 : poly->addRing(std::move(lr));
269 400 : geom = std::move(poly);
270 800 : geom->assignSpatialReference(
271 400 : m_defn->GetGeomFieldDefn(0)->GetSpatialRef());
272 : }
273 :
274 8584 : feature->SetGeometry(std::move(geom));
275 : }
276 :
277 8604 : m_col += 1;
278 8604 : if (m_col >= m_window.nXSize)
279 : {
280 282 : m_col = 0;
281 282 : m_row++;
282 : }
283 :
284 8604 : if (feature)
285 : {
286 8584 : return feature.release();
287 : }
288 : }
289 :
290 0 : return nullptr;
291 : }
292 :
293 : CPL_DISALLOW_COPY_ASSIGN(GDALRasterAsFeaturesLayer)
294 :
295 : private:
296 400 : bool IsNoData(double x) const
297 : {
298 400 : if (!m_noData.has_value())
299 : {
300 0 : return false;
301 : }
302 :
303 780 : return m_noData.value() == x ||
304 780 : (std::isnan(m_noData.value()) && std::isnan(x));
305 : }
306 :
307 26 : bool NextWindow()
308 : {
309 26 : int nBandCount = static_cast<int>(m_bands.size());
310 :
311 26 : if (m_it == m_end)
312 : {
313 13 : return false;
314 : }
315 :
316 13 : if (m_ds.GetRasterXSize() == 0 || m_ds.GetRasterYSize() == 0)
317 : {
318 1 : return false;
319 : }
320 :
321 12 : m_window = *m_it;
322 12 : ++m_it;
323 :
324 12 : if (!m_bands.empty())
325 : {
326 11 : const auto nBufTypeSize = GDALGetDataTypeSizeBytes(m_bufType);
327 : if constexpr (sizeof(int) < sizeof(size_t))
328 : {
329 22 : if (m_window.nYSize > 0 &&
330 11 : static_cast<size_t>(m_window.nXSize) >
331 11 : std::numeric_limits<size_t>::max() / m_window.nYSize)
332 : {
333 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
334 : "Failed to allocate buffer");
335 0 : return false;
336 : }
337 : }
338 11 : const size_t nPixelCount =
339 11 : static_cast<size_t>(m_window.nXSize) * m_window.nYSize;
340 22 : if (static_cast<size_t>(nBandCount) * nBufTypeSize >
341 11 : std::numeric_limits<size_t>::max() / nPixelCount)
342 : {
343 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
344 : "Failed to allocate buffer");
345 0 : return false;
346 : }
347 11 : const size_t nBufSize = nPixelCount * nBandCount * nBufTypeSize;
348 11 : if (m_buf.size() < nBufSize)
349 : {
350 : try
351 : {
352 11 : m_buf.resize(nBufSize);
353 : }
354 0 : catch (const std::exception &)
355 : {
356 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
357 : "Failed to allocate buffer");
358 0 : return false;
359 : }
360 : }
361 :
362 11 : const auto nPixelSpace =
363 11 : static_cast<GSpacing>(nBandCount) * nBufTypeSize;
364 11 : const auto eErr = m_ds.RasterIO(
365 : GF_Read, m_window.nXOff, m_window.nYOff, m_window.nXSize,
366 11 : m_window.nYSize, m_buf.data(), m_window.nXSize, m_window.nYSize,
367 11 : m_bufType, static_cast<int>(m_bands.size()), m_bands.data(),
368 11 : nPixelSpace, nPixelSpace * m_window.nXSize, nBufTypeSize,
369 : nullptr);
370 :
371 11 : if (eErr != CE_None)
372 : {
373 0 : CPLError(CE_Failure, CPLE_AppDefined,
374 : "Failed to read raster data");
375 0 : return false;
376 : }
377 : }
378 :
379 12 : m_row = 0;
380 12 : m_col = 0;
381 :
382 12 : return true;
383 : }
384 :
385 : GDALDataset &m_ds;
386 : std::vector<GByte> m_buf{};
387 : const GDALDataType m_bufType = GDT_Float64;
388 : GDALGeoTransform m_gt{};
389 : std::optional<double> m_noData{std::nullopt};
390 :
391 : std::vector<int> m_bands{};
392 : std::vector<int> m_bandFields{};
393 :
394 : GDALRasterBand::WindowIterator m_it;
395 : GDALRasterBand::WindowIterator m_end;
396 : GDALRasterWindow m_window{};
397 :
398 : int m_row{0};
399 : int m_col{0};
400 :
401 : OGRFeatureDefn *m_defn{nullptr};
402 : bool m_includeXY;
403 : bool m_includeRowCol;
404 : bool m_excludeNoDataPixels;
405 : };
406 :
407 28 : GDALRasterAsFeaturesLayer::~GDALRasterAsFeaturesLayer()
408 : {
409 14 : m_defn->Release();
410 28 : }
411 :
412 14 : bool GDALRasterAsFeaturesAlgorithm::RunStep(GDALPipelineStepRunContext &)
413 : {
414 14 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
415 :
416 28 : RasterAsFeaturesOptions options;
417 26 : options.geomType = m_geomTypeName == "point" ? wkbPoint
418 12 : : m_geomTypeName == "polygon" ? wkbPolygon
419 : : wkbNone;
420 14 : options.includeRowCol = m_includeRowCol;
421 14 : options.includeXY = m_includeXY;
422 14 : options.skipNoData = m_skipNoData;
423 :
424 14 : if (!m_bands.empty())
425 : {
426 1 : options.bands = std::move(m_bands);
427 : }
428 :
429 : auto poLayer =
430 28 : std::make_unique<GDALRasterAsFeaturesLayer>(*poSrcDS, options);
431 14 : auto poRetDS = std::make_unique<GDALVectorOutputDataset>();
432 14 : poRetDS->AddLayer(std::move(poLayer));
433 :
434 14 : m_outputDataset.Set(std::move(poRetDS));
435 :
436 28 : return true;
437 : }
438 :
439 : //! @endcond
|