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