/******************************************************************************
 *
 * Project:  GDAL Utilities
 * Purpose:  Rasterize OGR shapes into a GDAL raster.
 * Author:   Frank Warmerdam <warmerdam@pobox.com>
 *
 ******************************************************************************
 * Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
 * Copyright (c) 2008-2015, Even Rouault <even dot rouault at spatialys dot com>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ****************************************************************************/

#include "cpl_port.h"
#include "gdal_utils.h"
#include "gdal_utils_priv.h"

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <limits>
#include <vector>

#include "commonutils.h"
#include "cpl_conv.h"
#include "cpl_error.h"
#include "cpl_progress.h"
#include "cpl_string.h"
#include "gdal.h"
#include "gdal_alg.h"
#include "gdal_priv.h"
#include "ogr_api.h"
#include "ogr_core.h"
#include "ogr_srs_api.h"

/************************************************************************/
/*                       ArgIsNumericRasterize()                        */
/************************************************************************/

static bool ArgIsNumericRasterize(const char *pszArg)

{
    char *pszEnd = nullptr;
    CPLStrtod(pszArg, &pszEnd);
    return pszEnd != nullptr && pszEnd[0] == '\0';
}

/************************************************************************/
/*                          InvertGeometries()                          */
/************************************************************************/

static void InvertGeometries(GDALDatasetH hDstDS,
                             std::vector<OGRGeometryH> &ahGeometries)

{
    OGRMultiPolygon *poInvertMP = new OGRMultiPolygon();

    /* -------------------------------------------------------------------- */
    /*      Create a ring that is a bit outside the raster dataset.         */
    /* -------------------------------------------------------------------- */
    const int brx = GDALGetRasterXSize(hDstDS) + 2;
    const int bry = GDALGetRasterYSize(hDstDS) + 2;

    double adfGeoTransform[6] = {};
    GDALGetGeoTransform(hDstDS, adfGeoTransform);

    OGRLinearRing *poUniverseRing = new OGRLinearRing();

    poUniverseRing->addPoint(
        adfGeoTransform[0] + -2 * adfGeoTransform[1] + -2 * adfGeoTransform[2],
        adfGeoTransform[3] + -2 * adfGeoTransform[4] + -2 * adfGeoTransform[5]);

    poUniverseRing->addPoint(adfGeoTransform[0] + brx * adfGeoTransform[1] +
                                 -2 * adfGeoTransform[2],
                             adfGeoTransform[3] + brx * adfGeoTransform[4] +
                                 -2 * adfGeoTransform[5]);

    poUniverseRing->addPoint(adfGeoTransform[0] + brx * adfGeoTransform[1] +
                                 bry * adfGeoTransform[2],
                             adfGeoTransform[3] + brx * adfGeoTransform[4] +
                                 bry * adfGeoTransform[5]);

    poUniverseRing->addPoint(adfGeoTransform[0] + -2 * adfGeoTransform[1] +
                                 bry * adfGeoTransform[2],
                             adfGeoTransform[3] + -2 * adfGeoTransform[4] +
                                 bry * adfGeoTransform[5]);

    poUniverseRing->addPoint(
        adfGeoTransform[0] + -2 * adfGeoTransform[1] + -2 * adfGeoTransform[2],
        adfGeoTransform[3] + -2 * adfGeoTransform[4] + -2 * adfGeoTransform[5]);

    OGRPolygon *poUniversePoly = new OGRPolygon();
    poUniversePoly->addRingDirectly(poUniverseRing);
    poInvertMP->addGeometryDirectly(poUniversePoly);

    bool bFoundNonPoly = false;
    // If we have GEOS, use it to "subtract" each polygon from the universe
    // multipolygon
    if (OGRGeometryFactory::haveGEOS())
    {
        OGRGeometry *poInvertMPAsGeom = poInvertMP;
        poInvertMP = nullptr;
        CPL_IGNORE_RET_VAL(poInvertMP);
        for (unsigned int iGeom = 0; iGeom < ahGeometries.size(); iGeom++)
        {
            auto poGeom = OGRGeometry::FromHandle(ahGeometries[iGeom]);
            const auto eGType = OGR_GT_Flatten(poGeom->getGeometryType());
            if (eGType != wkbPolygon && eGType != wkbMultiPolygon)
            {
                if (!bFoundNonPoly)
                {
                    bFoundNonPoly = true;
                    CPLError(CE_Warning, CPLE_AppDefined,
                             "Ignoring non-polygon geometries in -i mode");
                }
            }
            else
            {
                auto poNewGeom = poInvertMPAsGeom->Difference(poGeom);
                if (poNewGeom)
                {
                    delete poInvertMPAsGeom;
                    poInvertMPAsGeom = poNewGeom;
                }
            }

            delete poGeom;
        }

        ahGeometries.resize(1);
        ahGeometries[0] = OGRGeometry::ToHandle(poInvertMPAsGeom);
        return;
    }

    /* -------------------------------------------------------------------- */
    /*      If we don't have GEOS, add outer rings of polygons as inner     */
    /*      rings of poUniversePoly and inner rings as sub-polygons. Note   */
    /*      that this only works properly if the polygons are disjoint, in  */
    /*      the sense that the outer ring of any polygon is not inside the  */
    /*      outer ring of another one. So the scenario of                   */
    /*      https://github.com/OSGeo/gdal/issues/8689 with an "island" in   */
    /*      the middle of a hole will not work properly.                    */
    /* -------------------------------------------------------------------- */
    for (unsigned int iGeom = 0; iGeom < ahGeometries.size(); iGeom++)
    {
        const auto eGType =
            OGR_GT_Flatten(OGR_G_GetGeometryType(ahGeometries[iGeom]));
        if (eGType != wkbPolygon && eGType != wkbMultiPolygon)
        {
            if (!bFoundNonPoly)
            {
                bFoundNonPoly = true;
                CPLError(CE_Warning, CPLE_AppDefined,
                         "Ignoring non-polygon geometries in -i mode");
            }
            OGR_G_DestroyGeometry(ahGeometries[iGeom]);
            continue;
        }

        const auto ProcessPoly =
            [poUniversePoly, poInvertMP](OGRPolygon *poPoly)
        {
            for (int i = poPoly->getNumInteriorRings() - 1; i >= 0; --i)
            {
                auto poNewPoly = new OGRPolygon();
                poNewPoly->addRingDirectly(poPoly->stealInteriorRing(i));
                poInvertMP->addGeometryDirectly(poNewPoly);
            }
            poUniversePoly->addRingDirectly(poPoly->stealExteriorRing());
        };

        if (eGType == wkbPolygon)
        {
            auto poPoly =
                OGRGeometry::FromHandle(ahGeometries[iGeom])->toPolygon();
            ProcessPoly(poPoly);
            delete poPoly;
        }
        else
        {
            auto poMulti =
                OGRGeometry::FromHandle(ahGeometries[iGeom])->toMultiPolygon();
            for (auto *poPoly : *poMulti)
            {
                ProcessPoly(poPoly);
            }
            delete poMulti;
        }
    }

    ahGeometries.resize(1);
    ahGeometries[0] = OGRGeometry::ToHandle(poInvertMP);
}

