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