/******************************************************************************

*

 * Project:  GDAL

 * Purpose:  "gdal vector check-geometry" subcommand

 * Author:   Daniel Baston

 *

 ******************************************************************************

 * Copyright (c) 2025, ISciences LLC

 *

 * SPDX-License-Identifier: MIT

 ****************************************************************************/

#include "gdalalg_vector_check_geometry.h"

#include "cpl_error.h"

#include "gdal_priv.h"

#include "gdalalg_vector_geom.h"

#include "ogr_geometry.h"

#include "ogr_geos.h"

#include <cinttypes>

#ifndef _

#define _(x) (x)

#endif

//! @cond Doxygen_Suppress

GDALVectorCheckGeometryAlgorithm::GDALVectorCheckGeometryAlgorithm(

    bool standaloneStep)

    : GDALVectorPipelineStepAlgorithm(

          NAME, DESCRIPTION, HELP_URL,

          ConstructorOptions()

              .SetStandaloneStep(standaloneStep)

              .SetNoCreateEmptyLayersArgument(standaloneStep))

{

    AddArg("include-field", 0,

           _("Fields from input layer to include in output (special values: "

             "ALL and NONE)"),

           &m_includeFields)

        .SetDefault("NONE");

    AddArg("include-valid", 0,

           _("Include valid inputs in output, with empty geometry"),

           &m_includeValid);

    AddArg("geometry-field", 0, _("Name of geometry field to check"),

           &m_geomField);

}

#ifdef HAVE_GEOS

class GDALInvalidLocationLayer final : public GDALVectorPipelineOutputLayer

{

  private:

    static constexpr const char *ERROR_DESCRIPTION_FIELD = "error";

  public:

    GDALInvalidLocationLayer(OGRLayer &layer,

                             const std::vector<int> &srcFieldIndices,

                             bool bSingleLayerOutput, int srcGeomField,

                             bool skipValid)

        : GDALVectorPipelineOutputLayer(layer),

          m_defn(OGRFeatureDefnRefCountedPtr::makeInstance(

              bSingleLayerOutput ? "error_location"

                                 : std::string("error_location_")

                                       .append(layer.GetDescription())

                                       .c_str())),

          m_geosContext(OGRGeometry::createGEOSContext()),

          m_srcGeomField(srcGeomField), m_skipValid(skipValid)

    {

        m_defn->SetGeomType(wkbMultiPoint);

        if (!srcFieldIndices.empty())

        {

            const OGRFeatureDefn &srcDefn = *layer.GetLayerDefn();

            m_srcFieldMap.resize(srcDefn.GetFieldCount(), -1);

            int iDstField = 0;

            for (int iSrcField : srcFieldIndices)

            {

                m_defn->AddFieldDefn(srcDefn.GetFieldDefn(iSrcField));

                m_srcFieldMap[iSrcField] = iDstField++;

            }

        }

        auto poDescriptionFieldDefn =

            std::make_unique<OGRFieldDefn>(ERROR_DESCRIPTION_FIELD, OFTString);

        m_defn->AddFieldDefn(std::move(poDescriptionFieldDefn));

        m_defn->GetGeomFieldDefn(0)->SetSpatialRef(

            m_srcLayer.GetLayerDefn()

                ->GetGeomFieldDefn(m_srcGeomField)

                ->GetSpatialRef());

    }

    ~GDALInvalidLocationLayer() override;

    bool TestCapability(const char *) const override

    {

        return false;

    }

    const char *GetDescription() const override

    {

        return GetName();

    }

    const OGRFeatureDefn *GetLayerDefn() const override

    {

        return m_defn.get();

    }

    std::unique_ptr<OGRFeature> CreateFeatureFromLastError() const

    {

        auto poErrorFeature = std::make_unique<OGRFeature>(m_defn.get());

        std::string msg = CPLGetLastErrorMsg();

        // Trim GEOS exception name

        const auto subMsgPos = msg.find(": ");

        if (subMsgPos != std::string::npos)

        {

            msg = msg.substr(subMsgPos + strlen(": "));

        }

        // Trim newline from end of GEOS exception message

        if (!msg.empty() && msg.back() == '\n')

        {

            msg.pop_back();

        }

        poErrorFeature->SetField(ERROR_DESCRIPTION_FIELD, msg.c_str());

        CPLErrorReset();

        return poErrorFeature;

    }

