/******************************************************************************
 *
 * Project:  GDAL Utilities
 * Purpose:  Create an OGR datasource from the values of a GDAL dataset
 *           May be useful to test gdal_grid and generate its input OGR file
 * Author:   Even Rouault, <even dot rouault at spatialys.com>
 *
 ******************************************************************************
 * Copyright (c) 2008, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "gdal.h"
#include "ogr_api.h"
#include "ogr_srs_api.h"
#include "cpl_string.h"

/************************************************************************/
/*                               Usage()                                */
/************************************************************************/

void Usage()
{
    int iDriver;
    int nDriverCount;

    printf("Usage: gdal2ogr [--help] [--help-general] [-f format_name]\n"
           "                [-b band_number] [-l dest_layer_name]\n"
           "                [-t type]\n"
           "                gdal_datasource_src_name ogr_datasource_dst_name\n"
           "\n"
           " -f format_name: output file format name, possible values are:\n");

    nDriverCount = OGRGetDriverCount();
    for (iDriver = 0; iDriver < nDriverCount; iDriver++)
    {
        OGRSFDriverH hDriver = OGRGetDriver(iDriver);

        if (OGR_Dr_TestCapability(hDriver, ODrCCreateDataSource))
            printf("     -f \"%s\"\n", OGR_Dr_GetName(hDriver));
    }

    printf(
        " -b band_number: band number of the GDAL datasource (1 by default)\n"
        " -l dest_layer_name : name of the layer created in the OGR "
        "datasource\n"
        "                      (basename of the OGR datasource by default)\n"
        " -t type: one of POINT, POINT25D (default), POLYGON\n"
        "\n"
        "Create an OGR datasource from the values of a GDAL dataset.\n\n");

    exit(1);
}

/************************************************************************/
/*                                main()                                */
/************************************************************************/

