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