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             : 
      15             : #include "cpl_conv.h"
      16             : #include "gdal_priv.h"
      17             : #include "gdal_alg.h"
      18             : #include "ogrsf_frmts.h"
      19             : 
      20             : #include <cmath>
      21             : #include <optional>
      22             : 
      23             : //! @cond Doxygen_Suppress
      24             : 
      25             : #ifndef _
      26             : #define _(x) (x)
      27             : #endif
      28             : 
      29          18 : GDALRasterAsFeaturesAlgorithm::GDALRasterAsFeaturesAlgorithm(
      30          18 :     bool standaloneStep)
      31             :     : GDALPipelineStepAlgorithm(
      32             :           NAME, DESCRIPTION, HELP_URL,
      33           0 :           ConstructorOptions()
      34          18 :               .SetStandaloneStep(standaloneStep)
      35          36 :               .SetOutputFormatCreateCapability(GDAL_DCAP_CREATE)),
      36             :       m_bands{}, m_geomTypeName("none"), m_skipNoData(false),
      37          54 :       m_includeXY(false), m_includeRowCol(false)
      38             : {
      39          18 :     m_outputLayerName = "pixels";
      40             : 
      41             :     // TODO: This block is copied from gdalalg_raster_polygonize. Avoid duplication?
      42          18 :     if (standaloneStep)
      43             :     {
      44          18 :         AddOutputFormatArg(&m_format).AddMetadataItem(
      45          54 :             GAAMDI_REQUIRED_CAPABILITIES, {GDAL_DCAP_VECTOR, GDAL_DCAP_CREATE});
      46          18 :         AddOpenOptionsArg(&m_openOptions);
      47          18 :         AddInputFormatsArg(&m_inputFormats)
      48          36 :             .AddMetadataItem(GAAMDI_REQUIRED_CAPABILITIES, {GDAL_DCAP_RASTER});
      49             : 
      50          18 :         AddInputDatasetArg(&m_inputDataset, GDAL_OF_RASTER);
      51          18 :         AddOutputDatasetArg(&m_outputDataset, GDAL_OF_VECTOR)
      52          18 :             .SetDatasetInputFlags(GADV_NAME | GADV_OBJECT);
      53          18 :         AddCreationOptionsArg(&m_creationOptions);
      54          18 :         AddLayerCreationOptionsArg(&m_layerCreationOptions);
      55          18 :         AddOverwriteArg(&m_overwrite);
      56          18 :         AddUpdateArg(&m_update);
      57          18 :         AddOverwriteLayerArg(&m_overwriteLayer);
      58          18 :         AddAppendLayerArg(&m_appendLayer);
      59          18 :         AddLayerNameArg(&m_outputLayerName)
      60          36 :             .AddAlias("nln")
      61          18 :             .SetDefault(m_outputLayerName);
      62             :     }
      63             : 
      64          18 :     AddBandArg(&m_bands);
      65          36 :     AddArg("geometry-type", 0, _("Geometry type"), &m_geomTypeName)
      66          18 :         .SetChoices("none", "point", "polygon")
      67          18 :         .SetDefault("none");
      68             :     AddArg("skip-nodata", 0, _("Omit NoData pixels from the result"),
      69          18 :            &m_skipNoData);
      70             :     AddArg("include-xy", 0, _("Include fields for cell center coordinates"),
      71          18 :            &m_includeXY);
      72             :     AddArg("include-row-col", 0, _("Include columns for row and column"),
      73          18 :            &m_includeRowCol);
      74          18 : }
      75             : 
      76             : GDALRasterAsFeaturesAlgorithm::~GDALRasterAsFeaturesAlgorithm() = default;
      77             : 
      78          14 : bool GDALRasterAsFeaturesAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
      79             :                                             void *pProgressData)
      80             : {
      81          14 :     GDALPipelineStepRunContext stepCtxt;
      82          14 :     stepCtxt.m_pfnProgress = pfnProgress;
      83          14 :     stepCtxt.m_pProgressData = pProgressData;
      84          14 :     return RunPreStepPipelineValidations() && RunStep(stepCtxt);
      85             : }
      86             : 
      87             : GDALRasterAsFeaturesAlgorithmStandalone::
      88             :     ~GDALRasterAsFeaturesAlgorithmStandalone() = default;
      89             : 
      90             : struct RasterAsFeaturesOptions
      91             : {
      92             :     OGRwkbGeometryType geomType{wkbNone};
      93             :     bool includeXY{false};
      94             :     bool includeRowCol{false};
      95             :     bool skipNoData{false};
      96             :     std::vector<int> bands{};
      97             : };
      98             : 
      99             : class GDALRasterToVectorPipelineOutputDataset final : public GDALDataset
     100             : {
     101             : 
     102             :   public:
     103           0 :     int GetLayerCount() const override
     104             :     {
     105           0 :         return static_cast<int>(m_layers.size());
     106             :     }
     107             : 
     108           0 :     const OGRLayer *GetLayer(int idx) const override
     109             :     {
     110           0 :         return m_layers[idx].get();
     111             :     }
     112             : 
     113             :     int TestCapability(const char *) const override;
     114             : 
     115             :     void AddLayer(std::unique_ptr<OGRLayer> layer)
     116             :     {
     117             :         m_layers.emplace_back(std::move(layer));
     118             :     }
     119             : 
     120             :   private:
     121             :     std::vector<std::unique_ptr<OGRLayer>> m_layers{};
     122             : };
     123             : 
     124           0 : int GDALRasterToVectorPipelineOutputDataset::TestCapability(const char *) const
     125             : {
     126           0 :     return 0;
     127             : }
     128             : 
     129             : class GDALRasterAsFeaturesLayer final : public OGRLayer
     130             : {
     131             :   public:
     132             :     static constexpr const char *ROW_FIELD = "ROW";
     133             :     static constexpr const char *COL_FIELD = "COL";
     134             :     static constexpr const char *X_FIELD = "CENTER_X";
     135             :     static constexpr const char *Y_FIELD = "CENTER_Y";
     136             : 
     137          14 :     GDALRasterAsFeaturesLayer(GDALDataset &ds, RasterAsFeaturesOptions options)
     138          14 :         : m_ds(ds), m_bufType(GDT_Float64),
     139             :           m_it(GDALRasterBand::WindowIterator(
     140          28 :               m_ds.GetRasterXSize(), m_ds.GetRasterYSize(),
     141          28 :               m_ds.GetRasterXSize(), m_ds.GetRasterYSize(), 0, 0)),
     142             :           m_end(GDALRasterBand::WindowIterator(
     143          28 :               m_ds.GetRasterXSize(), m_ds.GetRasterYSize(),
     144          28 :               m_ds.GetRasterXSize(), m_ds.GetRasterYSize(), 1, 0)),
     145          14 :           m_includeXY(options.includeXY),
     146          14 :           m_includeRowCol(options.includeRowCol),
     147          14 :           m_excludeNoDataPixels(options.skipNoData)
     148             :     {
     149             :         // TODO: Handle Int64, UInt64
     150          14 :         m_ds.GetGeoTransform(m_gt);
     151             : 
     152          14 :         int nBands = m_ds.GetRasterCount();
     153          14 :         m_bands.resize(nBands);
     154          31 :         for (int i = 1; i <= nBands; i++)
     155             :         {
     156          17 :             m_bands[i - 1] = i;
     157             :         }
     158             : 
     159             :         // TODO: Handle per-band NoData values
     160          14 :         if (nBands > 0)
     161             :         {
     162             :             int hasNoData;
     163          13 :             double noData = m_ds.GetRasterBand(1)->GetNoDataValue(&hasNoData);
     164          13 :             if (hasNoData)
     165             :             {
     166           1 :                 m_noData = noData;
     167             :             }
     168             :         }
     169             : 
     170          14 :         m_defn = new OGRFeatureDefn();
     171          14 :         m_defn->GetGeomFieldDefn(0)->SetType(options.geomType);
     172          14 :         m_defn->GetGeomFieldDefn(0)->SetSpatialRef(ds.GetSpatialRef());
     173          14 :         m_defn->Reference();
     174             : 
     175          14 :         if (m_includeXY)
     176             :         {
     177           2 :             auto xField = std::make_unique<OGRFieldDefn>(X_FIELD, OFTReal);
     178           2 :             auto yField = std::make_unique<OGRFieldDefn>(Y_FIELD, OFTReal);
     179           1 :             m_defn->AddFieldDefn(std::move(xField));
     180           1 :             m_defn->AddFieldDefn(std::move(yField));
     181             :         }
     182          14 :         if (m_includeRowCol)
     183             :         {
     184             :             auto rowField =
     185           8 :                 std::make_unique<OGRFieldDefn>(ROW_FIELD, OFTInteger);
     186             :             auto colField =
     187           8 :                 std::make_unique<OGRFieldDefn>(COL_FIELD, OFTInteger);
     188           4 :             m_defn->AddFieldDefn(std::move(rowField));
     189           4 :             m_defn->AddFieldDefn(std::move(colField));
     190             :         }
     191          31 :         for (int band : m_bands)
     192             :         {
     193          34 :             CPLString fieldName = CPLSPrintf("BAND_%d", band);
     194             :             auto bandField =
     195          17 :                 std::make_unique<OGRFieldDefn>(fieldName.c_str(), OFTReal);
     196          17 :             m_defn->AddFieldDefn(std::move(bandField));
     197          17 :             m_bandFields.push_back(m_defn->GetFieldIndex(fieldName));
     198             :         }
     199             : 
     200          14 :         GDALRasterAsFeaturesLayer::ResetReading();
     201          14 :     }
     202             : 
     203             :     ~GDALRasterAsFeaturesLayer() override;
     204             : 
     205          28 :     void ResetReading() override
     206             :     {
     207          28 :         if (m_ds.GetRasterCount() > 0)
     208             :         {
     209          26 :             GDALRasterBand *poFirstBand = m_ds.GetRasterBand(1);
     210          26 :             CPLAssert(poFirstBand);  // appease clang scan-build
     211          26 :             m_it = poFirstBand->IterateWindows().begin();
     212          26 :             m_end = poFirstBand->IterateWindows().end();
     213             :         }
     214          28 :     }
     215             : 
     216           0 :     int TestCapability(const char *) const override
     217             :     {
     218           0 :         return 0;
     219             :     }
     220             : 
     221          28 :     OGRFeatureDefn *GetLayerDefn() const override
     222             :     {
     223          28 :         return m_defn;
     224             :     }
     225             : 
     226        8598 :     OGRFeature *GetNextFeature() override
     227             :     {
     228        8598 :         if (m_row >= m_window.nYSize && !NextWindow())
     229             :         {
     230          14 :             return nullptr;
     231             :         }
     232             : 
     233       17168 :         std::unique_ptr<OGRFeature> feature;
     234             : 
     235        8604 :         while (m_row < m_window.nYSize)
     236             :         {
     237       17208 :             const double *pSrcVal = static_cast<double *>(m_buf) +
     238        8604 :                                     (m_bands.size() * m_row * m_window.nXSize +
     239        8604 :                                      m_col * m_bands.size());
     240             : 
     241             :             const bool emitFeature =
     242        8604 :                 !m_excludeNoDataPixels || !IsNoData(*pSrcVal);
     243             : 
     244        8604 :             if (emitFeature)
     245             :             {
     246        8584 :                 feature.reset(OGRFeature::CreateFeature(m_defn));
     247             : 
     248       27164 :                 for (int fieldPos : m_bandFields)
     249             :                 {
     250       18580 :                     feature->SetField(fieldPos, *pSrcVal);
     251       18580 :                     pSrcVal++;
     252             :                 }
     253             : 
     254        8584 :                 const double line = m_window.nYOff + m_row;
     255        8584 :                 const double pixel = m_window.nXOff + m_col;
     256             : 
     257        8584 :                 if (m_includeRowCol)
     258             :                 {
     259         404 :                     feature->SetField(ROW_FIELD, static_cast<GIntBig>(line));
     260         404 :                     feature->SetField(COL_FIELD, static_cast<GIntBig>(pixel));
     261             :                 }
     262        8584 :                 if (m_includeXY)
     263             :                 {
     264             :                     double x, y;
     265         400 :                     m_gt.Apply(pixel + 0.5, line + 0.5, &x, &y);
     266         400 :                     feature->SetField(X_FIELD, x);
     267         400 :                     feature->SetField(Y_FIELD, y);
     268             :                 }
     269             : 
     270           0 :                 std::unique_ptr<OGRGeometry> geom;
     271        8584 :                 const auto geomType = m_defn->GetGeomFieldDefn(0)->GetType();
     272        8584 :                 if (geomType == wkbPoint)
     273             :                 {
     274             :                     double x, y;
     275        2900 :                     m_gt.Apply(pixel + 0.5, line + 0.5, &x, &y);
     276             : 
     277        2900 :                     geom = std::make_unique<OGRPoint>(x, y);
     278             :                 }
     279        5684 :                 else if (geomType == wkbPolygon)
     280             :                 {
     281             :                     double x, y;
     282             : 
     283         800 :                     auto lr = std::make_unique<OGRLinearRing>();
     284             : 
     285         400 :                     m_gt.Apply(pixel, line, &x, &y);
     286         400 :                     lr->addPoint(x, y);
     287         400 :                     m_gt.Apply(pixel, line + 1, &x, &y);
     288         400 :                     lr->addPoint(x, y);
     289         400 :                     m_gt.Apply(pixel + 1, line + 1, &x, &y);
     290         400 :                     lr->addPoint(x, y);
     291         400 :                     m_gt.Apply(pixel + 1, line, &x, &y);
     292         400 :                     lr->addPoint(x, y);
     293         400 :                     m_gt.Apply(pixel, line, &x, &y);
     294         400 :                     lr->addPoint(x, y);
     295             : 
     296         800 :                     auto poly = std::make_unique<OGRPolygon>();
     297         400 :                     poly->addRing(std::move(lr));
     298         400 :                     geom = std::move(poly);
     299             :                 }
     300             : 
     301        8584 :                 feature->SetGeometry(std::move(geom));
     302             :             }
     303             : 
     304        8604 :             m_col += 1;
     305        8604 :             if (m_col >= m_window.nXSize)
     306             :             {
     307         282 :                 m_col = 0;
     308         282 :                 m_row++;
     309             :             }
     310             : 
     311        8604 :             if (feature)
     312             :             {
     313        8584 :                 return feature.release();
     314             :             }
     315             :         }
     316             : 
     317           0 :         return nullptr;
     318             :     }
     319             : 
     320             :     CPL_DISALLOW_COPY_ASSIGN(GDALRasterAsFeaturesLayer)
     321             : 
     322             :   private:
     323         400 :     bool IsNoData(double x) const
     324             :     {
     325         400 :         if (!m_noData.has_value())
     326             :         {
     327           0 :             return false;
     328             :         }
     329             : 
     330         780 :         return m_noData.value() == x ||
     331         780 :                (std::isnan(m_noData.value()) && std::isnan(x));
     332             :     }
     333             : 
     334          26 :     bool NextWindow()
     335             :     {
     336          26 :         int nBandCount = static_cast<int>(m_bands.size());
     337             : 
     338          26 :         if (m_it == m_end)
     339             :         {
     340          13 :             return false;
     341             :         }
     342             : 
     343          13 :         if (m_ds.GetRasterXSize() == 0 || m_ds.GetRasterYSize() == 0)
     344             :         {
     345           1 :             return false;
     346             :         }
     347             : 
     348          12 :         m_window = *m_it;
     349          12 :         ++m_it;
     350             : 
     351          12 :         if (m_buf == nullptr && !m_bands.empty())
     352             :         {
     353          11 :             m_buf =
     354          11 :                 VSI_MALLOC3_VERBOSE(m_window.nXSize, m_window.nYSize,
     355             :                                     static_cast<size_t>(nBandCount) *
     356             :                                         GDALGetDataTypeSizeBytes(m_bufType));
     357          11 :             if (m_buf == nullptr)
     358             :             {
     359           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     360             :                          "Failed to allocate buffer");
     361           0 :                 return false;
     362             :             }
     363             :         }
     364             : 
     365          27 :         for (size_t i = 0; i < m_bands.size(); i++)
     366             :         {
     367             :             auto eErr =
     368          15 :                 m_ds.GetRasterBand(m_bands[i])
     369          15 :                     ->RasterIO(GF_Read, m_window.nXOff, m_window.nYOff,
     370             :                                m_window.nXSize, m_window.nYSize,
     371          15 :                                static_cast<GByte *>(m_buf) +
     372          15 :                                    i * GDALGetDataTypeSizeBytes(m_bufType),
     373             :                                m_window.nXSize, m_window.nYSize, m_bufType,
     374          15 :                                static_cast<GSpacing>(nBandCount) *
     375          15 :                                    GDALGetDataTypeSizeBytes(m_bufType),
     376          30 :                                static_cast<GSpacing>(nBandCount) *
     377          15 :                                    GDALGetDataTypeSizeBytes(m_bufType) *
     378          15 :                                    m_window.nXSize,
     379             :                                nullptr);
     380             : 
     381          15 :             if (eErr != CE_None)
     382             :             {
     383           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     384             :                          "Failed to read raster data");
     385           0 :                 return false;
     386             :             }
     387             :         }
     388             : 
     389          12 :         m_row = 0;
     390          12 :         m_col = 0;
     391             : 
     392          12 :         return true;
     393             :     }
     394             : 
     395             :     GDALDataset &m_ds;
     396             :     void *m_buf{nullptr};
     397             :     GDALDataType m_bufType;
     398             :     GDALGeoTransform m_gt{};
     399             :     std::optional<double> m_noData{std::nullopt};
     400             : 
     401             :     std::vector<int> m_bands{};
     402             :     std::vector<int> m_bandFields{};
     403             : 
     404             :     GDALRasterBand::WindowIterator m_it;
     405             :     GDALRasterBand::WindowIterator m_end;
     406             :     GDALRasterWindow m_window{};
     407             : 
     408             :     int m_row{0};
     409             :     int m_col{0};
     410             : 
     411             :     OGRFeatureDefn *m_defn{nullptr};
     412             :     bool m_includeXY;
     413             :     bool m_includeRowCol;
     414             :     bool m_excludeNoDataPixels;
     415             : };
     416             : 
     417          28 : GDALRasterAsFeaturesLayer::~GDALRasterAsFeaturesLayer()
     418             : {
     419          14 :     m_defn->Release();
     420          14 :     if (m_buf != nullptr)
     421             :     {
     422          11 :         CPLFree(m_buf);
     423             :     }
     424          28 : }
     425             : 
     426          14 : bool GDALRasterAsFeaturesAlgorithm::RunStep(GDALPipelineStepRunContext &)
     427             : {
     428          14 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     429             : 
     430          14 :     GDALDataset *poDstDS = m_outputDataset.GetDatasetRef();
     431          28 :     std::string outputFilename = m_outputDataset.GetName();
     432             : 
     433          14 :     std::unique_ptr<GDALDataset> poRetDS;
     434          14 :     if (!poDstDS)
     435             :     {
     436          14 :         if (m_standaloneStep && m_format.empty())
     437             :         {
     438             :             const auto aosFormats =
     439             :                 CPLStringList(GDALGetOutputDriversForDatasetName(
     440          14 :                     m_outputDataset.GetName().c_str(), GDAL_OF_VECTOR,
     441             :                     /* bSingleMatch = */ true,
     442          14 :                     /* bWarn = */ true));
     443          14 :             if (aosFormats.size() != 1)
     444             :             {
     445           0 :                 ReportError(CE_Failure, CPLE_AppDefined,
     446             :                             "Cannot guess driver for %s",
     447           0 :                             m_outputDataset.GetName().c_str());
     448           0 :                 return false;
     449             :             }
     450          14 :             m_format = aosFormats[0];
     451             : 
     452             :             auto poDriver =
     453          14 :                 GetGDALDriverManager()->GetDriverByName(m_format.c_str());
     454             : 
     455          14 :             poRetDS.reset(
     456             :                 poDriver->Create(outputFilename.c_str(), 0, 0, 0, GDT_Unknown,
     457          28 :                                  CPLStringList(m_creationOptions).List()));
     458             :         }
     459             :         else
     460             :         {
     461             :             poRetDS =
     462           0 :                 std::make_unique<GDALRasterToVectorPipelineOutputDataset>();
     463             :         }
     464             : 
     465          14 :         if (!poRetDS)
     466           0 :             return false;
     467             : 
     468          14 :         poDstDS = poRetDS.get();
     469             :     }
     470             : 
     471          28 :     RasterAsFeaturesOptions options;
     472          26 :     options.geomType = m_geomTypeName == "point"     ? wkbPoint
     473          12 :                        : m_geomTypeName == "polygon" ? wkbPolygon
     474             :                                                      : wkbNone;
     475          14 :     options.includeRowCol = m_includeRowCol;
     476          14 :     options.includeXY = m_includeXY;
     477          14 :     options.skipNoData = m_skipNoData;
     478             : 
     479          14 :     if (!m_bands.empty())
     480             :     {
     481           1 :         options.bands = std::move(m_bands);
     482             :     }
     483             : 
     484             :     auto poLayer =
     485          14 :         std::make_unique<GDALRasterAsFeaturesLayer>(*poSrcDS, options);
     486             : 
     487          14 :     poDstDS->CopyLayer(poLayer.get(), m_outputLayerName.c_str());
     488             : 
     489          14 :     if (poRetDS)
     490             :     {
     491          14 :         m_outputDataset.Set(std::move(poRetDS));
     492             :     }
     493             : 
     494          14 :     return true;
     495             : }
     496             : 
     497             : //! @endcond