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