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

 *

 * Project:  GDAL

 * Purpose:  gdal "vector grid" subcommand

 * Author:   Even Rouault <even dot rouault at spatialys.com>

 *

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

 * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>

 *

 * SPDX-License-Identifier: MIT

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

#include "gdalalg_vector_grid.h"

#include "gdalalg_vector_grid_average.h"

#include "gdalalg_vector_grid_data_metrics.h"

#include "gdalalg_vector_grid_invdist.h"

#include "gdalalg_vector_grid_invdistnn.h"

#include "gdalalg_vector_grid_linear.h"

#include "gdalalg_vector_grid_nearest.h"

#include "cpl_conv.h"

#include "gdal_priv.h"

#include "gdal_utils.h"

#include "ogrsf_frmts.h"

//! @cond Doxygen_Suppress

#ifndef _

#define _(x) (x)

#endif

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

/*          GDALVectorGridAlgorithm::GDALVectorGridAlgorithm()          */

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

GDALVectorGridAlgorithm::GDALVectorGridAlgorithm(bool standaloneStep)

    : GDALPipelineStepAlgorithm(

          NAME, DESCRIPTION, HELP_URL,

          ConstructorOptions()

              .SetStandaloneStep(standaloneStep)

              .SetOutputFormatCreateCapability(GDAL_DCAP_CREATE))

{

    if (standaloneStep)

    {

        RegisterSubAlgorithm<GDALVectorGridAverageAlgorithmStandalone>();

        RegisterSubAlgorithm<GDALVectorGridInvdistAlgorithmStandalone>();

        RegisterSubAlgorithm<GDALVectorGridInvdistNNAlgorithmStandalone>();

        RegisterSubAlgorithm<GDALVectorGridLinearAlgorithmStandalone>();

        RegisterSubAlgorithm<GDALVectorGridNearestAlgorithmStandalone>();

        RegisterSubAlgorithm<GDALVectorGridMinimumAlgorithmStandalone>();

        RegisterSubAlgorithm<GDALVectorGridMaximumAlgorithmStandalone>();

        RegisterSubAlgorithm<GDALVectorGridRangeAlgorithmStandalone>();

        RegisterSubAlgorithm<GDALVectorGridCountAlgorithmStandalone>();

        RegisterSubAlgorithm<

            GDALVectorGridAverageDistanceAlgorithmStandalone>();

        RegisterSubAlgorithm<

            GDALVectorGridAverageDistancePointsAlgorithmStandalone>();

    }

    else

    {

        RegisterSubAlgorithm<GDALVectorGridAverageAlgorithm>();

        RegisterSubAlgorithm<GDALVectorGridInvdistAlgorithm>();

        RegisterSubAlgorithm<GDALVectorGridInvdistNNAlgorithm>();

        RegisterSubAlgorithm<GDALVectorGridLinearAlgorithm>();

        RegisterSubAlgorithm<GDALVectorGridNearestAlgorithm>();

        RegisterSubAlgorithm<GDALVectorGridMinimumAlgorithm>();

        RegisterSubAlgorithm<GDALVectorGridMaximumAlgorithm>();

        RegisterSubAlgorithm<GDALVectorGridRangeAlgorithm>();

        RegisterSubAlgorithm<GDALVectorGridCountAlgorithm>();

        RegisterSubAlgorithm<GDALVectorGridAverageDistanceAlgorithm>();

        RegisterSubAlgorithm<GDALVectorGridAverageDistancePointsAlgorithm>();

    }

}

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

/*                  GDALVectorGridAlgorithm::RunStep()                  */

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

bool GDALVectorGridAlgorithm::RunStep(GDALPipelineStepRunContext &)

{

    CPLError(CE_Failure, CPLE_AppDefined,

             "The Run() method should not be called directly on the \"gdal "

             "vector grid\" program.");

    return false;

}

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

/*                  GDALVectorGeomAlgorithm::RunImpl()                  */

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

bool GDALVectorGridAlgorithm::RunImpl(GDALProgressFunc, void *)

{

    CPLError(CE_Failure, CPLE_AppDefined,

             "The Run() method should not be called directly on the \"gdal "

             "vector grid\" program.");

    return false;

}

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

/*                  GDALVectorGridAbstractAlgorithm()                   */

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

