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

 *

 * Project:  GDAL Utilities

 * Purpose:  GDAL scattered data gridding (interpolation) tool

 * Author:   Andrey Kiselev, dron@ak4719.spb.edu

 *

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

 * Copyright (c) 2007, Andrey Kiselev <dron@ak4719.spb.edu>

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

 *

 * SPDX-License-Identifier: MIT

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

#include "cpl_port.h"

#include "gdal_utils.h"

#include "gdal_utils_priv.h"

#include "commonutils.h"

#include "gdalargumentparser.h"

#include <cmath>

#include <cstdint>

#include <cstdio>

#include <cstdlib>

#include <algorithm>

#include <vector>

#include "cpl_conv.h"

#include "cpl_error.h"

#include "cpl_progress.h"

#include "cpl_string.h"

#include "cpl_vsi.h"

#include "gdal.h"

#include "gdal_alg.h"

#include "gdal_priv.h"

#include "gdalgrid.h"

#include "ogr_api.h"

#include "ogr_core.h"

#include "ogr_feature.h"

#include "ogr_geometry.h"

#include "ogr_spatialref.h"

#include "ogr_srs_api.h"

#include "ogrsf_frmts.h"

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

/*                           GDALGridOptions                            */

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

/** Options for use with GDALGrid(). GDALGridOptions* must be allocated

 * and freed with GDALGridOptionsNew() and GDALGridOptionsFree() respectively.

 */

struct GDALGridOptions

{

    /*! output format. Use the short format name. */

    std::string osFormat{};

    /*! allow or suppress progress monitor and other non-error output */

    bool bQuiet = true;

    /*! the progress function to use */

    GDALProgressFunc pfnProgress = GDALDummyProgress;

    /*! pointer to the progress data variable */

    void *pProgressData = nullptr;

    CPLStringList aosLayers{};

    std::string osBurnAttribute{};

    double dfIncreaseBurnValue = 0.0;

    double dfMultiplyBurnValue = 1.0;

    std::string osWHERE{};

    std::string osSQL{};

    GDALDataType eOutputType = GDT_Float64;

    CPLStringList aosCreateOptions{};

    int nXSize = 0;

    int nYSize = 0;

    double dfXRes = 0;

    double dfYRes = 0;

    double dfXMin = 0;

    double dfXMax = 0;

    double dfYMin = 0;

    double dfYMax = 0;

    bool bIsXExtentSet = false;

    bool bIsYExtentSet = false;

    GDALGridAlgorithm eAlgorithm = GGA_InverseDistanceToAPower;

    std::unique_ptr<void, VSIFreeReleaser> pOptions{};

    std::string osOutputSRS{};

    std::unique_ptr<OGRGeometry> poSpatialFilter{};

    bool bClipSrc = false;

    std::unique_ptr<OGRGeometry> poClipSrc{};

    std::string osClipSrcDS{};

    std::string osClipSrcSQL{};

    std::string osClipSrcLayer{};

    std::string osClipSrcWhere{};

    bool bNoDataSet = false;

    double dfNoDataValue = 0;

    GDALGridOptions()

    {

        void *l_pOptions = nullptr;

        GDALGridParseAlgorithmAndOptions(szAlgNameInvDist, &eAlgorithm,

                                         &l_pOptions);

        pOptions.reset(l_pOptions);

    }

    CPL_DISALLOW_COPY_ASSIGN(GDALGridOptions)

};

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

/*                          GetAlgorithmName()                          */

/*                                                                      */

/*      Grids algorithm code into mnemonic name.                        */

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

static void PrintAlgorithmAndOptions(GDALGridAlgorithm eAlgorithm,

                                     void *pOptions)

