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