/************************************************************************/
/*                            ProcessLayer()                            */
/*                                                                      */
/*      Process all the features in a layer selection, collecting       */
/*      geometries and burn values.                                     */
/************************************************************************/

static CPLErr ProcessLayer(OGRLayerH hSrcLayer, bool bSRSIsSet,
                           GDALDatasetH hDstDS,
                           const std::vector<int> &anBandList,
                           const std::vector<double> &adfBurnValues, bool b3D,
                           bool bInverse, const char *pszBurnAttribute,
                           CSLConstList papszRasterizeOptions, char **papszTO,
                           GDALProgressFunc pfnProgress, void *pProgressData)

{
    /* -------------------------------------------------------------------- */
    /*      Checkout that SRS are the same.                                 */
    /*      If -a_srs is specified, skip the test                           */
    /* -------------------------------------------------------------------- */
    OGRCoordinateTransformationH hCT = nullptr;
    if (!bSRSIsSet)
    {
        OGRSpatialReferenceH hDstSRS = GDALGetSpatialRef(hDstDS);

        if (hDstSRS)
            hDstSRS = OSRClone(hDstSRS);
        else if (GDALGetMetadata(hDstDS, "RPC") != nullptr)
        {
            hDstSRS = OSRNewSpatialReference(nullptr);
            CPL_IGNORE_RET_VAL(
                OSRSetFromUserInput(hDstSRS, SRS_WKT_WGS84_LAT_LONG));
            OSRSetAxisMappingStrategy(hDstSRS, OAMS_TRADITIONAL_GIS_ORDER);
        }

        OGRSpatialReferenceH hSrcSRS = OGR_L_GetSpatialRef(hSrcLayer);
        if (hDstSRS != nullptr && hSrcSRS != nullptr)
        {
            if (OSRIsSame(hSrcSRS, hDstSRS) == FALSE)
            {
                hCT = OCTNewCoordinateTransformation(hSrcSRS, hDstSRS);
                if (hCT == nullptr)
                {
                    CPLError(CE_Warning, CPLE_AppDefined,
                             "The output raster dataset and the input vector "
                             "layer do not have the same SRS.\n"
                             "And reprojection of input data did not work. "
                             "Results might be incorrect.");
                }
            }
        }
        else if (hDstSRS != nullptr && hSrcSRS == nullptr)
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                     "The output raster dataset has a SRS, but the input "
                     "vector layer SRS is unknown.\n"
                     "Ensure input vector has the same SRS, otherwise results "
                     "might be incorrect.");
        }
        else if (hDstSRS == nullptr && hSrcSRS != nullptr)
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                     "The input vector layer has a SRS, but the output raster "
                     "dataset SRS is unknown.\n"
                     "Ensure output raster dataset has the same SRS, otherwise "
                     "results might be incorrect.");
        }

        if (hDstSRS != nullptr)
        {
            OSRDestroySpatialReference(hDstSRS);
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Get field index, and check.                                     */
    /* -------------------------------------------------------------------- */
    int iBurnField = -1;
    bool bUseInt64 = false;
    if (pszBurnAttribute)
    {
        OGRFeatureDefnH hLayerDefn = OGR_L_GetLayerDefn(hSrcLayer);
        iBurnField = OGR_FD_GetFieldIndex(hLayerDefn, pszBurnAttribute);
        if (iBurnField == -1)
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Failed to find field %s on layer %s.", pszBurnAttribute,
                     OGR_FD_GetName(OGR_L_GetLayerDefn(hSrcLayer)));
            if (hCT != nullptr)
                OCTDestroyCoordinateTransformation(hCT);
            return CE_Failure;
        }
        if (OGR_Fld_GetType(OGR_FD_GetFieldDefn(hLayerDefn, iBurnField)) ==
            OFTInteger64)
        {
            GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, anBandList[0]);
            if (hBand && GDALGetRasterDataType(hBand) == GDT_Int64)
            {
                bUseInt64 = true;
            }
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Collect the geometries from this layer, and build list of       */
    /*      burn values.                                                    */
    /* -------------------------------------------------------------------- */
    OGRFeatureH hFeat = nullptr;
    std::vector<OGRGeometryH> ahGeometries;
    std::vector<double> adfFullBurnValues;
    std::vector<int64_t> anFullBurnValues;

    OGR_L_ResetReading(hSrcLayer);

    while ((hFeat = OGR_L_GetNextFeature(hSrcLayer)) != nullptr)
    {
        OGRGeometryH hGeom = OGR_F_StealGeometry(hFeat);
        if (hGeom == nullptr)
        {
            OGR_F_Destroy(hFeat);
            continue;
        }

        if (hCT != nullptr)
        {
            if (OGR_G_Transform(hGeom, hCT) != OGRERR_NONE)
            {
                OGR_F_Destroy(hFeat);
                OGR_G_DestroyGeometry(hGeom);
                continue;
            }
        }
        ahGeometries.push_back(hGeom);

        for (unsigned int iBand = 0; iBand < anBandList.size(); iBand++)
        {
            if (!adfBurnValues.empty())
                adfFullBurnValues.push_back(adfBurnValues[std::min(
                    iBand,
                    static_cast<unsigned int>(adfBurnValues.size()) - 1)]);
            else if (pszBurnAttribute)
            {
                if (bUseInt64)
                    anFullBurnValues.push_back(
                        OGR_F_GetFieldAsInteger64(hFeat, iBurnField));
                else
                    adfFullBurnValues.push_back(
                        OGR_F_GetFieldAsDouble(hFeat, iBurnField));
            }
            else if (b3D)
            {
                /* Points and Lines will have their "z" values collected at the
                   point and line levels respectively. Not implemented for
                   polygons */
                adfFullBurnValues.push_back(0.0);
            }
        }

        OGR_F_Destroy(hFeat);
    }

    if (hCT != nullptr)
        OCTDestroyCoordinateTransformation(hCT);

    /* -------------------------------------------------------------------- */
    /*      If we are in inverse mode, we add one extra ring around the     */
    /*      whole dataset to invert the concept of insideness and then      */
    /*      merge everything into one geometry collection.                  */
    /* -------------------------------------------------------------------- */
    if (bInverse)
    {
        if (ahGeometries.empty())
        {
            for (unsigned int iBand = 0; iBand < anBandList.size(); iBand++)
            {
                if (!adfBurnValues.empty())
                    adfFullBurnValues.push_back(adfBurnValues[std::min(
                        iBand,
                        static_cast<unsigned int>(adfBurnValues.size()) - 1)]);
                else /* FIXME? Not sure what to do exactly in the else case, but
                        we must insert a value */
                {
                    adfFullBurnValues.push_back(0.0);
                    anFullBurnValues.push_back(0);
                }
            }
        }

        InvertGeometries(hDstDS, ahGeometries);
    }

    /* -------------------------------------------------------------------- */
    /*      If we have transformer options, create the transformer here     */
    /*      Coordinate transformation to the target SRS has already been    */
    /*      done, so we just need to convert to target raster space.        */
    /*      Note: this is somewhat identical to what is done in             */
    /*      GDALRasterizeGeometries() itself, except we can pass transformer*/
    /*      options.                                                        */
    /* -------------------------------------------------------------------- */

    void *pTransformArg = nullptr;
    GDALTransformerFunc pfnTransformer = nullptr;
    CPLErr eErr = CE_None;
    if (papszTO != nullptr)
    {
        GDALDataset *poDS = GDALDataset::FromHandle(hDstDS);
        char **papszTransformerOptions = CSLDuplicate(papszTO);
        double adfGeoTransform[6] = {0.0};
        if (poDS->GetGeoTransform(adfGeoTransform) != CE_None &&
            poDS->GetGCPCount() == 0 && poDS->GetMetadata("RPC") == nullptr)
        {
            papszTransformerOptions = CSLSetNameValue(
                papszTransformerOptions, "DST_METHOD", "NO_GEOTRANSFORM");
        }

        pTransformArg = GDALCreateGenImgProjTransformer2(
            nullptr, hDstDS, papszTransformerOptions);
        CSLDestroy(papszTransformerOptions);

        pfnTransformer = GDALGenImgProjTransform;
        if (pTransformArg == nullptr)
        {
            eErr = CE_Failure;
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Perform the burn.                                               */
    /* -------------------------------------------------------------------- */
    if (eErr == CE_None)
    {
        if (bUseInt64)
        {
            eErr = GDALRasterizeGeometriesInt64(
                hDstDS, static_cast<int>(anBandList.size()), anBandList.data(),
                static_cast<int>(ahGeometries.size()), ahGeometries.data(),
                pfnTransformer, pTransformArg, anFullBurnValues.data(),
                papszRasterizeOptions, pfnProgress, pProgressData);
        }
        else
        {
            eErr = GDALRasterizeGeometries(
                hDstDS, static_cast<int>(anBandList.size()), anBandList.data(),
                static_cast<int>(ahGeometries.size()), ahGeometries.data(),
                pfnTransformer, pTransformArg, adfFullBurnValues.data(),
                papszRasterizeOptions, pfnProgress, pProgressData);
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Cleanup                                                         */
    /* -------------------------------------------------------------------- */

    if (pTransformArg)
        GDALDestroyTransformer(pTransformArg);

    for (int iGeom = static_cast<int>(ahGeometries.size()) - 1; iGeom >= 0;
         iGeom--)
        OGR_G_DestroyGeometry(ahGeometries[iGeom]);

    return eErr;
}

/************************************************************************/
/*                  CreateOutputDataset()                               */
/************************************************************************/

static GDALDatasetH CreateOutputDataset(
    const std::vector<OGRLayerH> &ahLayers, OGRSpatialReferenceH hSRS,
    OGREnvelope sEnvelop, GDALDriverH hDriver, const char *pszDest, int nXSize,
    int nYSize, double dfXRes, double dfYRes, bool bTargetAlignedPixels,
    int nBandCount, GDALDataType eOutputType, char **papszCreationOptions,
    const std::vector<double> &adfInitVals, const char *pszNoData)
{
    bool bFirstLayer = true;
    char *pszWKT = nullptr;
    const bool bBoundsSpecifiedByUser = sEnvelop.IsInit();

    for (unsigned int i = 0; i < ahLayers.size(); i++)
    {
        OGRLayerH hLayer = ahLayers[i];

        if (!bBoundsSpecifiedByUser)
        {
            OGREnvelope sLayerEnvelop;

            if (OGR_L_GetExtent(hLayer, &sLayerEnvelop, TRUE) != OGRERR_NONE)
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Cannot get layer extent");
                return nullptr;
            }

            /* Voluntarily increase the extent by a half-pixel size to avoid */
            /* missing points on the border */
            if (!bTargetAlignedPixels && dfXRes != 0 && dfYRes != 0)
            {
                sLayerEnvelop.MinX -= dfXRes / 2;
                sLayerEnvelop.MaxX += dfXRes / 2;
                sLayerEnvelop.MinY -= dfYRes / 2;
                sLayerEnvelop.MaxY += dfYRes / 2;
            }

            sEnvelop.Merge(sLayerEnvelop);
        }

        if (bFirstLayer)
        {
            if (hSRS == nullptr)
                hSRS = OGR_L_GetSpatialRef(hLayer);

            bFirstLayer = false;
        }
    }

    if (!sEnvelop.IsInit())
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Could not determine bounds");
        return nullptr;
    }

    if (dfXRes == 0 && dfYRes == 0)
    {
        if (nXSize == 0 || nYSize == 0)
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Size and resolutions are missing");
            return nullptr;
        }
        dfXRes = (sEnvelop.MaxX - sEnvelop.MinX) / nXSize;
        dfYRes = (sEnvelop.MaxY - sEnvelop.MinY) / nYSize;
    }
    else if (bTargetAlignedPixels && dfXRes != 0 && dfYRes != 0)
    {
        sEnvelop.MinX = floor(sEnvelop.MinX / dfXRes) * dfXRes;
        sEnvelop.MaxX = ceil(sEnvelop.MaxX / dfXRes) * dfXRes;
        sEnvelop.MinY = floor(sEnvelop.MinY / dfYRes) * dfYRes;
        sEnvelop.MaxY = ceil(sEnvelop.MaxY / dfYRes) * dfYRes;
    }

    double adfProjection[6] = {sEnvelop.MinX, dfXRes, 0.0,
                               sEnvelop.MaxY, 0.0,    -dfYRes};

    if (nXSize == 0 && nYSize == 0)
    {
        const double dfXSize = 0.5 + (sEnvelop.MaxX - sEnvelop.MinX) / dfXRes;
        const double dfYSize = 0.5 + (sEnvelop.MaxY - sEnvelop.MinY) / dfYRes;
        if (dfXSize > std::numeric_limits<int>::max() ||
            dfXSize < std::numeric_limits<int>::min() ||
            dfYSize > std::numeric_limits<int>::max() ||
            dfYSize < std::numeric_limits<int>::min())
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Invalid computed output raster size: %f x %f", dfXSize,
                     dfYSize);
            return nullptr;
        }
        nXSize = static_cast<int>(dfXSize);
        nYSize = static_cast<int>(dfYSize);
    }

    GDALDatasetH hDstDS =
        GDALCreate(hDriver, pszDest, nXSize, nYSize, nBandCount, eOutputType,
                   papszCreationOptions);
    if (hDstDS == nullptr)
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Cannot create %s", pszDest);
        return nullptr;
    }

    GDALSetGeoTransform(hDstDS, adfProjection);

    if (hSRS)
        OSRExportToWkt(hSRS, &pszWKT);
    if (pszWKT)
        GDALSetProjection(hDstDS, pszWKT);
    CPLFree(pszWKT);

    /*if( nBandCount == 3 || nBandCount == 4 )
    {
        for( int iBand = 0; iBand < nBandCount; iBand++ )
        {
            GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);
            GDALSetRasterColorInterpretation(hBand,
    (GDALColorInterp)(GCI_RedBand + iBand));
        }
    }*/

    if (pszNoData)
    {
        for (int iBand = 0; iBand < nBandCount; iBand++)
        {
            GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);
            if (GDALGetRasterDataType(hBand) == GDT_Int64)
                GDALSetRasterNoDataValueAsInt64(hBand,
                                                CPLAtoGIntBig(pszNoData));
            else
                GDALSetRasterNoDataValue(hBand, CPLAtof(pszNoData));
        }
    }

    if (!adfInitVals.empty())
    {
        for (int iBand = 0;
             iBand < std::min(nBandCount, static_cast<int>(adfInitVals.size()));
             iBand++)
        {
            GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);
            GDALFillRaster(hBand, adfInitVals[iBand], 0);
        }
    }

    return hDstDS;
}