    bool TranslateFeature(

        std::unique_ptr<OGRFeature> poSrcFeature,

        std::vector<std::unique_ptr<OGRFeature>> &apoOutputFeatures) override

    {

        const OGRGeometry *poGeom =

            poSrcFeature->GetGeomFieldRef(m_srcGeomField);

        std::unique_ptr<OGRFeature> poErrorFeature;

        if (poGeom)

        {

            auto eType = wkbFlatten(poGeom->getGeometryType());

            GEOSGeometry *poGeosGeom = poGeom->exportToGEOS(m_geosContext);

            if (!poGeosGeom)

            {

                // Try to find a useful message / coordinate from

                // GEOS exception message.

                poErrorFeature = CreateFeatureFromLastError();

                if (eType == wkbPolygon)

                {

                    const OGRLinearRing *poRing =

                        poGeom->toPolygon()->getExteriorRing();

                    if (poRing != nullptr && !poRing->IsEmpty())

                    {

                        auto poPoint = std::make_unique<OGRPoint>();

                        poRing->StartPoint(poPoint.get());

                        auto poMultiPoint = std::make_unique<OGRMultiPoint>();

                        poMultiPoint->addGeometry(std::move(poPoint));

                        poErrorFeature->SetGeometry(std::move(poMultiPoint));

                    }

                    else

                    {

                        // TODO get a point from somewhere else?

                    }

                }

            }

            else

            {

                char *pszReason = nullptr;

                GEOSGeometry *location = nullptr;

                char ret = 1;

                bool warnAboutGeosVersion = false;

                bool checkedSimple = false;

                // check all geometry types for validity; IsSimple does not

                // detect NaN/Inf coordinates

                ret = GEOSisValidDetail_r(m_geosContext, poGeosGeom, 0,

                                          &pszReason, &location);

                if (ret == 1 &&

                    (eType == wkbLineString || eType == wkbMultiLineString ||

                     eType == wkbCircularString || eType == wkbCompoundCurve ||

                     eType == wkbGeometryCollection))

                {

                    checkedSimple = true;

#if GEOS_VERSION_MAJOR > 3 ||                                                  \

    (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 14)

                    ret = GEOSisSimpleDetail_r(m_geosContext, poGeosGeom, 1,

                                               &location);

#else

                    ret = GEOSisSimple_r(m_geosContext, poGeosGeom);

                    warnAboutGeosVersion = true;

#endif

                }

                GEOSGeom_destroy_r(m_geosContext, poGeosGeom);

                if (ret == 0)

                {

                    if (warnAboutGeosVersion)

                    {

                        CPLErrorOnce(

                            CE_Warning, CPLE_AppDefined,

                            "Detected a non-simple linear geometry, but "

                            "cannot output self-intersection points "

                            "because GEOS library version is < 3.14.");

                    }

                    poErrorFeature = std::make_unique<OGRFeature>(m_defn.get());

                    if (pszReason == nullptr)

                    {

                        if (checkedSimple)

                        {

                            poErrorFeature->SetField(ERROR_DESCRIPTION_FIELD,

                                                     "self-intersection");

                        }

                    }

                    else

                    {

                        poErrorFeature->SetField(ERROR_DESCRIPTION_FIELD,

                                                 pszReason);

                        GEOSFree_r(m_geosContext, pszReason);

                    }

                    if (location != nullptr)

                    {

                        std::unique_ptr<OGRGeometry> poErrorGeom(

                            OGRGeometryFactory::createFromGEOS(m_geosContext,

                                                               location));

                        GEOSGeom_destroy_r(m_geosContext, location);

                        if (poErrorGeom->getGeometryType() == wkbPoint)

                        {

                            auto poMultiPoint =

                                std::make_unique<OGRMultiPoint>();

                            poMultiPoint->addGeometry(std::move(poErrorGeom));

                            poErrorGeom = std::move(poMultiPoint);

                        }

                        poErrorGeom->assignSpatialReference(

                            m_srcLayer.GetLayerDefn()

                                ->GetGeomFieldDefn(m_srcGeomField)

                                ->GetSpatialRef());

                        poErrorFeature->SetGeometry(std::move(poErrorGeom));

                    }

                }

                else if (ret == 2)

                {

                    poErrorFeature = CreateFeatureFromLastError();

                }

            }

        }

        if (!poErrorFeature && !m_skipValid)

        {

            poErrorFeature = std::make_unique<OGRFeature>(m_defn.get());

            // TODO Set geometry to POINT EMPTY ?

        }

        if (poErrorFeature)

        {

            if (!m_srcFieldMap.empty())

            {

                poErrorFeature->SetFieldsFrom(

                    poSrcFeature.get(), m_srcFieldMap.data(), false, false);

            }

            poErrorFeature->SetFID(poSrcFeature->GetFID());

            if (PassesFilters(poErrorFeature.get()))

                apoOutputFeatures.push_back(std::move(poErrorFeature));

        }

        return true;

    }