{

    switch (eAlgorithm)

    {

        case GGA_InverseDistanceToAPower:

        {

            printf("Algorithm name: \"%s\".\n", szAlgNameInvDist);

            GDALGridInverseDistanceToAPowerOptions *pOptions2 =

                static_cast<GDALGridInverseDistanceToAPowerOptions *>(pOptions);

            CPLprintf("Options are "

                      "\"power=%f:smoothing=%f:radius1=%f:radius2=%f:angle=%f"

                      ":max_points=%u:min_points=%u:nodata=%f\"\n",

                      pOptions2->dfPower, pOptions2->dfSmoothing,

                      pOptions2->dfRadius1, pOptions2->dfRadius2,

                      pOptions2->dfAngle, pOptions2->nMaxPoints,

                      pOptions2->nMinPoints, pOptions2->dfNoDataValue);

            break;

        }

        case GGA_InverseDistanceToAPowerNearestNeighbor:

        {

            printf("Algorithm name: \"%s\".\n",

                   szAlgNameInvDistNearestNeighbor);

            GDALGridInverseDistanceToAPowerNearestNeighborOptions *pOptions2 =

                static_cast<

                    GDALGridInverseDistanceToAPowerNearestNeighborOptions *>(

                    pOptions);

            CPLString osStr;

            osStr.Printf("power=%f:smoothing=%f:radius=%f"

                         ":max_points=%u:min_points=%u:nodata=%f",

                         pOptions2->dfPower, pOptions2->dfSmoothing,

                         pOptions2->dfRadius, pOptions2->nMaxPoints,

                         pOptions2->nMinPoints, pOptions2->dfNoDataValue);

            if (pOptions2->nMinPointsPerQuadrant > 0)

                osStr += CPLSPrintf(":min_points_per_quadrant=%u",

                                    pOptions2->nMinPointsPerQuadrant);

            if (pOptions2->nMaxPointsPerQuadrant > 0)

                osStr += CPLSPrintf(":max_points_per_quadrant=%u",

                                    pOptions2->nMaxPointsPerQuadrant);

            printf("Options are: \"%s\n", osStr.c_str()); /* ok */

            break;

        }

        case GGA_MovingAverage:

        {

            printf("Algorithm name: \"%s\".\n", szAlgNameAverage);

            GDALGridMovingAverageOptions *pOptions2 =

                static_cast<GDALGridMovingAverageOptions *>(pOptions);

            CPLString osStr;

            osStr.Printf("radius1=%f:radius2=%f:angle=%f:min_points=%u"

                         ":nodata=%f",

                         pOptions2->dfRadius1, pOptions2->dfRadius2,

                         pOptions2->dfAngle, pOptions2->nMinPoints,

                         pOptions2->dfNoDataValue);

            if (pOptions2->nMinPointsPerQuadrant > 0)

                osStr += CPLSPrintf(":min_points_per_quadrant=%u",

                                    pOptions2->nMinPointsPerQuadrant);

            if (pOptions2->nMaxPointsPerQuadrant > 0)

                osStr += CPLSPrintf(":max_points_per_quadrant=%u",

                                    pOptions2->nMaxPointsPerQuadrant);

            if (pOptions2->nMaxPoints > 0)

                osStr += CPLSPrintf(":max_points=%u", pOptions2->nMaxPoints);

            printf("Options are: \"%s\n", osStr.c_str()); /* ok */

            break;

        }

        case GGA_NearestNeighbor:

        {

            printf("Algorithm name: \"%s\".\n", szAlgNameNearest);

            GDALGridNearestNeighborOptions *pOptions2 =

                static_cast<GDALGridNearestNeighborOptions *>(pOptions);

            CPLprintf("Options are "

                      "\"radius1=%f:radius2=%f:angle=%f:nodata=%f\"\n",

                      pOptions2->dfRadius1, pOptions2->dfRadius2,

                      pOptions2->dfAngle, pOptions2->dfNoDataValue);

            break;

        }

        case GGA_MetricMinimum:

        case GGA_MetricMaximum:

        case GGA_MetricRange:

        case GGA_MetricCount:

        case GGA_MetricAverageDistance:

        case GGA_MetricAverageDistancePts:

        {

            const char *pszAlgName = "";

            CPL_IGNORE_RET_VAL(pszAlgName);  // Make CSA happy

            switch (eAlgorithm)

            {

                case GGA_MetricMinimum:

                    pszAlgName = szAlgNameMinimum;

                    break;

                case GGA_MetricMaximum:

                    pszAlgName = szAlgNameMaximum;

                    break;

                case GGA_MetricRange:

                    pszAlgName = szAlgNameRange;

                    break;

                case GGA_MetricCount:

                    pszAlgName = szAlgNameCount;

                    break;

                case GGA_MetricAverageDistance:

                    pszAlgName = szAlgNameAverageDistance;

                    break;

                case GGA_MetricAverageDistancePts:

                    pszAlgName = szAlgNameAverageDistancePts;

                    break;

                default:

                    CPLAssert(false);

                    break;

            }

            printf("Algorithm name: \"%s\".\n", pszAlgName);

            GDALGridDataMetricsOptions *pOptions2 =

                static_cast<GDALGridDataMetricsOptions *>(pOptions);

            CPLString osStr;

            osStr.Printf("radius1=%f:radius2=%f:angle=%f:min_points=%u"

                         ":nodata=%f",

                         pOptions2->dfRadius1, pOptions2->dfRadius2,

                         pOptions2->dfAngle, pOptions2->nMinPoints,

                         pOptions2->dfNoDataValue);

            if (pOptions2->nMinPointsPerQuadrant > 0)

                osStr += CPLSPrintf(":min_points_per_quadrant=%u",

                                    pOptions2->nMinPointsPerQuadrant);

            if (pOptions2->nMaxPointsPerQuadrant > 0)

                osStr += CPLSPrintf(":max_points_per_quadrant=%u",

                                    pOptions2->nMaxPointsPerQuadrant);

            printf("Options are: \"%s\n", osStr.c_str()); /* ok */

            break;

        }

        case GGA_Linear:

        {

            printf("Algorithm name: \"%s\".\n", szAlgNameLinear);

            GDALGridLinearOptions *pOptions2 =

                static_cast<GDALGridLinearOptions *>(pOptions);

            CPLprintf("Options are "

                      "\"radius=%f:nodata=%f\"\n",

                      pOptions2->dfRadius, pOptions2->dfNoDataValue);

            break;

        }

        default:

        {

            printf("Algorithm is unknown.\n");

            break;

        }

    }

}

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