struct GDALRasterizeOptions
{
    /*! output format. Use the short format name. */
    char *pszFormat;

    /*! the progress function to use */
    GDALProgressFunc pfnProgress;

    /*! pointer to the progress data variable */
    void *pProgressData;

    bool bCreateOutput;
    bool b3D;
    bool bInverse;
    char **papszLayers;
    char *pszSQL;
    char *pszDialect;
    char *pszBurnAttribute;
    char *pszWHERE;
    std::vector<int> anBandList;
    std::vector<double> adfBurnValues;
    char **papszRasterizeOptions;
    char **papszTO;
    double dfXRes;
    double dfYRes;
    char **papszCreationOptions;
    GDALDataType eOutputType;
    std::vector<double> adfInitVals;
    char *pszNoData;
    OGREnvelope sEnvelop;
    int nXSize, nYSize;
    OGRSpatialReferenceH hSRS;
    bool bTargetAlignedPixels;
};

/************************************************************************/
/*                             GDALRasterize()                          */
/************************************************************************/

/* clang-format off */
/**
 * Burns vector geometries into a raster
 *
 * This is the equivalent of the
 * <a href="/programs/gdal_rasterize.html">gdal_rasterize</a> utility.
 *
 * GDALRasterizeOptions* must be allocated and freed with
 * GDALRasterizeOptionsNew() and GDALRasterizeOptionsFree() respectively.
 * pszDest and hDstDS cannot be used at the same time.
 *
 * @param pszDest the destination dataset path or NULL.
 * @param hDstDS the destination dataset or NULL.
 * @param hSrcDataset the source dataset handle.
 * @param psOptionsIn the options struct returned by GDALRasterizeOptionsNew()
 * 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 hDstDS is not NULL) or NULL in case of error.
 *
 * @since GDAL 2.1
 */
