Line data    Source code

       1             : /******************************************************************************
       2             : *
       3             :  * Project:  GDAL
       4             :  * Purpose:  "gdal vector check-coverage" subcommand
       5             :  * Author:   Daniel Baston
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2025, ISciences LLC
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "gdalalg_vector_check_coverage.h"
      14             : 
      15             : #include "cpl_error.h"
      16             : #include "gdal_priv.h"
      17             : #include "gdalalg_vector_geom.h"
      18             : #include "ogr_geometry.h"
      19             : #include "ogr_geos.h"
      20             : 
      21             : #include <cinttypes>
      22             : 
      23             : #ifndef _
      24             : #define _(x) (x)
      25             : #endif
      26             : 
      27             : //! @cond Doxygen_Suppress
      28             : 
      29          89 : GDALVectorCheckCoverageAlgorithm::GDALVectorCheckCoverageAlgorithm(
      30          89 :     bool standaloneStep)
      31             :     : GDALVectorPipelineStepAlgorithm(
      32             :           NAME, DESCRIPTION, HELP_URL,
      33           0 :           ConstructorOptions()
      34          89 :               .SetStandaloneStep(standaloneStep)
      35         178 :               .SetNoCreateEmptyLayersArgument(standaloneStep))
      36             : {
      37             :     AddArg("include-valid", 0,
      38             :            _("Include valid inputs in output, with empty geometry"),
      39          89 :            &m_includeValid);
      40             : 
      41             :     AddArg("geometry-field", 0, _("Name of geometry field to check"),
      42          89 :            &m_geomField);
      43             : 
      44             :     AddArg("maximum-gap-width", 0, _("Maximum width of a gap to be flagged"),
      45         178 :            &m_maximumGapWidth)
      46          89 :         .SetMinValueIncluded(0);
      47          89 : }
      48             : 
      49             : #if defined HAVE_GEOS &&                                                       \
      50             :     (GEOS_VERSION_MAJOR > 3 ||                                                 \
      51             :      (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 12))
      52             : 
      53             : class GDALVectorCheckCoverageOutputLayer final
      54             :     : public GDALGeosNonStreamingAlgorithmLayer
      55             : {
      56             :   public:
      57          14 :     explicit GDALVectorCheckCoverageOutputLayer(OGRLayer &srcLayer,
      58             :                                                 int geomFieldIndex,
      59             :                                                 const std::string &name,
      60             :                                                 double maximumGapWidth,
      61             :                                                 bool includeValid)
      62          14 :         : GDALGeosNonStreamingAlgorithmLayer(srcLayer, geomFieldIndex),
      63          14 :           m_defn(OGRFeatureDefn::CreateFeatureDefn(name.c_str())),
      64          14 :           m_maximumGapWidth(maximumGapWidth), m_includeValid(includeValid)
      65             :     {
      66          14 :         m_defn->Reference();
      67             : 
      68          14 :         const OGRFeatureDefn *poSrcLayerDefn = srcLayer.GetLayerDefn();
      69          14 :         m_defn->SetGeomType(wkbMultiLineString);
      70          28 :         m_defn->GetGeomFieldDefn(0)->SetSpatialRef(
      71          14 :             poSrcLayerDefn->GetGeomFieldDefn(geomFieldIndex)->GetSpatialRef());
      72          14 :     }
      73             : 
      74             :     ~GDALVectorCheckCoverageOutputLayer() override;
      75             : 
      76         297 :     const OGRFeatureDefn *GetLayerDefn() const override
      77             :     {
      78         297 :         return m_defn;
      79             :     }
      80             : 
      81          12 :     int TestCapability(const char *) const override
      82             :     {
      83          12 :         return false;
      84             :     }
      85             : 
      86          84 :     bool PolygonsOnly() const override
      87             :     {
      88          84 :         return true;
      89             :     }
      90             : 
      91          67 :     bool SkipEmpty() const override
      92             :     {
      93          67 :         return !m_includeValid;
      94             :     }
      95             : 
      96          14 :     bool ProcessGeos() override
      97             :     {
      98             :         // Perform coverage checking
      99          14 :         GEOSGeometry *coll = GEOSGeom_createCollection_r(
     100             :             m_poGeosContext, GEOS_GEOMETRYCOLLECTION, m_apoGeosInputs.data(),
     101          14 :             static_cast<unsigned int>(m_apoGeosInputs.size()));
     102             : 
     103          14 :         if (coll == nullptr)
     104             :         {
     105           0 :             return false;
     106             :         }
     107             : 
     108          14 :         m_apoGeosInputs.clear();
     109             : 
     110             :         int geos_result =
     111          14 :             GEOSCoverageIsValid_r(m_poGeosContext, coll, m_maximumGapWidth,
     112             :                                   &m_poGeosResultAsCollection);
     113          14 :         GEOSGeom_destroy_r(m_poGeosContext, coll);
     114             : 
     115          14 :         CPLDebug("CoverageIsValid", "%d", geos_result);
     116             : 
     117          14 :         return geos_result != 2;
     118             :     }
     119             : 
     120             :   private:
     121             :     OGRFeatureDefn *m_defn;
     122             :     const double m_maximumGapWidth;
     123             :     const bool m_includeValid;
     124             : 
     125             :     CPL_DISALLOW_COPY_ASSIGN(GDALVectorCheckCoverageOutputLayer)
     126             : };
     127             : 
     128          28 : GDALVectorCheckCoverageOutputLayer::~GDALVectorCheckCoverageOutputLayer()
     129             : {
     130          14 :     if (m_defn != nullptr)
     131             :     {
     132          14 :         m_defn->Release();
     133             :     }
     134          28 : }
     135             : 
     136          12 : bool GDALVectorCheckCoverageAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
     137             : {
     138          12 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     139             :     auto poDstDS =
     140          24 :         std::make_unique<GDALVectorNonStreamingAlgorithmDataset>(*poSrcDS);
     141             : 
     142          12 :     const bool bSingleLayerOutput = m_inputLayerNames.empty()
     143          15 :                                         ? poSrcDS->GetLayerCount() == 1
     144           3 :                                         : m_inputLayerNames.size() == 1;
     145             : 
     146          24 :     GDALVectorAlgorithmLayerProgressHelper progressHelper(ctxt);
     147             : 
     148          28 :     for (auto &&poSrcLayer : poSrcDS->GetLayers())
     149             :     {
     150          21 :         if (m_inputLayerNames.empty() ||
     151           0 :             std::find(m_inputLayerNames.begin(), m_inputLayerNames.end(),
     152          21 :                       poSrcLayer->GetDescription()) != m_inputLayerNames.end())
     153             :         {
     154          17 :             const auto poSrcLayerDefn = poSrcLayer->GetLayerDefn();
     155          17 :             if (poSrcLayerDefn->GetGeomFieldCount() == 0)
     156             :             {
     157           2 :                 if (m_inputLayerNames.empty())
     158           1 :                     continue;
     159           1 :                 ReportError(CE_Failure, CPLE_AppDefined,
     160             :                             "Specified layer '%s' has no geometry field",
     161           1 :                             poSrcLayer->GetDescription());
     162           1 :                 return false;
     163             :             }
     164             : 
     165          15 :             progressHelper.AddProcessedLayer(*poSrcLayer);
     166             :         }
     167             :     }
     168             : 
     169          29 :     for ([[maybe_unused]] auto [poSrcLayer, bProcessed, layerProgressFunc,
     170          54 :                                 layerProgressData] : progressHelper)
     171             :     {
     172          15 :         const auto poSrcLayerDefn = poSrcLayer->GetLayerDefn();
     173             :         const int geomFieldIndex =
     174          15 :             m_geomField.empty()
     175          15 :                 ? 0
     176           2 :                 : poSrcLayerDefn->GetGeomFieldIndex(m_geomField.c_str());
     177             : 
     178          15 :         if (geomFieldIndex == -1)
     179             :         {
     180           1 :             ReportError(CE_Failure, CPLE_AppDefined,
     181             :                         "Specified geometry field '%s' does not exist in "
     182             :                         "layer '%s'",
     183           1 :                         m_geomField.c_str(), poSrcLayer->GetDescription());
     184           1 :             return false;
     185             :         }
     186             : 
     187             :         std::string layerName = bSingleLayerOutput
     188             :                                     ? "invalid_edge"
     189          34 :                                     : std::string("invalid_edge_")
     190          38 :                                           .append(poSrcLayer->GetDescription());
     191             : 
     192             :         auto poLayer = std::make_unique<GDALVectorCheckCoverageOutputLayer>(
     193          14 :             *poSrcLayer, geomFieldIndex, layerName, m_maximumGapWidth,
     194          14 :             m_includeValid);
     195             : 
     196          14 :         if (!poDstDS->AddProcessedLayer(std::move(poLayer), layerProgressFunc,
     197             :                                         layerProgressData.get()))
     198             :         {
     199           0 :             return false;
     200             :         }
     201             :     }
     202             : 
     203          10 :     m_outputDataset.Set(std::move(poDstDS));
     204             : 
     205          10 :     return true;
     206             : }
     207             : 
     208             : #else
     209             : 
     210             : bool GDALVectorCheckCoverageAlgorithm::RunStep(GDALPipelineStepRunContext &)
     211             : {
     212             :     ReportError(CE_Failure, CPLE_AppDefined,
     213             :                 "%s requires GDAL to be built against version 3.12 or later of "
     214             :                 "the GEOS library.",
     215             :                 NAME);
     216             :     return false;
     217             : }
     218             : #endif  // HAVE_GEOS
     219             : 
     220             : GDALVectorCheckCoverageAlgorithmStandalone::
     221             :     ~GDALVectorCheckCoverageAlgorithmStandalone() = default;
     222             : 
     223             : //! @endcond