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