/*  Extract point coordinates from the geometry reference and set the   */

/*  Z value as requested. Test whether we are in the clipped region     */

/*  before processing.                                                  */

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

class GDALGridGeometryVisitor final : public OGRDefaultConstGeometryVisitor

{

  public:

    const OGRGeometry *poClipSrc = nullptr;

    int iBurnField = 0;

    double dfBurnValue = 0;

    double dfIncreaseBurnValue = 0;

    double dfMultiplyBurnValue = 1;

    std::vector<double> adfX{};

    std::vector<double> adfY{};

    std::vector<double> adfZ{};

    using OGRDefaultConstGeometryVisitor::visit;

    void visit(const OGRPoint *p) override;

};

void GDALGridGeometryVisitor::visit(const OGRPoint *p)

{

    if (poClipSrc && !p->Within(poClipSrc))

        return;

    if (iBurnField < 0 && std::isnan(p->getZ()))

        return;

    adfX.push_back(p->getX());

    adfY.push_back(p->getY());

    if (iBurnField < 0)

        adfZ.push_back((p->getZ() + dfIncreaseBurnValue) * dfMultiplyBurnValue);

    else

        adfZ.push_back((dfBurnValue + dfIncreaseBurnValue) *

                       dfMultiplyBurnValue);

}

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

/*                            ProcessLayer()                            */

/*                                                                      */

/*      Process all the features in a layer selection, collecting       */

/*      geometries and burn values.                                     */

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

static CPLErr ProcessLayer(OGRLayer *poSrcLayer, GDALDataset *poDstDS,

                           const OGRGeometry *poClipSrc, int nXSize, int nYSize,

                           int nBand, bool &bIsXExtentSet, bool &bIsYExtentSet,

                           double &dfXMin, double &dfXMax, double &dfYMin,

                           double &dfYMax, const std::string &osBurnAttribute,

                           const double dfIncreaseBurnValue,

                           const double dfMultiplyBurnValue, GDALDataType eType,

                           GDALGridAlgorithm eAlgorithm, void *pOptions,

                           bool bQuiet, GDALProgressFunc pfnProgress,

                           void *pProgressData)