int main(int argc, char *argv[])
{
    const char *pszFormat = "ESRI Shapefile";
    char *pszLayerName = NULL;
    const char *pszSrcFilename = NULL, *pszDstFilename = NULL;
    int iBand = 1;
    GDALDatasetH hDS;
    GDALRasterBandH hBand;
    int nXSize, nYSize;
    int i, j;
    FILE *fOut = NULL;
    double *padfBuffer;
    double adfGeotransform[6];
    OGRSFDriverH hOGRDriver;
    OGRDataSourceH hOGRDS;
    OGRLayerH hOGRLayer;
    OGRwkbGeometryType eType = wkbPoint25D;
    int xStep = 1, yStep = 1;

    OGRRegisterAll();
    GDALAllRegister();

    argc = GDALGeneralCmdLineProcessor(argc, &argv, 0);
    if (argc < 1)
        exit(-argc);

    /* -------------------------------------------------------------------- */
    /*      Parse arguments.                                                */
    /* -------------------------------------------------------------------- */
    for (i = 1; i < argc; i++)
    {
        if (EQUAL(argv[i], "--help"))
            Usage();
        else if (EQUAL(argv[i], "-b") && i < argc - 1)
            iBand = atoi(argv[++i]);
        else if (EQUAL(argv[i], "-f") && i < argc - 1)
            pszFormat = argv[++i];
        else if (EQUAL(argv[i], "-l") && i < argc - 1)
            pszLayerName = CPLStrdup(argv[++i]);
        else if (EQUAL(argv[i], "-t") && i < argc - 1)
        {
            i++;
            if (EQUAL(argv[i], "POLYGON"))
                eType = wkbPolygon;
            else if (EQUAL(argv[i], "POINT"))
                eType = wkbPoint;
            else if (EQUAL(argv[i], "POINT25D"))
                eType = wkbPoint25D;
            else
            {
                fprintf(stderr, "unhandled geometry type : %s\n", argv[i]);
            }
        }
        else if (EQUAL(argv[i], "-step") && i < argc - 1)
            xStep = yStep = atoi(argv[++i]);
        else if (argv[i][0] == '-')
            Usage();
        else if (pszSrcFilename == NULL)
            pszSrcFilename = argv[i];
        else if (pszDstFilename == NULL)
            pszDstFilename = argv[i];
        else
            Usage();
    }

    if (pszSrcFilename == NULL || pszDstFilename == NULL)
        Usage();

    /* -------------------------------------------------------------------- */
    /*      Open GDAL source dataset                                        */
    /* -------------------------------------------------------------------- */
    hDS = GDALOpen(pszSrcFilename, GA_ReadOnly);
    if (hDS == NULL)
    {
        fprintf(stderr, "Can't open %s\n", pszSrcFilename);
        exit(1);
    }

    hBand = GDALGetRasterBand(hDS, iBand);
    if (hBand == NULL)
    {
        fprintf(stderr, "Can't get band %d\n", iBand);
        exit(1);
    }

    if (GDALGetGeoTransform(hDS, adfGeotransform) != CE_None)
    {
        fprintf(stderr, "Can't get geotransform\n");
        exit(1);
    }

    nXSize = GDALGetRasterXSize(hDS);
    nYSize = GDALGetRasterYSize(hDS);

    /* -------------------------------------------------------------------- */
    /*     Create OGR destination dataset                                   */
    /* -------------------------------------------------------------------- */
    /* Special case for CSV : we generate the appropriate VRT file in the same
     * time */
    if (EQUAL(pszFormat, "CSV") &&
        EQUAL(CPLGetExtension(pszDstFilename), "CSV"))
    {
        FILE *fOutCSVT;
        FILE *fOutVRT;
        char *pszDstFilenameCSVT;
        char *pszDstFilenameVRT;

        fOut = fopen(pszDstFilename, "wt");
        if (fOut == NULL)
        {
            fprintf(stderr, "Can't open %s for writing\n", pszDstFilename);
            exit(1);
        }
        fprintf(fOut, "x,y,z\n");

        pszDstFilenameCSVT = CPLMalloc(strlen(pszDstFilename) + 2);
        strcpy(pszDstFilenameCSVT, pszDstFilename);
        strcat(pszDstFilenameCSVT, "t");
        fOutCSVT = fopen(pszDstFilenameCSVT, "wt");
        if (fOutCSVT == NULL)
        {
            fprintf(stderr, "Can't open %s for writing\n", pszDstFilenameCSVT);
            exit(1);
        }
        CPLFree(pszDstFilenameCSVT);
        fprintf(fOutCSVT, "Real,Real,Real\n");
        fclose(fOutCSVT);
        fOutCSVT = NULL;

        pszDstFilenameVRT = CPLStrdup(pszDstFilename);
        strcpy(pszDstFilenameVRT + strlen(pszDstFilename) - 3, "vrt");
        fOutVRT = fopen(pszDstFilenameVRT, "wt");
        if (fOutVRT == NULL)
        {
            fprintf(stderr, "Can't open %s for writing\n", pszDstFilenameVRT);
            exit(1);
        }
        CPLFree(pszDstFilenameVRT);
        fprintf(fOutVRT, "<OGRVRTDataSource>\n");
        fprintf(fOutVRT, "  <OGRVRTLayer name=\"%s\">\n",
                CPLGetBasename(pszDstFilename));
        fprintf(fOutVRT, "    <SrcDataSource>%s</SrcDataSource> \n",
                pszDstFilename);
        fprintf(fOutVRT, "    <GeometryType>wkbPoint</GeometryType>\n");
        fprintf(fOutVRT, "    <GeometryField encoding=\"PointFromColumns\" "
                         "x=\"x\" y=\"y\" z=\"z\"/>\n");
        fprintf(fOutVRT, "  </OGRVRTLayer>\n");
        fprintf(fOutVRT, "</OGRVRTDataSource>\n");
        fclose(fOutVRT);
        fOutVRT = NULL;
    }
    else
    {
        OGRSpatialReferenceH hSRS = NULL;
        const char *pszWkt;

        hOGRDriver = OGRGetDriverByName(pszFormat);
        if (hOGRDriver == NULL)
        {
            fprintf(stderr, "Can't find OGR driver %s\n", pszFormat);
            exit(1);
        }

        hOGRDS = OGR_Dr_CreateDataSource(hOGRDriver, pszDstFilename, NULL);
        if (hOGRDS == NULL)
        {
            fprintf(stderr, "Can't create OGR datasource %s\n", pszDstFilename);
            exit(1);
        }

        pszWkt = GDALGetProjectionRef(hDS);
        if (pszWkt && pszWkt[0])
        {
            hSRS = OSRNewSpatialReference(pszWkt);
        }

        if (pszLayerName == NULL)
            pszLayerName = CPLStrdup(CPLGetBasename(pszDstFilename));

        hOGRLayer = OGR_DS_CreateLayer(hOGRDS, pszLayerName, hSRS, eType, NULL);

        if (hSRS)
            OSRDestroySpatialReference(hSRS);

        if (hOGRLayer == NULL)
        {
            fprintf(stderr, "Can't create layer %s\n", pszLayerName);
            exit(1);
        }

        if (eType != wkbPoint25D)
        {
            OGRFieldDefnH hFieldDefn = OGR_Fld_Create("z", OFTReal);
            OGR_L_CreateField(hOGRLayer, hFieldDefn, 0);
            OGR_Fld_Destroy(hFieldDefn);
        }
    }

    padfBuffer = (double *)CPLMalloc(nXSize * sizeof(double));

#define GET_X(j, i)                                                            \
    adfGeotransform[0] + (j)*adfGeotransform[1] + (i)*adfGeotransform[2]
#define GET_Y(j, i)                                                            \
    adfGeotransform[3] + (j)*adfGeotransform[4] + (i)*adfGeotransform[5]
#define GET_XY(j, i) GET_X(j, i), GET_Y(j, i)

    /* -------------------------------------------------------------------- */
    /*     "Translate" the source dataset                                   */
    /* -------------------------------------------------------------------- */
    for (i = 0; i < nYSize; i += yStep)
    {
        GDALRasterIO(hBand, GF_Read, 0, i, nXSize, 1, padfBuffer, nXSize, 1,
                     GDT_Float64, 0, 0);
        for (j = 0; j < nXSize; j += xStep)
        {
            if (fOut)
            {
                fprintf(fOut, "%f,%f,%f\n", GET_XY(j + .5, i + .5),
                        padfBuffer[j]);
            }
            else
            {
                OGRFeatureH hFeature =
                    OGR_F_Create(OGR_L_GetLayerDefn(hOGRLayer));
                OGRGeometryH hGeometry = OGR_G_CreateGeometry(eType);
                if (eType == wkbPoint25D)
                {
                    OGR_G_SetPoint(hGeometry, 0, GET_XY(j + .5, i + .5),
                                   padfBuffer[j]);
                }
                else if (eType == wkbPoint)
                {
                    OGR_G_SetPoint_2D(hGeometry, 0, GET_XY(j + .5, i + .5));
                    OGR_F_SetFieldDouble(hFeature, 0, padfBuffer[j]);
                }
                else
                {
                    OGRGeometryH hLinearRing =
                        OGR_G_CreateGeometry(wkbLinearRing);
                    OGR_G_SetPoint_2D(hLinearRing, 0, GET_XY(j + 0, i + 0));
                    OGR_G_SetPoint_2D(hLinearRing, 1, GET_XY(j + 1, i + 0));
                    OGR_G_SetPoint_2D(hLinearRing, 2, GET_XY(j + 1, i + 1));
                    OGR_G_SetPoint_2D(hLinearRing, 3, GET_XY(j + 0, i + 1));
                    OGR_G_SetPoint_2D(hLinearRing, 4, GET_XY(j + 0, i + 0));
                    OGR_G_AddGeometryDirectly(hGeometry, hLinearRing);
                    OGR_F_SetFieldDouble(hFeature, 0, padfBuffer[j]);
                }
                OGR_F_SetGeometryDirectly(hFeature, hGeometry);
                OGR_L_CreateFeature(hOGRLayer, hFeature);
                OGR_F_Destroy(hFeature);
            }
        }
    }

    /* -------------------------------------------------------------------- */
    /*     Cleanup                                                          */
    /* -------------------------------------------------------------------- */
    if (fOut)
        fclose(fOut);
    else
        OGR_DS_Destroy(hOGRDS);
    GDALClose(hDS);

    CPLFree(padfBuffer);
    CPLFree(pszLayerName);

    GDALDumpOpenDatasets(stderr);
    GDALDestroyDriverManager();
    OGRCleanupAll();
    CSLDestroy(argv);

    return 0;
}