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          49 : GDALVectorCheckCoverageAlgorithm::GDALVectorCheckCoverageAlgorithm(
      30          49 :     bool standaloneStep)
      31             :     : GDALVectorPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
      32          49 :                                       standaloneStep)
      33             : {
      34             :     AddArg("include-valid", 0,
      35             :            _("Include valid inputs in output, with empty geometry"),
      36          49 :            &m_includeValid);
      37             : 
      38             :     AddArg("geometry-field", 0, _("Name of geometry field to check"),
      39          49 :            &m_geomField);
      40             : 
      41             :     AddArg("maximum-gap-width", 0, _("Maximum width of a gap to be flagged"),
      42          98 :            &m_maximumGapWidth)
      43          49 :         .SetMinValueIncluded(0);
      44          49 : }
      45             : 
      46             : #if defined HAVE_GEOS &&                                                       \
      47             :     (GEOS_VERSION_MAJOR > 3 ||                                                 \
      48             :      (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 12))
      49             : 
      50          18 : class GDALVectorCheckCoverageOutputDataset final
      51             :     : public GDALGeosNonStreamingAlgorithmDataset
      52             : {
      53             :   public:
      54           9 :     explicit GDALVectorCheckCoverageOutputDataset(double maximumGapWidth,
      55             :                                                   bool includeValid)
      56           9 :         : m_maximumGapWidth(maximumGapWidth), m_includeValid(includeValid)
      57             :     {
      58           9 :     }
      59             : 
      60             :     ~GDALVectorCheckCoverageOutputDataset() override;
      61             : 
      62          32 :     bool PolygonsOnly() const override
      63             :     {
      64          32 :         return true;
      65             :     }
      66             : 
      67          32 :     bool SkipEmpty() const override
      68             :     {
      69          32 :         return !m_includeValid;
      70             :     }
      71             : 
      72           8 :     bool ProcessGeos() override
      73             :     {
      74             :         // Perform coverage checking
      75           8 :         GEOSGeometry *coll = GEOSGeom_createCollection_r(
      76             :             m_poGeosContext, GEOS_GEOMETRYCOLLECTION, m_apoGeosInputs.data(),
      77           8 :             static_cast<unsigned int>(m_apoGeosInputs.size()));
      78             : 
      79           8 :         if (coll == nullptr)
      80             :         {
      81           0 :             return false;
      82             :         }
      83             : 
      84           8 :         m_apoGeosInputs.clear();
      85             : 
      86             :         int geos_result =
      87           8 :             GEOSCoverageIsValid_r(m_poGeosContext, coll, m_maximumGapWidth,
      88             :                                   &m_poGeosResultAsCollection);
      89           8 :         GEOSGeom_destroy_r(m_poGeosContext, coll);
      90             : 
      91           8 :         CPLDebug("CoverageIsValid", "%d", geos_result);
      92             : 
      93           8 :         return geos_result != 2;
      94             :     }
      95             : 
      96             :   private:
      97             :     const double m_maximumGapWidth;
      98             :     const bool m_includeValid;
      99             : };
     100             : 
     101             : GDALVectorCheckCoverageOutputDataset::~GDALVectorCheckCoverageOutputDataset() =
     102             :     default;
     103             : 
     104           9 : bool GDALVectorCheckCoverageAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
     105             : {
     106           9 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     107             :     auto poDstDS = std::make_unique<GDALVectorCheckCoverageOutputDataset>(
     108          18 :         m_maximumGapWidth, m_includeValid);
     109             : 
     110           9 :     const bool bSingleLayerOutput = m_inputLayerNames.empty()
     111          12 :                                         ? poSrcDS->GetLayerCount() == 1
     112           3 :                                         : m_inputLayerNames.size() == 1;
     113             : 
     114          18 :     GDALVectorAlgorithmLayerProgressHelper progressHelper(ctxt);
     115             : 
     116          19 :     for (auto &&poSrcLayer : poSrcDS->GetLayers())
     117             :     {
     118          15 :         if (m_inputLayerNames.empty() ||
     119           0 :             std::find(m_inputLayerNames.begin(), m_inputLayerNames.end(),
     120          15 :                       poSrcLayer->GetDescription()) != m_inputLayerNames.end())
     121             :         {
     122          11 :             const auto poSrcLayerDefn = poSrcLayer->GetLayerDefn();
     123          11 :             if (poSrcLayerDefn->GetGeomFieldCount() == 0)
     124             :             {
     125           2 :                 if (m_inputLayerNames.empty())
     126           1 :                     continue;
     127           1 :                 ReportError(CE_Failure, CPLE_AppDefined,
     128             :                             "Specified layer '%s' has no geometry field",
     129           1 :                             poSrcLayer->GetDescription());
     130           1 :                 return false;
     131             :             }
     132             : 
     133           9 :             progressHelper.AddProcessedLayer(*poSrcLayer);
     134             :         }
     135             :     }
     136             : 
     137          17 :     for ([[maybe_unused]] auto [poSrcLayer, bProcessed, layerProgressFunc,
     138          33 :                                 layerProgressData] : progressHelper)
     139             :     {
     140           9 :         const auto poSrcLayerDefn = poSrcLayer->GetLayerDefn();
     141             :         const int geomFieldIndex =
     142           9 :             m_geomField.empty()
     143           9 :                 ? 0
     144           2 :                 : poSrcLayerDefn->GetGeomFieldIndex(m_geomField.c_str());
     145             : 
     146           9 :         if (geomFieldIndex == -1)
     147             :         {
     148           1 :             ReportError(CE_Failure, CPLE_AppDefined,
     149             :                         "Specified geometry field '%s' does not exist in "
     150             :                         "layer '%s'",
     151           1 :                         m_geomField.c_str(), poSrcLayer->GetDescription());
     152           1 :             return false;
     153             :         }
     154             : 
     155             :         OGRFeatureDefn defn(bSingleLayerOutput
     156             :                                 ? "invalid_edge"
     157          16 :                                 : std::string("invalid_edge_")
     158           4 :                                       .append(poSrcLayer->GetDescription())
     159          16 :                                       .c_str());
     160           8 :         defn.SetGeomType(wkbMultiLineString);
     161          16 :         defn.GetGeomFieldDefn(0)->SetSpatialRef(
     162           8 :             poSrcLayerDefn->GetGeomFieldDefn(geomFieldIndex)->GetSpatialRef());
     163             : 
     164           8 :         if (!poDstDS->AddProcessedLayer(*poSrcLayer, defn, geomFieldIndex,
     165             :                                         layerProgressFunc,
     166             :                                         layerProgressData.get()))
     167             :         {
     168           0 :             return false;
     169             :         }
     170             :     }
     171             : 
     172           7 :     m_outputDataset.Set(std::move(poDstDS));
     173             : 
     174           7 :     return true;
     175             : }
     176             : 
     177             : #else
     178             : 
     179             : bool GDALVectorCheckCoverageAlgorithm::RunStep(GDALPipelineStepRunContext &)
     180             : {
     181             :     ReportError(CE_Failure, CPLE_AppDefined,
     182             :                 "%s requires GDAL to be built against version 3.12 or later of "
     183             :                 "the GEOS library.",
     184             :                 NAME);
     185             :     return false;
     186             : }
     187             : #endif  // HAVE_GEOS
     188             : 
     189             : GDALVectorCheckCoverageAlgorithmStandalone::
     190             :     ~GDALVectorCheckCoverageAlgorithmStandalone() = default;
     191             : 
     192             : //! @endcond