{

    /* -------------------------------------------------------------------- */

    /*      Get field index, and check.                                     */

    /* -------------------------------------------------------------------- */

    int iBurnField = -1;

    if (!osBurnAttribute.empty())

    {

        iBurnField =

            poSrcLayer->GetLayerDefn()->GetFieldIndex(osBurnAttribute.c_str());

        if (iBurnField == -1)

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "Failed to find field %s on layer %s.",

                     osBurnAttribute.c_str(), poSrcLayer->GetName());

            return CE_Failure;

        }

    }

    /* -------------------------------------------------------------------- */

    /*      Collect the geometries from this layer, and build list of       */

    /*      values to be interpolated.                                      */

    /* -------------------------------------------------------------------- */

    GDALGridGeometryVisitor oVisitor;

    oVisitor.poClipSrc = poClipSrc;

    oVisitor.iBurnField = iBurnField;

    oVisitor.dfIncreaseBurnValue = dfIncreaseBurnValue;

    oVisitor.dfMultiplyBurnValue = dfMultiplyBurnValue;

    for (auto &&poFeat : poSrcLayer)

    {

        const OGRGeometry *poGeom = poFeat->GetGeometryRef();

        if (poGeom)

        {

            if (iBurnField >= 0)

            {

                if (!poFeat->IsFieldSetAndNotNull(iBurnField))

                {

                    continue;

                }

                oVisitor.dfBurnValue = poFeat->GetFieldAsDouble(iBurnField);

            }

            poGeom->accept(&oVisitor);

        }

    }

    if (oVisitor.adfX.empty())

    {

        CPLError(CE_Warning, CPLE_AppDefined,

                 "No point geometry found on layer %s, skipping.",

                 poSrcLayer->GetName());

        return CE_None;

    }

    /* -------------------------------------------------------------------- */

    /*      Compute grid geometry.                                          */

    /* -------------------------------------------------------------------- */

    if (!bIsXExtentSet || !bIsYExtentSet)

    {

        OGREnvelope sEnvelope;

        if (poSrcLayer->GetExtent(&sEnvelope, TRUE) == OGRERR_FAILURE)

        {

            return CE_Failure;

        }

        if (!bIsXExtentSet)

        {

            dfXMin = sEnvelope.MinX;

            dfXMax = sEnvelope.MaxX;

            bIsXExtentSet = true;

        }

        if (!bIsYExtentSet)

        {

            dfYMin = sEnvelope.MinY;

            dfYMax = sEnvelope.MaxY;

            bIsYExtentSet = true;

        }

    }

    // Produce north-up images

    if (dfYMin < dfYMax)

        std::swap(dfYMin, dfYMax);

    /* -------------------------------------------------------------------- */

    /*      Perform gridding.                                               */

    /* -------------------------------------------------------------------- */

    const double dfDeltaX = (dfXMax - dfXMin) / nXSize;

    const double dfDeltaY = (dfYMax - dfYMin) / nYSize;

    if (!bQuiet)

    {

        printf("Grid data type is \"%s\"\n", GDALGetDataTypeName(eType));

        printf("Grid size = (%d %d).\n", nXSize, nYSize);

        CPLprintf("Corner coordinates = (%f %f)-(%f %f).\n", dfXMin, dfYMin,

                  dfXMax, dfYMax);

        CPLprintf("Grid cell size = (%f %f).\n", dfDeltaX, dfDeltaY);

        printf("Source point count = %lu.\n",

               static_cast<unsigned long>(oVisitor.adfX.size()));

        PrintAlgorithmAndOptions(eAlgorithm, pOptions);

        printf("\n");

    }

    GDALRasterBand *poBand = poDstDS->GetRasterBand(nBand);

    int nBlockXSize = 0;

    int nBlockYSize = 0;

    const int nDataTypeSize = GDALGetDataTypeSizeBytes(eType);

    // Try to grow the work buffer up to 16 MB if it is smaller

    poBand->GetBlockSize(&nBlockXSize, &nBlockYSize);

    if (nXSize == 0 || nYSize == 0 || nBlockXSize == 0 || nBlockYSize == 0)

        return CE_Failure;

    const int nDesiredBufferSize = 16 * 1024 * 1024;

    if (nBlockXSize < nXSize && nBlockYSize < nYSize &&

        nBlockXSize < nDesiredBufferSize / (nBlockYSize * nDataTypeSize))

    {

        const int nNewBlockXSize =

            nDesiredBufferSize / (nBlockYSize * nDataTypeSize);

        nBlockXSize = (nNewBlockXSize / nBlockXSize) * nBlockXSize;

        if (nBlockXSize > nXSize)

            nBlockXSize = nXSize;

    }

    else if (nBlockXSize == nXSize && nBlockYSize < nYSize &&

             nBlockYSize < nDesiredBufferSize / (nXSize * nDataTypeSize))

    {

        const int nNewBlockYSize =

            nDesiredBufferSize / (nXSize * nDataTypeSize);

        nBlockYSize = (nNewBlockYSize / nBlockYSize) * nBlockYSize;

        if (nBlockYSize > nYSize)

            nBlockYSize = nYSize;

    }

    CPLDebug("GDAL_GRID", "Work buffer: %d * %d", nBlockXSize, nBlockYSize);

    std::unique_ptr<void, VSIFreeReleaser> pData(

        VSIMalloc3(nBlockXSize, nBlockYSize, nDataTypeSize));

    if (!pData)

    {

        CPLError(CE_Failure, CPLE_OutOfMemory, "Cannot allocate work buffer");

        return CE_Failure;

    }

    GIntBig nBlock = 0;

    const double dfBlockCount =

        static_cast<double>(DIV_ROUND_UP(nXSize, nBlockXSize)) *

        DIV_ROUND_UP(nYSize, nBlockYSize);

    struct GDALGridContextReleaser

    {

        void operator()(GDALGridContext *psContext)

        {

            GDALGridContextFree(psContext);

        }

    };

    std::unique_ptr<GDALGridContext, GDALGridContextReleaser> psContext(

        GDALGridContextCreate(eAlgorithm, pOptions,

                              static_cast<int>(oVisitor.adfX.size()),

                              &(oVisitor.adfX[0]), &(oVisitor.adfY[0]),

                              &(oVisitor.adfZ[0]), TRUE));

    if (!psContext)

    {

        return CE_Failure;

    }

    CPLErr eErr = CE_None;

    for (int nYOffset = 0; nYOffset < nYSize && eErr == CE_None;

         nYOffset += nBlockYSize)

    {

        for (int nXOffset = 0; nXOffset < nXSize && eErr == CE_None;

             nXOffset += nBlockXSize)

        {

            std::unique_ptr<void, GDALScaledProgressReleaser> pScaledProgress(

                GDALCreateScaledProgress(

                    static_cast<double>(nBlock) / dfBlockCount,

                    static_cast<double>(nBlock + 1) / dfBlockCount, pfnProgress,

                    pProgressData));

            nBlock++;

            int nXRequest = nBlockXSize;

            if (nXOffset > nXSize - nXRequest)

                nXRequest = nXSize - nXOffset;

            int nYRequest = nBlockYSize;

            if (nYOffset > nYSize - nYRequest)

                nYRequest = nYSize - nYOffset;

            eErr = GDALGridContextProcess(

                psContext.get(), dfXMin + dfDeltaX * nXOffset,

                dfXMin + dfDeltaX * (nXOffset + nXRequest),

                dfYMin + dfDeltaY * nYOffset,

                dfYMin + dfDeltaY * (nYOffset + nYRequest), nXRequest,

                nYRequest, eType, pData.get(), GDALScaledProgress,

                pScaledProgress.get());

            if (eErr == CE_None)

                eErr = poBand->RasterIO(GF_Write, nXOffset, nYOffset, nXRequest,

                                        nYRequest, pData.get(), nXRequest,

                                        nYRequest, eType, 0, 0, nullptr);

        }

    }

    if (eErr == CE_None && pfnProgress)

        pfnProgress(1.0, "", pProgressData);

    return eErr;

}

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