GDALVectorGridAbstractAlgorithm::GDALVectorGridAbstractAlgorithm(

    const std::string &name, const std::string &description,

    const std::string &helpURL, bool standaloneStep)

    : GDALPipelineStepAlgorithm(

          name, description, helpURL,

          ConstructorOptions()

              .SetStandaloneStep(standaloneStep)

              .SetOutputFormatCreateCapability(GDAL_DCAP_CREATE))

{

    AddProgressArg();

    if (standaloneStep)

    {

        AddOutputFormatArg(&m_format).AddMetadataItem(

            GAAMDI_REQUIRED_CAPABILITIES, {GDAL_DCAP_RASTER, GDAL_DCAP_CREATE});

        AddOpenOptionsArg(&m_openOptions);

        AddInputFormatsArg(&m_inputFormats)

            .AddMetadataItem(GAAMDI_REQUIRED_CAPABILITIES, {GDAL_DCAP_VECTOR});

        AddInputDatasetArg(&m_inputDataset, GDAL_OF_VECTOR);

        AddOutputDatasetArg(&m_outputDataset, GDAL_OF_RASTER);

        AddCreationOptionsArg(&m_creationOptions);

    }

    else

    {

        AddVectorHiddenInputDatasetArg();

    }

    AddArg("extent", 0, _("Set the target georeferenced extent"),

           &m_targetExtent)

        .SetMinCount(4)

        .SetMaxCount(4)

        .SetRepeatedArgAllowed(false)

        .SetMetaVar("<xmin>,<ymin>,<xmax>,<ymax>");

    AddArg("resolution", 0, _("Set the target resolution"), &m_targetResolution)

        .SetMinCount(2)

        .SetMaxCount(2)

        .SetRepeatedArgAllowed(false)

        .SetMetaVar("<xres>,<yres>")

        .SetMutualExclusionGroup("size-or-resolution");

    AddArg("size", 0, _("Set the target size in pixels and lines"),

           &m_targetSize)

        .SetMinCount(2)

        .SetMaxCount(2)

        .SetRepeatedArgAllowed(false)

        .SetMetaVar("<xsize>,<ysize>")

        .SetMutualExclusionGroup("size-or-resolution");

    AddOutputDataTypeArg(&m_outputType).SetDefault(m_outputType);

    AddArg("crs", 0, _("Override the projection for the output file"), &m_crs)

        .AddHiddenAlias("srs")

        .SetIsCRSArg(/*noneAllowed=*/false);

    AddOverwriteArg(&m_overwrite);

    AddLayerNameArg(&m_layers)

        .SetMutualExclusionGroup("layer-sql")

        .AddAlias("layer");

    AddArg("sql", 0, _("SQL statement"), &m_sql)

        .SetReadFromFileAtSyntaxAllowed()

        .SetMetaVar("<statement>|@<filename>")

        .SetRemoveSQLCommentsEnabled()

        .SetMutualExclusionGroup("layer-sql");

    AddBBOXArg(

        &m_bbox,

        _("Select only points contained within the specified bounding box"));

    AddArg("zfield", 0, _("Field name from which to get Z values."), &m_zField)

        .AddHiddenAlias("z-field");

    AddArg("zoffset", 0,

           _("Value to add to the Z field value (applied before zmultiply)"),

           &m_zOffset)

        .SetDefault(m_zOffset)

        .AddHiddenAlias("z-offset");

    AddArg("zmultiply", 0,

           _("Multiplication factor for the Z field value (applied after "

             "zoffset)"),

           &m_zMultiply)

        .SetDefault(m_zMultiply)

        .AddHiddenAlias("z-multiply");

    AddValidationAction(

        [this]()

        {

            bool ret = true;

            if (!m_targetResolution.empty() && m_targetExtent.empty())

            {

                ReportError(CE_Failure, CPLE_AppDefined,

                            "'resolution' should be defined when 'extent' is.");

                ret = false;

            }

            return ret;

        });

}

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

/*              GDALVectorGridAbstractAlgorithm::RunStep()              */

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

bool GDALVectorGridAbstractAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)

