Line data    Source code

       1             : /******************************************************************************
       2             : *
       3             :  * Project:  GDAL
       4             :  * Purpose:  "gdal vector check-geometry" 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_geometry.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         122 : GDALVectorCheckGeometryAlgorithm::GDALVectorCheckGeometryAlgorithm(
      30         122 :     bool standaloneStep)
      31             :     : GDALVectorPipelineStepAlgorithm(
      32             :           NAME, DESCRIPTION, HELP_URL,
      33           0 :           ConstructorOptions()
      34         122 :               .SetStandaloneStep(standaloneStep)
      35         244 :               .SetNoCreateEmptyLayersArgument(standaloneStep))
      36             : {
      37             :     AddArg("include-field", 0,
      38             :            _("Fields from input layer to include in output (special values: "
      39             :              "ALL and NONE)"),
      40         244 :            &m_includeFields)
      41         122 :         .SetDefault("NONE");
      42             : 
      43             :     AddArg("include-valid", 0,
      44             :            _("Include valid inputs in output, with empty geometry"),
      45         122 :            &m_includeValid);
      46             : 
      47             :     AddArg("geometry-field", 0, _("Name of geometry field to check"),
      48         122 :            &m_geomField);
      49         122 : }
      50             : 
      51             : #ifdef HAVE_GEOS
      52             : 
      53             : class GDALInvalidLocationLayer final : public GDALVectorPipelineOutputLayer
      54             : {
      55             :   private:
      56             :     static constexpr const char *ERROR_DESCRIPTION_FIELD = "error";
      57             : 
      58             :   public:
      59          41 :     GDALInvalidLocationLayer(OGRLayer &layer,
      60             :                              const std::vector<int> &srcFieldIndices,
      61             :                              bool bSingleLayerOutput, int srcGeomField,
      62             :                              bool skipValid)
      63          41 :         : GDALVectorPipelineOutputLayer(layer),
      64             :           m_defn(OGRFeatureDefnRefCountedPtr::makeInstance(
      65             :               bSingleLayerOutput ? "error_location"
      66          49 :                                  : std::string("error_location_")
      67           8 :                                        .append(layer.GetDescription())
      68           8 :                                        .c_str())),
      69          41 :           m_geosContext(OGRGeometry::createGEOSContext()),
      70          57 :           m_srcGeomField(srcGeomField), m_skipValid(skipValid)
      71             :     {
      72          41 :         m_defn->SetGeomType(wkbMultiPoint);
      73             : 
      74          41 :         if (!srcFieldIndices.empty())
      75             :         {
      76           4 :             const OGRFeatureDefn &srcDefn = *layer.GetLayerDefn();
      77           4 :             m_srcFieldMap.resize(srcDefn.GetFieldCount(), -1);
      78           4 :             int iDstField = 0;
      79          13 :             for (int iSrcField : srcFieldIndices)
      80             :             {
      81           9 :                 m_defn->AddFieldDefn(srcDefn.GetFieldDefn(iSrcField));
      82           9 :                 m_srcFieldMap[iSrcField] = iDstField++;
      83             :             }
      84             :         }
      85             : 
      86             :         auto poDescriptionFieldDefn =
      87          82 :             std::make_unique<OGRFieldDefn>(ERROR_DESCRIPTION_FIELD, OFTString);
      88          41 :         m_defn->AddFieldDefn(std::move(poDescriptionFieldDefn));
      89             : 
      90          82 :         m_defn->GetGeomFieldDefn(0)->SetSpatialRef(
      91          41 :             m_srcLayer.GetLayerDefn()
      92          41 :                 ->GetGeomFieldDefn(m_srcGeomField)
      93          41 :                 ->GetSpatialRef());
      94          41 :     }
      95             : 
      96             :     ~GDALInvalidLocationLayer() override;
      97             : 
      98          13 :     bool TestCapability(const char *) const override
      99             :     {
     100          13 :         return false;
     101             :     }
     102             : 
     103           1 :     const char *GetDescription() const override
     104             :     {
     105           1 :         return GetName();
     106             :     }
     107             : 
     108         187 :     const OGRFeatureDefn *GetLayerDefn() const override
     109             :     {
     110         187 :         return m_defn.get();
     111             :     }
     112             : 
     113           6 :     std::unique_ptr<OGRFeature> CreateFeatureFromLastError() const
     114             :     {
     115           6 :         auto poErrorFeature = std::make_unique<OGRFeature>(m_defn.get());
     116             : 
     117          12 :         std::string msg = CPLGetLastErrorMsg();
     118             : 
     119             :         // Trim GEOS exception name
     120           6 :         const auto subMsgPos = msg.find(": ");
     121           6 :         if (subMsgPos != std::string::npos)
     122             :         {
     123           6 :             msg = msg.substr(subMsgPos + strlen(": "));
     124             :         }
     125             : 
     126             :         // Trim newline from end of GEOS exception message
     127           6 :         if (!msg.empty() && msg.back() == '\n')
     128             :         {
     129           3 :             msg.pop_back();
     130             :         }
     131             : 
     132           6 :         poErrorFeature->SetField(ERROR_DESCRIPTION_FIELD, msg.c_str());
     133             : 
     134           6 :         CPLErrorReset();
     135             : 
     136          12 :         return poErrorFeature;
     137             :     }
     138             : 
     139         354 :     bool TranslateFeature(
     140             :         std::unique_ptr<OGRFeature> poSrcFeature,
     141             :         std::vector<std::unique_ptr<OGRFeature>> &apoOutputFeatures) override
     142             :     {
     143             :         const OGRGeometry *poGeom =
     144         354 :             poSrcFeature->GetGeomFieldRef(m_srcGeomField);
     145           0 :         std::unique_ptr<OGRFeature> poErrorFeature;
     146             : 
     147         354 :         if (poGeom)
     148             :         {
     149         354 :             auto eType = wkbFlatten(poGeom->getGeometryType());
     150         354 :             GEOSGeometry *poGeosGeom = poGeom->exportToGEOS(m_geosContext);
     151             : 
     152         354 :             if (!poGeosGeom)
     153             :             {
     154             :                 // Try to find a useful message / coordinate from
     155             :                 // GEOS exception message.
     156           6 :                 poErrorFeature = CreateFeatureFromLastError();
     157             : 
     158           6 :                 if (eType == wkbPolygon)
     159             :                 {
     160             :                     const OGRLinearRing *poRing =
     161           4 :                         poGeom->toPolygon()->getExteriorRing();
     162           4 :                     if (poRing != nullptr && !poRing->IsEmpty())
     163             :                     {
     164           6 :                         auto poPoint = std::make_unique<OGRPoint>();
     165           3 :                         poRing->StartPoint(poPoint.get());
     166           3 :                         auto poMultiPoint = std::make_unique<OGRMultiPoint>();
     167           3 :                         poMultiPoint->addGeometry(std::move(poPoint));
     168           3 :                         poErrorFeature->SetGeometry(std::move(poMultiPoint));
     169             :                     }
     170             :                     else
     171             :                     {
     172             :                         // TODO get a point from somewhere else?
     173             :                     }
     174             :                 }
     175             :             }
     176             :             else
     177             :             {
     178         348 :                 char *pszReason = nullptr;
     179         348 :                 GEOSGeometry *location = nullptr;
     180         348 :                 char ret = 1;
     181         348 :                 bool warnAboutGeosVersion = false;
     182         348 :                 bool checkedSimple = false;
     183             : 
     184             :                 // check all geometry types for validity; IsSimple does not
     185             :                 // detect NaN/Inf coordinates
     186         348 :                 ret = GEOSisValidDetail_r(m_geosContext, poGeosGeom, 0,
     187             :                                           &pszReason, &location);
     188             : 
     189         348 :                 if (ret == 1 &&
     190         158 :                     (eType == wkbLineString || eType == wkbMultiLineString ||
     191         158 :                      eType == wkbCircularString || eType == wkbCompoundCurve ||
     192             :                      eType == wkbGeometryCollection))
     193             :                 {
     194          12 :                     checkedSimple = true;
     195             : #if GEOS_VERSION_MAJOR > 3 ||                                                  \
     196             :     (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 14)
     197          12 :                     ret = GEOSisSimpleDetail_r(m_geosContext, poGeosGeom, 1,
     198             :                                                &location);
     199             : #else
     200             :                     ret = GEOSisSimple_r(m_geosContext, poGeosGeom);
     201             :                     warnAboutGeosVersion = true;
     202             : #endif
     203             :                 }
     204             : 
     205         348 :                 GEOSGeom_destroy_r(m_geosContext, poGeosGeom);
     206         348 :                 if (ret == 0)
     207             :                 {
     208         190 :                     if (warnAboutGeosVersion)
     209             :                     {
     210           0 :                         CPLErrorOnce(
     211             :                             CE_Warning, CPLE_AppDefined,
     212             :                             "Detected a non-simple linear geometry, but "
     213             :                             "cannot output self-intersection points "
     214             :                             "because GEOS library version is < 3.14.");
     215             :                     }
     216             : 
     217         190 :                     poErrorFeature = std::make_unique<OGRFeature>(m_defn.get());
     218         190 :                     if (pszReason == nullptr)
     219             :                     {
     220           8 :                         if (checkedSimple)
     221             :                         {
     222           8 :                             poErrorFeature->SetField(ERROR_DESCRIPTION_FIELD,
     223             :                                                      "self-intersection");
     224             :                         }
     225             :                     }
     226             :                     else
     227             :                     {
     228         182 :                         poErrorFeature->SetField(ERROR_DESCRIPTION_FIELD,
     229             :                                                  pszReason);
     230         182 :                         GEOSFree_r(m_geosContext, pszReason);
     231             :                     }
     232             : 
     233         190 :                     if (location != nullptr)
     234             :                     {
     235             :                         std::unique_ptr<OGRGeometry> poErrorGeom(
     236         190 :                             OGRGeometryFactory::createFromGEOS(m_geosContext,
     237         190 :                                                                location));
     238         190 :                         GEOSGeom_destroy_r(m_geosContext, location);
     239             : 
     240         190 :                         if (poErrorGeom->getGeometryType() == wkbPoint)
     241             :                         {
     242             :                             auto poMultiPoint =
     243         372 :                                 std::make_unique<OGRMultiPoint>();
     244         186 :                             poMultiPoint->addGeometry(std::move(poErrorGeom));
     245         186 :                             poErrorGeom = std::move(poMultiPoint);
     246             :                         }
     247             : 
     248         380 :                         poErrorGeom->assignSpatialReference(
     249         380 :                             m_srcLayer.GetLayerDefn()
     250         190 :                                 ->GetGeomFieldDefn(m_srcGeomField)
     251         380 :                                 ->GetSpatialRef());
     252             : 
     253         190 :                         poErrorFeature->SetGeometry(std::move(poErrorGeom));
     254             :                     }
     255             :                 }
     256         158 :                 else if (ret == 2)
     257             :                 {
     258           0 :                     poErrorFeature = CreateFeatureFromLastError();
     259             :                 }
     260             :             }
     261             :         }
     262             : 
     263         354 :         if (!poErrorFeature && !m_skipValid)
     264             :         {
     265           6 :             poErrorFeature = std::make_unique<OGRFeature>(m_defn.get());
     266             :             // TODO Set geometry to POINT EMPTY ?
     267             :         }
     268             : 
     269         354 :         if (poErrorFeature)
     270             :         {
     271         202 :             if (!m_srcFieldMap.empty())
     272             :             {
     273           8 :                 poErrorFeature->SetFieldsFrom(
     274           4 :                     poSrcFeature.get(), m_srcFieldMap.data(), false, false);
     275             :             }
     276         202 :             poErrorFeature->SetFID(poSrcFeature->GetFID());
     277         202 :             apoOutputFeatures.push_back(std::move(poErrorFeature));
     278             :         }
     279             : 
     280         708 :         return true;
     281             :     }
     282             : 
     283             :     CPL_DISALLOW_COPY_ASSIGN(GDALInvalidLocationLayer)
     284             : 
     285             :   private:
     286             :     std::vector<int> m_srcFieldMap{};
     287             :     const OGRFeatureDefnRefCountedPtr m_defn;
     288             :     const GEOSContextHandle_t m_geosContext;
     289             :     const int m_srcGeomField;
     290             :     const bool m_skipValid;
     291             : };
     292             : 
     293          82 : GDALInvalidLocationLayer::~GDALInvalidLocationLayer()
     294             : {
     295          41 :     finishGEOS_r(m_geosContext);
     296          82 : }
     297             : 
     298             : #endif
     299             : 
     300          41 : bool GDALVectorCheckGeometryAlgorithm::RunStep(GDALPipelineStepRunContext &)
     301             : {
     302             : #ifdef HAVE_GEOS
     303          41 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     304          41 :     CPLAssert(poSrcDS);
     305          41 :     CPLAssert(m_outputDataset.GetName().empty());
     306          41 :     CPLAssert(!m_outputDataset.GetDatasetRef());
     307             : 
     308          41 :     const bool bSingleLayerOutput = m_inputLayerNames.empty()
     309          43 :                                         ? poSrcDS->GetLayerCount() == 1
     310           2 :                                         : m_inputLayerNames.size() == 1;
     311             : 
     312          82 :     auto outDS = std::make_unique<GDALVectorPipelineOutputDataset>(*poSrcDS);
     313          83 :     for (auto &&poSrcLayer : poSrcDS->GetLayers())
     314             :     {
     315          47 :         if (m_inputLayerNames.empty() ||
     316           0 :             std::find(m_inputLayerNames.begin(), m_inputLayerNames.end(),
     317          47 :                       poSrcLayer->GetDescription()) != m_inputLayerNames.end())
     318             :         {
     319          45 :             const auto poSrcLayerDefn = poSrcLayer->GetLayerDefn();
     320          45 :             if (poSrcLayerDefn->GetGeomFieldCount() == 0)
     321             :             {
     322           2 :                 if (m_inputLayerNames.empty())
     323           1 :                     continue;
     324           1 :                 ReportError(CE_Failure, CPLE_AppDefined,
     325             :                             "Specified layer '%s' has no geometry field",
     326           1 :                             poSrcLayer->GetDescription());
     327           3 :                 return false;
     328             :             }
     329             : 
     330             :             const int geomFieldIndex =
     331          43 :                 m_geomField.empty()
     332          43 :                     ? 0
     333           2 :                     : poSrcLayerDefn->GetGeomFieldIndex(m_geomField.c_str());
     334             : 
     335          43 :             if (geomFieldIndex == -1)
     336             :             {
     337           1 :                 ReportError(CE_Failure, CPLE_AppDefined,
     338             :                             "Specified geometry field '%s' does not exist in "
     339             :                             "layer '%s'",
     340           1 :                             m_geomField.c_str(), poSrcLayer->GetDescription());
     341           1 :                 return false;
     342             :             }
     343             : 
     344          42 :             std::vector<int> includeFieldIndices;
     345          42 :             if (!GetFieldIndices(m_includeFields,
     346             :                                  OGRLayer::ToHandle(poSrcLayer),
     347             :                                  includeFieldIndices))
     348             :             {
     349           1 :                 return false;
     350             :             }
     351             : 
     352          82 :             outDS->AddLayer(*poSrcLayer,
     353          41 :                             std::make_unique<GDALInvalidLocationLayer>(
     354             :                                 *poSrcLayer, includeFieldIndices,
     355             :                                 bSingleLayerOutput, geomFieldIndex,
     356          82 :                                 !m_includeValid));
     357             :         }
     358             :     }
     359             : 
     360          38 :     m_outputDataset.Set(std::move(outDS));
     361             : 
     362          38 :     return true;
     363             : #else
     364             :     ReportError(CE_Failure, CPLE_AppDefined,
     365             :                 "%s requires GDAL to be built against the GEOS library.", NAME);
     366             :     return false;
     367             : #endif
     368             : }
     369             : 
     370             : GDALVectorCheckGeometryAlgorithmStandalone::
     371             :     ~GDALVectorCheckGeometryAlgorithmStandalone() = default;
     372             : 
     373             : //! @endcond