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