{

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

    CPLAssert(poSrcDS);

    CPLAssert(!m_outputDataset.GetDatasetRef());

    CPLStringList aosOptions;

    std::string outputFilename;

    if (m_standaloneStep)

    {

        outputFilename = m_outputDataset.GetName();

        if (!m_format.empty())

        {

            aosOptions.AddString("-of");

            aosOptions.AddString(m_format.c_str());

        }

        for (const std::string &co : m_creationOptions)

        {

            aosOptions.AddString("-co");

            aosOptions.AddString(co.c_str());

        }

    }

    else

    {

        outputFilename = CPLGenerateTempFilenameSafe("_grid.tif");

        aosOptions.AddString("-of");

        aosOptions.AddString("GTiff");

        aosOptions.AddString("-co");

        aosOptions.AddString("TILED=YES");

    }

    if (!m_targetExtent.empty())

    {

        aosOptions.AddString("-txe");

        aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[0]));

        aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[2]));

        aosOptions.AddString("-tye");

        aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[1]));

        aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[3]));

    }

    if (!m_bbox.empty())

    {

        aosOptions.AddString("-clipsrc");

        for (double v : m_bbox)

        {

            aosOptions.AddString(CPLSPrintf("%.17g", v));

        }

    }

    if (!m_targetResolution.empty())

    {

        aosOptions.AddString("-tr");

        for (double targetResolution : m_targetResolution)

        {

            aosOptions.AddString(CPLSPrintf("%.17g", targetResolution));

        }

    }

    if (!m_targetSize.empty())

    {

        aosOptions.AddString("-outsize");

        for (int targetSize : m_targetSize)

        {

            aosOptions.AddString(CPLSPrintf("%d", targetSize));

        }

    }

    if (!m_outputType.empty())

    {

        aosOptions.AddString("-ot");

        aosOptions.AddString(m_outputType.c_str());

    }

    if (!m_crs.empty())

    {

        aosOptions.AddString("-a_srs");

        aosOptions.AddString(m_crs.c_str());

    }

    if (m_sql.empty())

    {

        for (const std::string &layer : m_layers)

        {

            aosOptions.AddString("-l");

            aosOptions.AddString(layer.c_str());

        }

    }

    else

    {

        aosOptions.AddString("-sql");

        aosOptions.AddString(m_sql.c_str());

    }

    if (m_zOffset != 0)

    {

        aosOptions.AddString("-z_increase");

        aosOptions.AddString(CPLSPrintf("%.17g", m_zOffset));

    }

    if (m_zMultiply != 0)

    {

        aosOptions.AddString("-z_multiply");

        aosOptions.AddString(CPLSPrintf("%.17g", m_zMultiply));

    }

    if (!m_zField.empty())

    {

        aosOptions.AddString("-zfield");

        aosOptions.AddString(m_zField.c_str());

    }

    else if (m_sql.empty())

    {

        const auto CheckLayer = [this](OGRLayer *poLayer)

        {

            auto poFeat =

                std::unique_ptr<OGRFeature>(poLayer->GetNextFeature());

            poLayer->ResetReading();

            if (poFeat)

            {

                const auto poGeom = poFeat->GetGeometryRef();

                if (!poGeom || !poGeom->Is3D())

                {

                    ReportError(

                        CE_Warning, CPLE_AppDefined,

                        "At least one geometry of layer '%s' lacks a Z "

                        "component. You may need to set the 'zfield' argument",

                        poLayer->GetName());

                    return false;

                }

            }

            return true;

        };

        if (m_layers.empty())

        {

            for (auto poLayer : poSrcDS->GetLayers())

            {

                if (!CheckLayer(poLayer))

                    break;

            }

        }

        else

        {

            for (const std::string &layerName : m_layers)

            {

                auto poLayer = poSrcDS->GetLayerByName(layerName.c_str());

                if (poLayer && !CheckLayer(poLayer))

                    break;

            }

        }

    }

    aosOptions.AddString("-a");

    aosOptions.AddString(GetGridAlgorithm().c_str());

    bool bOK = false;

    std::unique_ptr<GDALGridOptions, decltype(&GDALGridOptionsFree)> psOptions{

        GDALGridOptionsNew(aosOptions.List(), nullptr), GDALGridOptionsFree};

    if (psOptions)

    {

        GDALGridOptionsSetProgress(psOptions.get(), ctxt.m_pfnProgress,

                                   ctxt.m_pProgressData);

        auto poRetDS = std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(

            GDALGrid(outputFilename.c_str(), GDALDataset::ToHandle(poSrcDS),

                     psOptions.get(), nullptr)));

        if (poRetDS)

        {

            bOK = true;

            if (!m_standaloneStep)

            {

                VSIUnlink(outputFilename.c_str());

                poRetDS->MarkSuppressOnClose();

            }

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

        }

    }

    return bOK;

}

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

