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