    CPL_DISALLOW_COPY_ASSIGN(GDALInvalidLocationLayer)

  private:

    std::vector<int> m_srcFieldMap{};

    const OGRFeatureDefnRefCountedPtr m_defn;

    const GEOSContextHandle_t m_geosContext;

    const int m_srcGeomField;

    const bool m_skipValid;

};

GDALInvalidLocationLayer::~GDALInvalidLocationLayer()

{

    finishGEOS_r(m_geosContext);

}

#endif

bool GDALVectorCheckGeometryAlgorithm::RunStep(GDALPipelineStepRunContext &)

{

#ifdef HAVE_GEOS

    auto poSrcDS = m_inputDataset[0].GetDatasetRef();

    CPLAssert(poSrcDS);

    CPLAssert(m_outputDataset.GetName().empty());

    CPLAssert(!m_outputDataset.GetDatasetRef());

    const bool bSingleLayerOutput = m_inputLayerNames.empty()

                                        ? poSrcDS->GetLayerCount() == 1

                                        : m_inputLayerNames.size() == 1;

    auto outDS = std::make_unique<GDALVectorPipelineOutputDataset>(*poSrcDS);

    for (auto &&poSrcLayer : poSrcDS->GetLayers())

    {

        if (m_inputLayerNames.empty() ||

            std::find(m_inputLayerNames.begin(), m_inputLayerNames.end(),

                      poSrcLayer->GetDescription()) != m_inputLayerNames.end())

        {

            const auto poSrcLayerDefn = poSrcLayer->GetLayerDefn();

            if (poSrcLayerDefn->GetGeomFieldCount() == 0)

            {

                if (m_inputLayerNames.empty())

                    continue;

                ReportError(CE_Failure, CPLE_AppDefined,

                            "Specified layer '%s' has no geometry field",

                            poSrcLayer->GetDescription());

                return false;

            }

            const int geomFieldIndex =

                m_geomField.empty()

                    ? 0

                    : poSrcLayerDefn->GetGeomFieldIndex(m_geomField.c_str());

            if (geomFieldIndex == -1)

            {

                ReportError(CE_Failure, CPLE_AppDefined,

                            "Specified geometry field '%s' does not exist in "

                            "layer '%s'",

                            m_geomField.c_str(), poSrcLayer->GetDescription());

                return false;

            }

            std::vector<int> includeFieldIndices;

            if (!GetFieldIndices(m_includeFields,

                                 OGRLayer::ToHandle(poSrcLayer),

                                 includeFieldIndices))

            {

                return false;

            }

            outDS->AddLayer(*poSrcLayer,

                            std::make_unique<GDALInvalidLocationLayer>(

                                *poSrcLayer, includeFieldIndices,

                                bSingleLayerOutput, geomFieldIndex,

                                !m_includeValid));

        }

    }

    m_outputDataset.Set(std::move(outDS));

    return true;

#else

    ReportError(CE_Failure, CPLE_AppDefined,

                "%s requires GDAL to be built against the GEOS library.", NAME);

    return false;

#endif

}

GDALVectorCheckGeometryAlgorithmStandalone::

    ~GDALVectorCheckGeometryAlgorithmStandalone() = default;

//! @endcond