/*              GDALVectorGridAbstractAlgorithm::RunImpl()              */

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

bool GDALVectorGridAbstractAlgorithm::RunImpl(GDALProgressFunc pfnProgress,

                                              void *pProgressData)

{

    GDALPipelineStepRunContext stepCtxt;

    stepCtxt.m_pfnProgress = pfnProgress;

    stepCtxt.m_pProgressData = pProgressData;

    return RunPreStepPipelineValidations() && RunStep(stepCtxt);

}

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

/*           GDALVectorGridAbstractAlgorithm::AddRadiusArg()            */

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

GDALInConstructionAlgorithmArg &GDALVectorGridAbstractAlgorithm::AddRadiusArg()

{

    return AddArg("radius", 0, _("Radius of the search circle"), &m_radius)

        .SetMutualExclusionGroup("radius");

}

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

/*      GDALVectorGridAbstractAlgorithm::AddRadius1AndRadius2Arg()      */

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

void GDALVectorGridAbstractAlgorithm::AddRadius1AndRadius2Arg()

{

    AddArg("radius1", 0, _("First axis of the search ellipse"), &m_radius1)

        .SetMutualExclusionGroup("radius");

    AddArg("radius2", 0, _("Second axis of the search ellipse"), &m_radius2);

    AddValidationAction(

        [this]()

        {

            bool ret = true;

            if (m_radius1 > 0 && m_radius2 == 0)

            {

                ReportError(CE_Failure, CPLE_AppDefined,

                            "'radius2' should be defined when 'radius1' is.");

                ret = false;

            }

            else if (m_radius2 > 0 && m_radius1 == 0)

            {

                ReportError(CE_Failure, CPLE_AppDefined,

                            "'radius1' should be defined when 'radius2' is.");

                ret = false;

            }

            return ret;

        });

}

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

/*            GDALVectorGridAbstractAlgorithm::AddAngleArg()            */

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

GDALInConstructionAlgorithmArg &GDALVectorGridAbstractAlgorithm::AddAngleArg()

{

    return AddArg("angle", 0,

                  _("Angle of search ellipse rotation in degrees (counter "

                    "clockwise)"),

                  &m_angle)

        .SetDefault(m_angle);

}

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

/*          GDALVectorGridAbstractAlgorithm::AddMinPointsArg()          */

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

GDALInConstructionAlgorithmArg &

GDALVectorGridAbstractAlgorithm::AddMinPointsArg()

{

    return AddArg("min-points", 0, _("Minimum number of data points to use"),

                  &m_minPoints)

        .SetDefault(m_minPoints);

}

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

/*          GDALVectorGridAbstractAlgorithm::AddMaxPointsArg()          */

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

GDALInConstructionAlgorithmArg &

GDALVectorGridAbstractAlgorithm::AddMaxPointsArg()

{

    return AddArg("max-points", 0, _("Maximum number of data points to use"),

                  &m_maxPoints)

        .SetDefault(m_maxPoints);

}

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

/*   GDALVectorGridAbstractAlgorithm::AddMinMaxPointsPerQuadrantArg()   */

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

void GDALVectorGridAbstractAlgorithm::AddMinMaxPointsPerQuadrantArg()

{

    AddArg("min-points-per-quadrant", 0,

           _("Minimum number of data points to use per quadrant"),

           &m_minPointsPerQuadrant)

        .SetDefault(m_minPointsPerQuadrant);

    AddArg("max-points-per-quadrant", 0,

           _("Maximum number of data points to use per quadrant"),

           &m_maxPointsPerQuadrant)

        .SetDefault(m_maxPointsPerQuadrant);

}

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

/*           GDALVectorGridAbstractAlgorithm::AddNodataArg()            */

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

GDALInConstructionAlgorithmArg &GDALVectorGridAbstractAlgorithm::AddNodataArg()

{

    return AddArg("nodata", 0, _("Target nodata value"), &m_nodata)

        .SetDefault(m_nodata);

}

GDALVectorGridAlgorithmStandalone::~GDALVectorGridAlgorithmStandalone() =

    default;

//! @endcond