/*                            LoadGeometry()                            */

/*                                                                      */

/*  Read geometries from the given dataset using specified filters and  */

/*  returns a collection of read geometries.                            */

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

static std::unique_ptr<OGRGeometry> LoadGeometry(const std::string &osDS,

                                                 const std::string &osSQL,

                                                 const std::string &osLyr,

                                                 const std::string &osWhere)

{

    auto poDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(

        osDS.c_str(), GDAL_OF_VECTOR, nullptr, nullptr, nullptr));

    if (!poDS)

        return nullptr;

    OGRLayer *poLyr = nullptr;

    if (!osSQL.empty())

        poLyr = poDS->ExecuteSQL(osSQL.c_str(), nullptr, nullptr);

    else if (!osLyr.empty())

        poLyr = poDS->GetLayerByName(osLyr.c_str());

    else

        poLyr = poDS->GetLayer(0);

    if (poLyr == nullptr)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "Failed to identify source layer from datasource.");

        return nullptr;

    }

    if (!osWhere.empty())

        poLyr->SetAttributeFilter(osWhere.c_str());

    std::unique_ptr<OGRGeometryCollection> poGeom;

    for (auto &poFeat : poLyr)

    {

        const OGRGeometry *poSrcGeom = poFeat->GetGeometryRef();

        if (poSrcGeom)

        {

            const OGRwkbGeometryType eType =

                wkbFlatten(poSrcGeom->getGeometryType());

            if (!poGeom)

                poGeom = std::make_unique<OGRMultiPolygon>();

            if (eType == wkbPolygon)

            {

                poGeom->addGeometry(poSrcGeom);

            }

            else if (eType == wkbMultiPolygon)

            {

                const int nGeomCount =

                    poSrcGeom->toMultiPolygon()->getNumGeometries();

                for (int iGeom = 0; iGeom < nGeomCount; iGeom++)

                {

                    poGeom->addGeometry(

                        poSrcGeom->toMultiPolygon()->getGeometryRef(iGeom));

                }

            }

            else

            {

                CPLError(CE_Failure, CPLE_AppDefined,

                         "Geometry not of polygon type.");

                if (!osSQL.empty())

                    poDS->ReleaseResultSet(poLyr);

                return nullptr;

            }

        }

    }

    if (!osSQL.empty())

        poDS->ReleaseResultSet(poLyr);

    return poGeom;

}

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