/* clang-format on */

GDALDatasetH GDALRasterize(const char *pszDest, GDALDatasetH hDstDS,
                           GDALDatasetH hSrcDataset,
                           const GDALRasterizeOptions *psOptionsIn,
                           int *pbUsageError)
{
    if (pszDest == nullptr && hDstDS == nullptr)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "pszDest == NULL && hDstDS == NULL");

        if (pbUsageError)
            *pbUsageError = TRUE;
        return nullptr;
    }
    if (hSrcDataset == nullptr)
    {
        CPLError(CE_Failure, CPLE_AppDefined, "hSrcDataset== NULL");

        if (pbUsageError)
            *pbUsageError = TRUE;
        return nullptr;
    }
    if (hDstDS != nullptr && psOptionsIn && psOptionsIn->bCreateOutput)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "hDstDS != NULL but options that imply creating a new dataset "
                 "have been set.");

        if (pbUsageError)
            *pbUsageError = TRUE;
        return nullptr;
    }

    GDALRasterizeOptions *psOptionsToFree = nullptr;
    const GDALRasterizeOptions *psOptions = psOptionsIn;
    if (psOptions == nullptr)
    {
        psOptionsToFree = GDALRasterizeOptionsNew(nullptr, nullptr);
        psOptions = psOptionsToFree;
    }

    const bool bCloseOutDSOnError = hDstDS == nullptr;
    if (pszDest == nullptr)
        pszDest = GDALGetDescription(hDstDS);

    if (psOptions->pszSQL == nullptr && psOptions->papszLayers == nullptr &&
        GDALDatasetGetLayerCount(hSrcDataset) != 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;
        GDALRasterizeOptionsFree(psOptionsToFree);
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*      Open target raster file.  Eventually we will add optional       */
    /*      creation.                                                       */
    /* -------------------------------------------------------------------- */
    const bool bCreateOutput = psOptions->bCreateOutput || hDstDS == nullptr;

    GDALDriverH hDriver = nullptr;
    if (bCreateOutput)
    {
        CPLString osFormat;
        if (psOptions->pszFormat == nullptr)
        {
            osFormat = GetOutputDriverForRaster(pszDest);
            if (osFormat.empty())
            {
                GDALRasterizeOptionsFree(psOptionsToFree);
                return nullptr;
            }
        }
        else
        {
            osFormat = psOptions->pszFormat;
        }

        /* --------------------------------------------------------------------
         */
        /*      Find the output driver. */
        /* --------------------------------------------------------------------
         */
        hDriver = GDALGetDriverByName(osFormat);
        char **papszDriverMD =
            hDriver ? GDALGetMetadata(hDriver, nullptr) : nullptr;
        if (hDriver == nullptr ||
            !CPLTestBool(CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_RASTER,
                                              "FALSE")) ||
            !CPLTestBool(
                CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_CREATE, "FALSE")))
        {
            CPLError(CE_Failure, CPLE_NotSupported,
                     "Output driver `%s' not recognised or does not support "
                     "direct output file creation.",
                     osFormat.c_str());
            GDALRasterizeOptionsFree(psOptionsToFree);
            return nullptr;
        }
    }

    const auto GetOutputDataType = [&](OGRLayerH hLayer)
    {
        CPLAssert(bCreateOutput);
        CPLAssert(hDriver);
        GDALDataType eOutputType = psOptions->eOutputType;
        if (eOutputType == GDT_Unknown &&
            psOptions->pszBurnAttribute != nullptr)
        {
            OGRFeatureDefnH hLayerDefn = OGR_L_GetLayerDefn(hLayer);
            const int iBurnField =
                OGR_FD_GetFieldIndex(hLayerDefn, psOptions->pszBurnAttribute);
            if (iBurnField >= 0 && OGR_Fld_GetType(OGR_FD_GetFieldDefn(
                                       hLayerDefn, iBurnField)) == OFTInteger64)
            {
                const char *pszMD = GDALGetMetadataItem(
                    hDriver, GDAL_DMD_CREATIONDATATYPES, nullptr);
                if (pszMD && CPLStringList(CSLTokenizeString2(pszMD, " ", 0))
                                     .FindString("Int64") >= 0)
                {
                    eOutputType = GDT_Int64;
                }
            }
        }
        if (eOutputType == GDT_Unknown)
        {
            eOutputType = GDT_Float64;
        }
        return eOutputType;
    };

    /* -------------------------------------------------------------------- */
    /*      Process SQL request.                                            */
    /* -------------------------------------------------------------------- */
    CPLErr eErr = CE_Failure;

    if (psOptions->pszSQL != nullptr)
    {
        OGRLayerH hLayer = GDALDatasetExecuteSQL(
            hSrcDataset, psOptions->pszSQL, nullptr, psOptions->pszDialect);
        if (hLayer != nullptr)
        {
            if (bCreateOutput)
            {
                std::vector<OGRLayerH> ahLayers;
                ahLayers.push_back(hLayer);

                const GDALDataType eOutputType = GetOutputDataType(hLayer);
                hDstDS = CreateOutputDataset(
                    ahLayers, psOptions->hSRS, psOptions->sEnvelop, hDriver,
                    pszDest, psOptions->nXSize, psOptions->nYSize,
                    psOptions->dfXRes, psOptions->dfYRes,
                    psOptions->bTargetAlignedPixels,
                    static_cast<int>(psOptions->anBandList.size()), eOutputType,
                    psOptions->papszCreationOptions, psOptions->adfInitVals,
                    psOptions->pszNoData);
                if (hDstDS == nullptr)
                {
                    GDALDatasetReleaseResultSet(hSrcDataset, hLayer);
                    GDALRasterizeOptionsFree(psOptionsToFree);
                    return nullptr;
                }
            }

            eErr = ProcessLayer(
                hLayer, psOptions->hSRS != nullptr, hDstDS,
                psOptions->anBandList, psOptions->adfBurnValues, psOptions->b3D,
                psOptions->bInverse, psOptions->pszBurnAttribute,
                psOptions->papszRasterizeOptions, psOptions->papszTO,
                psOptions->pfnProgress, psOptions->pProgressData);

            GDALDatasetReleaseResultSet(hSrcDataset, hLayer);
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Create output file if necessary.                                */
    /* -------------------------------------------------------------------- */
    const int nLayerCount =
        (psOptions->pszSQL == nullptr && psOptions->papszLayers == nullptr)
            ? 1
            : CSLCount(psOptions->papszLayers);

    if (bCreateOutput && hDstDS == nullptr)
    {
        std::vector<OGRLayerH> ahLayers;

        GDALDataType eOutputType = psOptions->eOutputType;

        for (int i = 0; i < nLayerCount; i++)
        {
            OGRLayerH hLayer;
            if (psOptions->papszLayers)
                hLayer = GDALDatasetGetLayerByName(hSrcDataset,
                                                   psOptions->papszLayers[i]);
            else
                hLayer = GDALDatasetGetLayer(hSrcDataset, 0);
            if (hLayer == nullptr)
            {
                CPLError(
                    CE_Failure, CPLE_AppDefined, "Unable to find layer \"%s\".",
                    psOptions->papszLayers ? psOptions->papszLayers[i] : "0");
                GDALRasterizeOptionsFree(psOptionsToFree);
                return nullptr;
            }
            if (eOutputType == GDT_Unknown)
            {
                if (GetOutputDataType(hLayer) == GDT_Int64)
                    eOutputType = GDT_Int64;
            }

            ahLayers.push_back(hLayer);
        }

        if (eOutputType == GDT_Unknown)
        {
            eOutputType = GDT_Float64;
        }

        hDstDS = CreateOutputDataset(
            ahLayers, psOptions->hSRS, psOptions->sEnvelop, hDriver, pszDest,
            psOptions->nXSize, psOptions->nYSize, psOptions->dfXRes,
            psOptions->dfYRes, psOptions->bTargetAlignedPixels,
            static_cast<int>(psOptions->anBandList.size()), eOutputType,
            psOptions->papszCreationOptions, psOptions->adfInitVals,
            psOptions->pszNoData);
        if (hDstDS == nullptr)
        {
            GDALRasterizeOptionsFree(psOptionsToFree);
            return nullptr;
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Process each layer.                                             */
    /* -------------------------------------------------------------------- */

    for (int i = 0; i < nLayerCount; i++)
    {
        OGRLayerH hLayer;
        if (psOptions->papszLayers)
            hLayer = GDALDatasetGetLayerByName(hSrcDataset,
                                               psOptions->papszLayers[i]);
        else
            hLayer = GDALDatasetGetLayer(hSrcDataset, 0);
        if (hLayer == nullptr)
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Unable to find layer \"%s\".",
                     psOptions->papszLayers ? psOptions->papszLayers[i] : "0");
            eErr = CE_Failure;
            break;
        }

        if (psOptions->pszWHERE)
        {
            if (OGR_L_SetAttributeFilter(hLayer, psOptions->pszWHERE) !=
                OGRERR_NONE)
            {
                eErr = CE_Failure;
                break;
            }
        }

        void *pScaledProgress = GDALCreateScaledProgress(
            0.0, 1.0 * (i + 1) / nLayerCount, psOptions->pfnProgress,
            psOptions->pProgressData);

        eErr = ProcessLayer(
            hLayer, psOptions->hSRS != nullptr, hDstDS, psOptions->anBandList,
            psOptions->adfBurnValues, psOptions->b3D, psOptions->bInverse,
            psOptions->pszBurnAttribute, psOptions->papszRasterizeOptions,
            psOptions->papszTO, GDALScaledProgress, pScaledProgress);

        GDALDestroyScaledProgress(pScaledProgress);
        if (eErr != CE_None)
            break;
    }

    GDALRasterizeOptionsFree(psOptionsToFree);

    if (eErr != CE_None)
    {
        if (bCloseOutDSOnError)
            GDALClose(hDstDS);
        return nullptr;
    }

    return hDstDS;
}

/************************************************************************/
/*                           GDALRasterizeOptionsNew()                  */
/************************************************************************/

/**
 * Allocates a GDALRasterizeOptions struct.
 *
 * @param papszArgv NULL terminated list of options (potentially including
 * filename and open options too), or NULL. The accepted options are the ones of
 * the <a href="/programs/gdal_rasterize.html">gdal_rasterize</a> utility.
 * @param psOptionsForBinary (output) may be NULL (and should generally be
 * NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with
 *                           GDALRasterizeOptionsForBinaryNew() prior to this
 * function. Will be filled with potentially present filename, open options,...
 * @return pointer to the allocated GDALRasterizeOptions struct. Must be freed
 * with GDALRasterizeOptionsFree().
 *
 * @since GDAL 2.1
 */

GDALRasterizeOptions *
GDALRasterizeOptionsNew(char **papszArgv,
                        GDALRasterizeOptionsForBinary *psOptionsForBinary)
{
    GDALRasterizeOptions *psOptions = new GDALRasterizeOptions;

    psOptions->pszFormat = nullptr;
    psOptions->pfnProgress = GDALDummyProgress;
    psOptions->pProgressData = nullptr;
    psOptions->bCreateOutput = false;
    psOptions->b3D = false;
    psOptions->bInverse = false;
    // sEnvelop implicitly initialized
    psOptions->papszCreationOptions = nullptr;
    psOptions->papszLayers = nullptr;
    psOptions->pszSQL = nullptr;
    psOptions->pszDialect = nullptr;
    psOptions->pszBurnAttribute = nullptr;
    psOptions->pszWHERE = nullptr;
    psOptions->papszRasterizeOptions = nullptr;
    psOptions->papszTO = nullptr;
    psOptions->dfXRes = 0;
    psOptions->dfYRes = 0;
    psOptions->eOutputType = GDT_Unknown;
    psOptions->pszNoData = nullptr;
    psOptions->nXSize = 0;
    psOptions->nYSize = 0;
    psOptions->hSRS = nullptr;
    psOptions->bTargetAlignedPixels = false;

    bool bGotSourceFilename = false;
    bool bGotDestFilename = false;
    /* -------------------------------------------------------------------- */
    /*      Handle command line arguments.                                  */
    /* -------------------------------------------------------------------- */
    const int argc = CSLCount(papszArgv);
    for (int i = 0; papszArgv != nullptr && i < argc; i++)
    {
        if (i < argc - 1 &&
            (EQUAL(papszArgv[i], "-of") || EQUAL(papszArgv[i], "-f")))
        {
            ++i;
            CPLFree(psOptions->pszFormat);
            psOptions->pszFormat = CPLStrdup(papszArgv[i]);
            psOptions->bCreateOutput = true;
        }

        else if (EQUAL(papszArgv[i], "-q") || EQUAL(papszArgv[i], "-quiet"))
        {
            if (psOptionsForBinary)
            {
                psOptionsForBinary->bQuiet = true;
            }
            else
            {
                CPLError(CE_Failure, CPLE_NotSupported,
                         "%s switch only supported from gdal_rasterize binary.",
                         papszArgv[i]);
            }
        }

        else if (i < argc - 1 && EQUAL(papszArgv[i], "-a"))
        {
            CPLFree(psOptions->pszBurnAttribute);
            psOptions->pszBurnAttribute = CPLStrdup(papszArgv[++i]);
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-b"))
        {
            if (strchr(papszArgv[i + 1], ' '))
            {
                const CPLStringList aosTokens(
                    CSLTokenizeString(papszArgv[i + 1]));
                for (const char *pszToken : aosTokens)
                {
                    psOptions->anBandList.push_back(atoi(pszToken));
                }
                i += 1;
            }
            else
            {
                while (i < argc - 1 && ArgIsNumericRasterize(papszArgv[i + 1]))
                {
                    psOptions->anBandList.push_back(atoi(papszArgv[i + 1]));
                    i += 1;
                }
            }
        }
        else if (EQUAL(papszArgv[i], "-3d"))
        {
            psOptions->b3D = true;
            psOptions->papszRasterizeOptions = CSLSetNameValue(
                psOptions->papszRasterizeOptions, "BURN_VALUE_FROM", "Z");
        }
        else if (EQUAL(papszArgv[i], "-add"))
        {
            psOptions->papszRasterizeOptions = CSLSetNameValue(
                psOptions->papszRasterizeOptions, "MERGE_ALG", "ADD");
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-chunkysize"))
        {
            psOptions->papszRasterizeOptions = CSLSetNameValue(
                psOptions->papszRasterizeOptions, "CHUNKYSIZE", papszArgv[++i]);
        }
        else if (EQUAL(papszArgv[i], "-i"))
        {
            psOptions->bInverse = true;
        }
        else if (EQUAL(papszArgv[i], "-at"))
        {
            psOptions->papszRasterizeOptions = CSLSetNameValue(
                psOptions->papszRasterizeOptions, "ALL_TOUCHED", "TRUE");
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-optim"))
        {
            psOptions->papszRasterizeOptions = CSLSetNameValue(
                psOptions->papszRasterizeOptions, "OPTIM", papszArgv[++i]);
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-burn"))
        {
            if (strchr(papszArgv[i + 1], ' '))
            {
                const CPLStringList aosTokens(
                    CSLTokenizeString(papszArgv[i + 1]));
                for (const char *pszToken : aosTokens)
                {
                    psOptions->adfBurnValues.push_back(CPLAtof(pszToken));
                }
                i += 1;
            }
            else
            {
                while (i < argc - 1 && ArgIsNumericRasterize(papszArgv[i + 1]))
                {
                    psOptions->adfBurnValues.push_back(
                        CPLAtof(papszArgv[i + 1]));
                    i += 1;
                }
            }
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-where"))
        {
            CPLFree(psOptions->pszWHERE);
            psOptions->pszWHERE = CPLStrdup(papszArgv[++i]);
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-l"))
        {
            psOptions->papszLayers =
                CSLAddString(psOptions->papszLayers, papszArgv[++i]);
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-sql"))
        {
            ++i;
            CPLFree(psOptions->pszSQL);
            GByte *pabyRet = nullptr;
            if (papszArgv[i][0] == '@' &&
                VSIIngestFile(nullptr, papszArgv[i] + 1, &pabyRet, nullptr,
                              1024 * 1024))
            {
                GDALRemoveBOM(pabyRet);
                char *pszSQLStatement = reinterpret_cast<char *>(pabyRet);
                psOptions->pszSQL =
                    CPLStrdup(GDALRemoveSQLComments(pszSQLStatement).c_str());
                VSIFree(pszSQLStatement);
            }
            else
            {
                psOptions->pszSQL = CPLStrdup(papszArgv[i]);
            }
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-dialect"))
        {
            CPLFree(psOptions->pszDialect);
            psOptions->pszDialect = CPLStrdup(papszArgv[++i]);
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-init"))
        {
            if (strchr(papszArgv[i + 1], ' '))
            {
                const CPLStringList aosTokens(
                    CSLTokenizeString(papszArgv[i + 1]));
                for (const char *pszToken : aosTokens)
                {
                    psOptions->adfInitVals.push_back(CPLAtof(pszToken));
                }
                i += 1;
            }
            else
            {
                while (i < argc - 1 && ArgIsNumericRasterize(papszArgv[i + 1]))
                {
                    psOptions->adfInitVals.push_back(CPLAtof(papszArgv[i + 1]));
                    i += 1;
                }
            }
            psOptions->bCreateOutput = true;
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-a_nodata"))
        {
            CPLFree(psOptions->pszNoData);
            psOptions->pszNoData = CPLStrdup(papszArgv[i + 1]);
            i += 1;
            psOptions->bCreateOutput = true;
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-a_srs"))
        {
            OSRDestroySpatialReference(psOptions->hSRS);
            psOptions->hSRS = OSRNewSpatialReference(nullptr);

            if (OSRSetFromUserInput(psOptions->hSRS, papszArgv[i + 1]) !=
                OGRERR_NONE)
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Failed to process SRS definition: %s",
                         papszArgv[i + 1]);
                GDALRasterizeOptionsFree(psOptions);
                return nullptr;
            }

            i++;
            psOptions->bCreateOutput = true;
        }

        else if (i < argc - 4 && EQUAL(papszArgv[i], "-te"))
        {
            psOptions->sEnvelop.MinX = CPLAtof(papszArgv[++i]);
            psOptions->sEnvelop.MinY = CPLAtof(papszArgv[++i]);
            psOptions->sEnvelop.MaxX = CPLAtof(papszArgv[++i]);
            psOptions->sEnvelop.MaxY = CPLAtof(papszArgv[++i]);
            psOptions->bCreateOutput = true;
        }
        else if (i < argc - 4 && EQUAL(papszArgv[i], "-a_ullr"))
        {
            psOptions->sEnvelop.MinX = CPLAtof(papszArgv[++i]);
            psOptions->sEnvelop.MaxY = CPLAtof(papszArgv[++i]);
            psOptions->sEnvelop.MaxX = CPLAtof(papszArgv[++i]);
            psOptions->sEnvelop.MinY = CPLAtof(papszArgv[++i]);
            psOptions->bCreateOutput = true;
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-co"))
        {
            psOptions->papszCreationOptions =
                CSLAddString(psOptions->papszCreationOptions, papszArgv[++i]);
            psOptions->bCreateOutput = true;
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-ot"))
        {
            for (int iType = 1; iType < GDT_TypeCount; iType++)
            {
                const GDALDataType eType = static_cast<GDALDataType>(iType);
                if (GDALGetDataTypeName(eType) != nullptr &&
                    EQUAL(GDALGetDataTypeName(eType), papszArgv[i + 1]))
                {
                    psOptions->eOutputType = eType;
                }
            }

            if (psOptions->eOutputType == GDT_Unknown)
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Unknown output pixel type: %s", papszArgv[i + 1]);
                GDALRasterizeOptionsFree(psOptions);
                return nullptr;
            }
            i++;
            psOptions->bCreateOutput = true;
        }
        else if (i < argc - 2 && (EQUAL(papszArgv[i], "-ts") ||
                                  EQUAL(papszArgv[i], "-outsize")))
        {
            psOptions->nXSize = atoi(papszArgv[++i]);
            psOptions->nYSize = atoi(papszArgv[++i]);
            if (psOptions->nXSize <= 0 || psOptions->nYSize <= 0)
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Wrong value for -ts parameter.");
                GDALRasterizeOptionsFree(psOptions);
                return nullptr;
            }
            psOptions->bCreateOutput = true;
        }
        else if (i < argc - 2 && EQUAL(papszArgv[i], "-tr"))
        {
            psOptions->dfXRes = CPLAtof(papszArgv[++i]);
            psOptions->dfYRes = fabs(CPLAtof(papszArgv[++i]));
            if (psOptions->dfXRes == 0 || psOptions->dfYRes == 0)
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Wrong value for -tr parameter.");
                GDALRasterizeOptionsFree(psOptions);
                return nullptr;
            }
            psOptions->bCreateOutput = true;
        }
        else if (EQUAL(papszArgv[i], "-tap"))
        {
            psOptions->bTargetAlignedPixels = true;
            psOptions->bCreateOutput = true;
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-to"))
        {
            psOptions->papszTO =
                CSLAddString(psOptions->papszTO, papszArgv[++i]);
        }
        else if (i < argc - 1 && EQUAL(papszArgv[i], "-oo"))
        {
            i++;
            if (psOptionsForBinary)
            {
                psOptionsForBinary->aosOpenOptions.AddString(papszArgv[i]);
            }
            else
            {
                CPLError(
                    CE_Failure, CPLE_NotSupported,
                    "-oo switch only supported from gdal_rasterize binary.");
            }
        }
        else if (papszArgv[i][0] == '-')
        {
            CPLError(CE_Failure, CPLE_NotSupported, "Unknown option name '%s'",
                     papszArgv[i]);
            GDALRasterizeOptionsFree(psOptions);
            return nullptr;
        }
        else if (!bGotSourceFilename)
        {
            bGotSourceFilename = true;
            if (psOptionsForBinary)
            {
                psOptionsForBinary->osSource = papszArgv[i];
            }
            else
            {
                CPLError(CE_Failure, CPLE_NotSupported,
                         "{source_filename} only supported from gdal_rasterize "
                         "binary.");
            }
        }
        else if (!bGotDestFilename)
        {
            bGotDestFilename = true;
            if (psOptionsForBinary)
            {
                psOptionsForBinary->bDestSpecified = true;
                psOptionsForBinary->osDest = papszArgv[i];
            }
            else
            {
                CPLError(CE_Failure, CPLE_NotSupported,
                         "{dest_filename} only supported from gdal_rasterize "
                         "binary.");
            }
        }
        else
        {
            CPLError(CE_Failure, CPLE_NotSupported,
                     "Too many command options '%s'", papszArgv[i]);
            GDALRasterizeOptionsFree(psOptions);
            return nullptr;
        }
    }

    int nExclusiveOptionsCount = 0;
    nExclusiveOptionsCount += !psOptions->adfBurnValues.empty() ? 1 : 0;
    nExclusiveOptionsCount += psOptions->pszBurnAttribute != nullptr ? 1 : 0;
    nExclusiveOptionsCount += psOptions->b3D ? 1 : 0;
    if (nExclusiveOptionsCount != 1)
    {
        if (nExclusiveOptionsCount == 0 && psOptionsForBinary == nullptr)
        {
            psOptions->adfBurnValues.push_back(255);
        }
        else
        {
            CPLError(CE_Failure, CPLE_NotSupported,
                     "One and only one of -3d, -burn or -a is required.");
            GDALRasterizeOptionsFree(psOptions);
            return nullptr;
        }
    }

    if (psOptions->bCreateOutput)
    {
        if (psOptions->dfXRes == 0 && psOptions->dfYRes == 0 &&
            psOptions->nXSize == 0 && psOptions->nYSize == 0)
        {
            CPLError(CE_Failure, CPLE_NotSupported,
                     "'-tr xres yres' or '-ts xsize ysize' is required.");
            GDALRasterizeOptionsFree(psOptions);
            return nullptr;
        }

        if (psOptions->bTargetAlignedPixels && psOptions->dfXRes == 0 &&
            psOptions->dfYRes == 0)
        {
            CPLError(CE_Failure, CPLE_NotSupported,
                     "-tap option cannot be used without using -tr.");
            GDALRasterizeOptionsFree(psOptions);
            return nullptr;
        }

        if (!psOptions->anBandList.empty())
        {
            CPLError(CE_Failure, CPLE_NotSupported,
                     "-b option cannot be used when creating a GDAL dataset.");
            GDALRasterizeOptionsFree(psOptions);
            return nullptr;
        }

        int nBandCount = 1;

        if (!psOptions->adfBurnValues.empty())
            nBandCount = static_cast<int>(psOptions->adfBurnValues.size());

        if (static_cast<int>(psOptions->adfInitVals.size()) > nBandCount)
            nBandCount = static_cast<int>(psOptions->adfInitVals.size());

        if (psOptions->adfInitVals.size() == 1)
        {
            for (int i = 1; i <= nBandCount - 1; i++)
                psOptions->adfInitVals.push_back(psOptions->adfInitVals[0]);
        }

        for (int i = 1; i <= nBandCount; i++)
            psOptions->anBandList.push_back(i);
    }
    else
    {
        if (psOptions->anBandList.empty())
            psOptions->anBandList.push_back(1);
    }

    if (psOptions->pszDialect != nullptr && psOptions->pszWHERE != nullptr &&
        psOptions->pszSQL == nullptr)
    {
        CPLError(CE_Warning, CPLE_AppDefined,
                 "-dialect is ignored with -where. Use -sql instead");
    }

    if (psOptionsForBinary)
    {
        psOptionsForBinary->bCreateOutput = psOptions->bCreateOutput;
        if (psOptions->pszFormat)
            psOptionsForBinary->osFormat = psOptions->pszFormat;
    }

    return psOptions;
}

/************************************************************************/
/*                       GDALRasterizeOptionsFree()                     */
/************************************************************************/

/**
 * Frees the GDALRasterizeOptions struct.
 *
 * @param psOptions the options struct for GDALRasterize().
 *
 * @since GDAL 2.1
 */

void GDALRasterizeOptionsFree(GDALRasterizeOptions *psOptions)
{
    if (psOptions == nullptr)
        return;

    CPLFree(psOptions->pszFormat);
    CSLDestroy(psOptions->papszCreationOptions);
    CSLDestroy(psOptions->papszLayers);
    CSLDestroy(psOptions->papszRasterizeOptions);
    CSLDestroy(psOptions->papszTO);
    CPLFree(psOptions->pszSQL);
    CPLFree(psOptions->pszDialect);
    CPLFree(psOptions->pszBurnAttribute);
    CPLFree(psOptions->pszWHERE);
    CPLFree(psOptions->pszNoData);
    OSRDestroySpatialReference(psOptions->hSRS);

    delete psOptions;
}

/************************************************************************/
/*                  GDALRasterizeOptionsSetProgress()                   */
/************************************************************************/

/**
 * Set a progress function.
 *
 * @param psOptions the options struct for GDALRasterize().
 * @param pfnProgress the progress callback.
 * @param pProgressData the user data for the progress callback.
 *
 * @since GDAL 2.1
 */

void GDALRasterizeOptionsSetProgress(GDALRasterizeOptions *psOptions,
                                     GDALProgressFunc pfnProgress,
                                     void *pProgressData)
{
    psOptions->pfnProgress = pfnProgress ? pfnProgress : GDALDummyProgress;
    psOptions->pProgressData = pProgressData;
}