/*                              GDALGrid()                              */

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

/* clang-format off */

/**

 * Create raster from the scattered data.

 *

 * This is the equivalent of the

 * <a href="/programs/gdal_grid.html">gdal_grid</a> utility.

 *

 * GDALGridOptions* must be allocated and freed with GDALGridOptionsNew()

 * and GDALGridOptionsFree() respectively.

 *

 * @param pszDest the destination dataset path.

 * @param hSrcDataset the source dataset handle.

 * @param psOptionsIn the options struct returned by GDALGridOptionsNew() or

 * NULL.

 * @param pbUsageError pointer to a integer output variable to store if any

 * usage error has occurred or NULL.

 * @return the output dataset (new dataset that must be closed using

 * GDALClose()) or NULL in case of error.

 *

 * @since GDAL 2.1

 */

/* clang-format on */

GDALDatasetH GDALGrid(const char *pszDest, GDALDatasetH hSrcDataset,

                      const GDALGridOptions *psOptionsIn, int *pbUsageError)

{

    if (hSrcDataset == nullptr)

    {

        CPLError(CE_Failure, CPLE_AppDefined, "No source dataset specified.");

        if (pbUsageError)

            *pbUsageError = TRUE;

        return nullptr;

    }

    if (pszDest == nullptr)

    {

        CPLError(CE_Failure, CPLE_AppDefined, "No target dataset specified.");

        if (pbUsageError)

            *pbUsageError = TRUE;

        return nullptr;

    }

    std::unique_ptr<GDALGridOptions> psOptionsToFree;

    const GDALGridOptions *psOptions = psOptionsIn;

    if (psOptions == nullptr)

    {

        psOptionsToFree = std::make_unique<GDALGridOptions>();

        psOptions = psOptionsToFree.get();

    }

    GDALDataset *poSrcDS = GDALDataset::FromHandle(hSrcDataset);

    if (psOptions->osSQL.empty() && psOptions->aosLayers.empty() &&

        poSrcDS->GetLayerCount() != 1)

    {

        CPLError(CE_Failure, CPLE_NotSupported,

                 "Neither -sql nor -l are specified, but the source dataset "

                 "has not one single layer.");

        if (pbUsageError)

            *pbUsageError = TRUE;

        return nullptr;

    }

    if ((psOptions->nXSize != 0 || psOptions->nYSize != 0) &&

        (psOptions->dfXRes != 0 || psOptions->dfYRes != 0))

    {

        CPLError(CE_Failure, CPLE_IllegalArg,

                 "-outsize and -tr options cannot be used at the same time.");

        return nullptr;

    }

    /* -------------------------------------------------------------------- */

    /*      Find the output driver.                                         */

    /* -------------------------------------------------------------------- */

    std::string osFormat;

    if (psOptions->osFormat.empty())

    {

        osFormat = GetOutputDriverForRaster(pszDest);

        if (osFormat.empty())

        {

            return nullptr;

        }

    }

    else

    {

        osFormat = psOptions->osFormat;

    }

    GDALDriverH hDriver = GDALGetDriverByName(osFormat.c_str());

    if (hDriver == nullptr)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "Output driver `%s' not recognised.", osFormat.c_str());

        fprintf(stderr, "The following format drivers are enabled and "

                        "support writing:\n");

        for (int iDr = 0; iDr < GDALGetDriverCount(); iDr++)

        {

            hDriver = GDALGetDriver(iDr);

            if (GDALGetMetadataItem(hDriver, GDAL_DCAP_RASTER, nullptr) !=

                    nullptr &&

                (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, nullptr) !=

                     nullptr ||

                 GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATECOPY, nullptr) !=

                     nullptr))

            {

                fprintf(stderr, "  %s: %s\n", GDALGetDriverShortName(hDriver),

                        GDALGetDriverLongName(hDriver));

            }

        }

        printf("\n");

        return nullptr;

    }

    /* -------------------------------------------------------------------- */

    /*      Create target raster file.                                      */

    /* -------------------------------------------------------------------- */

    int nLayerCount = psOptions->aosLayers.size();

    if (nLayerCount == 0 && psOptions->osSQL.empty())

        nLayerCount = 1; /* due to above check */

    int nBands = nLayerCount;

    if (!psOptions->osSQL.empty())

        nBands++;

    int nXSize;

    int nYSize;

    if (psOptions->dfXRes != 0 && psOptions->dfYRes != 0)

    {

        double dfXSize = (std::fabs(psOptions->dfXMax - psOptions->dfXMin) +

                          (psOptions->dfXRes / 2.0)) /

                         psOptions->dfXRes;

        double dfYSize = (std::fabs(psOptions->dfYMax - psOptions->dfYMin) +

                          (psOptions->dfYRes / 2.0)) /

                         psOptions->dfYRes;

        if (dfXSize >= 1 && dfXSize <= INT_MAX && dfYSize >= 1 &&

            dfYSize <= INT_MAX)

        {

            nXSize = static_cast<int>(dfXSize);

            nYSize = static_cast<int>(dfYSize);

        }

        else

        {

            CPLError(

                CE_Failure, CPLE_IllegalArg,

                "Invalid output size detected. Please check your -tr argument");

            if (pbUsageError)

                *pbUsageError = TRUE;

            return nullptr;

        }

    }

    else

    {

        // FIXME

        nXSize = psOptions->nXSize;

        if (nXSize == 0)

            nXSize = 256;

        nYSize = psOptions->nYSize;

        if (nYSize == 0)

            nYSize = 256;

    }

    std::unique_ptr<GDALDataset> poDstDS(GDALDataset::FromHandle(GDALCreate(

        hDriver, pszDest, nXSize, nYSize, nBands, psOptions->eOutputType,

        psOptions->aosCreateOptions.List())));

    if (!poDstDS)

    {

        return nullptr;

    }

    if (psOptions->bNoDataSet)

    {

        for (int i = 1; i <= nBands; i++)

        {

            poDstDS->GetRasterBand(i)->SetNoDataValue(psOptions->dfNoDataValue);

        }

    }

    double dfXMin = psOptions->dfXMin;

    double dfYMin = psOptions->dfYMin;

    double dfXMax = psOptions->dfXMax;

    double dfYMax = psOptions->dfYMax;

    bool bIsXExtentSet = psOptions->bIsXExtentSet;

    bool bIsYExtentSet = psOptions->bIsYExtentSet;

    CPLErr eErr = CE_None;

    const bool bCloseReportsProgress = poDstDS->GetCloseReportsProgress();

    /* -------------------------------------------------------------------- */

    /*      Process SQL request.                                            */

    /* -------------------------------------------------------------------- */

    if (!psOptions->osSQL.empty())

    {

        OGRLayer *poLayer =

            poSrcDS->ExecuteSQL(psOptions->osSQL.c_str(),

                                psOptions->poSpatialFilter.get(), nullptr);

        if (poLayer == nullptr)

        {

            return nullptr;

        }

        std::unique_ptr<void, decltype(&GDALDestroyScaledProgress)>

            pScaledProgressArg(

                GDALCreateScaledProgress(0.0, bCloseReportsProgress ? 0.5 : 1.0,

                                         psOptions->pfnProgress,

                                         psOptions->pProgressData),

                GDALDestroyScaledProgress);

        // Custom layer will be rasterized in the first band.

        eErr = ProcessLayer(

            poLayer, poDstDS.get(), psOptions->poSpatialFilter.get(), nXSize,

            nYSize, 1, bIsXExtentSet, bIsYExtentSet, dfXMin, dfXMax, dfYMin,

            dfYMax, psOptions->osBurnAttribute, psOptions->dfIncreaseBurnValue,

            psOptions->dfMultiplyBurnValue, psOptions->eOutputType,

            psOptions->eAlgorithm, psOptions->pOptions.get(), psOptions->bQuiet,

            GDALScaledProgress, pScaledProgressArg.get());

        poSrcDS->ReleaseResultSet(poLayer);

    }

    /* -------------------------------------------------------------------- */

    /*      Process each layer.                                             */

    /* -------------------------------------------------------------------- */

    std::string osOutputSRS(psOptions->osOutputSRS);

    for (int i = 0; i < nLayerCount; i++)

    {

        auto poLayer = psOptions->aosLayers.empty()

                           ? poSrcDS->GetLayer(0)

                           : poSrcDS->GetLayerByName(psOptions->aosLayers[i]);

        if (!poLayer)

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "Unable to find layer \"%s\".",

                     !psOptions->aosLayers.empty() && psOptions->aosLayers[i]

                         ? psOptions->aosLayers[i]

                         : "null");

            eErr = CE_Failure;

            break;

        }

        if (!psOptions->osWHERE.empty())

        {

            if (poLayer->SetAttributeFilter(psOptions->osWHERE.c_str()) !=

                OGRERR_NONE)

            {

                eErr = CE_Failure;

                break;

            }

        }

        if (psOptions->poSpatialFilter)

            poLayer->SetSpatialFilter(psOptions->poSpatialFilter.get());

        // Fetch the first meaningful SRS definition

        if (osOutputSRS.empty())

        {

            auto poSRS = poLayer->GetSpatialRef();

            if (poSRS)

                osOutputSRS = poSRS->exportToWkt();

        }

        std::unique_ptr<void, decltype(&GDALDestroyScaledProgress)>

            pScaledProgressArg(

                GDALCreateScaledProgress(0.0, bCloseReportsProgress ? 0.5 : 1.0,

                                         psOptions->pfnProgress,

                                         psOptions->pProgressData),

                GDALDestroyScaledProgress);

        eErr = ProcessLayer(

            poLayer, poDstDS.get(), psOptions->poSpatialFilter.get(), nXSize,

            nYSize, i + 1 + nBands - nLayerCount, bIsXExtentSet, bIsYExtentSet,

            dfXMin, dfXMax, dfYMin, dfYMax, psOptions->osBurnAttribute,

            psOptions->dfIncreaseBurnValue, psOptions->dfMultiplyBurnValue,

            psOptions->eOutputType, psOptions->eAlgorithm,

            psOptions->pOptions.get(), psOptions->bQuiet, GDALScaledProgress,

            pScaledProgressArg.get());

        if (eErr != CE_None)

            break;

    }

    /* -------------------------------------------------------------------- */

    /*      Apply geotransformation matrix.                                 */