/******************************************************************************
 *
 * Project:  GDAL
 * Purpose:  Command line raster query tool.
 * Author:   Frank Warmerdam <warmerdam@pobox.com>
 *
 ******************************************************************************
 * Copyright (c) 2010, Frank Warmerdam <warmerdam@pobox.com>
 * Copyright (c) 2010-2013, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_string.h"
#include "cpl_minixml.h"
#include "gdal_version.h"
#include "gdal.h"
#include "gdal_priv.h"
#include "commonutils.h"
#include "ogr_spatialref.h"
#include "gdalargumentparser.h"

#include <cmath>
#include <limits>
#include <vector>

#include <cctype>

/************************************************************************/
/*                             GetSRSAsWKT                              */
/************************************************************************/

static std::string GetSRSAsWKT(const char *pszUserInput)

{
    OGRSpatialReference oSRS;
    oSRS.SetFromUserInput(pszUserInput);
    return oSRS.exportToWkt();
}

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

MAIN_START(argc, argv)

{
    double dfGeoX = std::numeric_limits<double>::quiet_NaN();
    double dfGeoY = std::numeric_limits<double>::quiet_NaN();
    std::string osSrcFilename;
    std::string osSourceSRS;
    std::vector<int> anBandList;
    bool bAsXML = false, bLIFOnly = false;
    bool bQuiet = false, bValOnly = false;
    int nOverview = 0;
    CPLStringList aosOpenOptions;
    std::string osFieldSep;
    bool bIgnoreExtraInput = false;
    bool bEcho = false;

    GDALAllRegister();
    argc = GDALGeneralCmdLineProcessor(argc, &argv, 0);
    if (argc < 1)
        exit(-argc);
    CPLStringList aosArgv;
    aosArgv.Assign(argv, /* bAssign = */ true);

    GDALArgumentParser argParser(aosArgv[0], /* bForBinary=*/true);

    argParser.add_description(_("Raster query tool."));

    const char *pszEpilog =
        _("For more details, consult "
          "https://gdal.org/programs/gdallocationinfo.html");
    argParser.add_epilog(pszEpilog);

    argParser.add_argument("-xml").flag().store_into(bAsXML).help(
        _("Format the output report as XML."));

    argParser.add_argument("-lifonly")
        .flag()
        .store_into(bLIFOnly)
        .help(_("Only outputs filenames from the LocationInfo request against "
                "the database."));

    argParser.add_argument("-valonly")
        .flag()
        .store_into(bValOnly)
        .help(_("Only outputs pixel values of the selected pixel on each of "
                "the selected bands."));

    argParser.add_argument("-E").flag().store_into(bEcho).help(
        _("Enable Echo mode, where input coordinates are prepended to the "
          "output lines in -valonly mode."));

    argParser.add_argument("-field_sep")
        .metavar("<sep>")
        .store_into(osFieldSep)
        .help(_("Defines the field separator, used in -valonly mode, to "
                "separate different values."));

    argParser.add_argument("-ignore_extra_input")
        .flag()
        .store_into(bIgnoreExtraInput)
        .help(_("Set this flag to avoid extra non-numeric content at end of "
                "input lines."));

    argParser.add_argument("-b")
        .append()
        .metavar("<band>")
        .store_into(anBandList)
        .help(_("Select band(s)."));

    argParser.add_argument("-overview")
        .metavar("<overview_level>")
        .store_into(nOverview)
        .help(_("Query the (overview_level)th overview (overview_level=1 is "
                "the 1st overview)."));

    std::string osResampling;
    argParser.add_argument("-r")
        .store_into(osResampling)
        .metavar("nearest|bilinear|cubic|cubicspline")
        .help(_("Select an interpolation algorithm."));

    {
        auto &group = argParser.add_mutually_exclusive_group();

        group.add_argument("-l_srs")
            .metavar("<srs_def>")
            .store_into(osSourceSRS)
            .help(_("Coordinate system of the input x, y location."));

        group.add_argument("-geoloc")
            .flag()
            .action([&osSourceSRS](const std::string &)
                    { osSourceSRS = "-geoloc"; })
            .help(_("Indicates input x,y points are in the georeferencing "
                    "system of the image."));

        group.add_argument("-wgs84")
            .flag()
            .action([&osSourceSRS](const std::string &)
                    { osSourceSRS = GetSRSAsWKT("WGS84"); })
            .help(_("Indicates input x,y points are WGS84 long, lat."));
    }

    argParser.add_open_options_argument(&aosOpenOptions);

    argParser.add_argument("srcfile")
        .metavar("<srcfile>")
        .nargs(1)
        .store_into(osSrcFilename)
        .help(_("The source GDAL raster datasource name."));

    argParser.add_argument("x")
        .metavar("<x>")
        .nargs(argparse::nargs_pattern::optional)
        .store_into(dfGeoX)
        .help(_("X location of target pixel."));

    argParser.add_argument("y")
        .metavar("<y>")
        .nargs(argparse::nargs_pattern::optional)
        .store_into(dfGeoY)
        .help(_("Y location of target pixel."));

    const auto displayUsage = [&argParser]()
    {
        std::stringstream usageStringStream;
        usageStringStream << argParser.usage();
        std::cerr << CPLString(usageStringStream.str())
                         .replaceAll("<x> <y>", "[<x> <y>]")
                  << std::endl
                  << std::endl;
        std::cout << _("Note: ") << "gdallocationinfo"
                  << _(" --long-usage for full help.") << std::endl;
    };

    try
    {
        argParser.parse_args(aosArgv);
    }
    catch (const std::exception &err)
    {
        std::cerr << _("Error: ") << err.what() << std::endl;
        displayUsage();
        std::exit(1);
    }

    if (bLIFOnly || bValOnly)
        bQuiet = true;

    // User specifies with 1-based index, but internally we use 0-based index
    --nOverview;

    // Deal with special characters
    osFieldSep = CPLString(osFieldSep)
                     .replaceAll("\\t", '\t')
                     .replaceAll("\\r", '\r')
                     .replaceAll("\\n", '\n');

    if (!std::isnan(dfGeoX) && std::isnan(dfGeoY))
    {
        fprintf(stderr, "<y> should be specified when <x> is specified\n\n");
        displayUsage();
        exit(1);
    }

    const bool bIsXYSpecifiedAsArgument = !std::isnan(dfGeoX);

    if (bEcho && !bValOnly)
    {
        fprintf(stderr, "-E can only be used with -valonly\n");
        exit(1);
    }
    if (bEcho && osFieldSep.empty())
    {
        fprintf(stderr, "-E can only be used if -field_sep is specified (to a "
                        "non-newline value)\n");
        exit(1);
    }

    if (osFieldSep.empty())
    {
        osFieldSep = "\n";
    }
    else if (!bValOnly)
    {
        fprintf(stderr, "-field_sep can only be used with -valonly\n");
        exit(1);
    }

    const GDALRIOResampleAlg eInterpolation =
        osResampling.empty() ? GRIORA_NearestNeighbour
                             : GDALRasterIOGetResampleAlg(osResampling.c_str());
    if (eInterpolation != GRIORA_NearestNeighbour &&
        eInterpolation != GRIORA_Bilinear && eInterpolation != GRIORA_Cubic &&
        eInterpolation != GRIORA_CubicSpline)
    {
        fprintf(stderr, "-r can only be used with values nearest, bilinear, "
                        "cubic and cubicspline\n");
        exit(1);
    }

    /* -------------------------------------------------------------------- */
    /*      Open source file.                                               */
    /* -------------------------------------------------------------------- */
    GDALDatasetH hSrcDS = GDALOpenEx(osSrcFilename.c_str(),
                                     GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
                                     nullptr, aosOpenOptions.List(), nullptr);
    if (hSrcDS == nullptr)
        exit(1);

    /* -------------------------------------------------------------------- */
    /*      Setup coordinate transformation, if required                    */
    /* -------------------------------------------------------------------- */
    OGRSpatialReferenceH hSrcSRS = nullptr;
    OGRCoordinateTransformationH hCT = nullptr;
    if (!osSourceSRS.empty() && !EQUAL(osSourceSRS.c_str(), "-geoloc"))
    {
        hSrcSRS = OSRNewSpatialReference(nullptr);
        OGRErr err = OSRSetFromUserInput(hSrcSRS, osSourceSRS.c_str());
        if (err != OGRERR_NONE)
            exit(1);
        OSRSetAxisMappingStrategy(hSrcSRS, OAMS_TRADITIONAL_GIS_ORDER);
        auto hTrgSRS = GDALGetSpatialRef(hSrcDS);
        if (!hTrgSRS)
            exit(1);

        hCT = OCTNewCoordinateTransformation(hSrcSRS, hTrgSRS);
        if (hCT == nullptr)
            exit(1);
    }

    /* -------------------------------------------------------------------- */
    /*      If no bands were requested, we will query them all.             */
    /* -------------------------------------------------------------------- */
    if (anBandList.empty())
    {
        for (int i = 0; i < GDALGetRasterCount(hSrcDS); i++)
            anBandList.push_back(i + 1);
    }

    /* -------------------------------------------------------------------- */
    /*      Turn the location into a pixel and line location.               */
    /* -------------------------------------------------------------------- */
    bool inputAvailable = true;
    CPLString osXML;
    char szLine[1024];
    int nLine = 0;
    std::string osExtraContent;

    if (std::isnan(dfGeoX))
    {
        // Is it an interactive terminal ?
        if (CPLIsInteractive(stdin))
        {
            if (!osSourceSRS.empty())
            {
                fprintf(stderr, "Enter X Y values separated by space, and "
                                "press Return.\n");
            }
            else
            {
                fprintf(stderr, "Enter pixel line values separated by space, "
                                "and press Return.\n");
            }
        }

        if (fgets(szLine, sizeof(szLine) - 1, stdin))
        {
            const CPLStringList aosTokens(CSLTokenizeString(szLine));
            const int nCount = aosTokens.size();

            ++nLine;
            if (nCount < 2)
            {
                fprintf(stderr, "Not enough values at line %d\n", nLine);
                inputAvailable = false;
            }
            else
            {
                dfGeoX = CPLAtof(aosTokens[0]);
                dfGeoY = CPLAtof(aosTokens[1]);
                if (!bIgnoreExtraInput)
                {
                    for (int i = 2; i < nCount; ++i)
                    {
                        if (!osExtraContent.empty())
                            osExtraContent += ' ';
                        osExtraContent += aosTokens[i];
                    }
                    while (!osExtraContent.empty() &&
                           isspace(static_cast<int>(osExtraContent.back())))
                    {
                        osExtraContent.pop_back();
                    }
                }
            }
        }
        else
        {
            inputAvailable = false;
        }
    }

    int nRetCode = 0;
    while (inputAvailable)
    {
        int iPixel, iLine;
        double dfPixel{0}, dfLine{0};
        const double dfXIn = dfGeoX;
        const double dfYIn = dfGeoY;

        if (hCT)
        {
            if (!OCTTransform(hCT, 1, &dfGeoX, &dfGeoY, nullptr))
                exit(1);
        }

        if (!osSourceSRS.empty())
        {
            double adfGeoTransform[6] = {};
            if (GDALGetGeoTransform(hSrcDS, adfGeoTransform) != CE_None)
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Cannot get geotransform");
                exit(1);
            }

            double adfInvGeoTransform[6] = {};
            if (!GDALInvGeoTransform(adfGeoTransform, adfInvGeoTransform))
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Cannot invert geotransform");
                exit(1);
            }

            dfPixel = adfInvGeoTransform[0] + adfInvGeoTransform[1] * dfGeoX +
                      adfInvGeoTransform[2] * dfGeoY;
            dfLine = adfInvGeoTransform[3] + adfInvGeoTransform[4] * dfGeoX +
                     adfInvGeoTransform[5] * dfGeoY;
        }
        else
        {
            dfPixel = dfGeoX;
            dfLine = dfGeoY;
        }
        iPixel = static_cast<int>(floor(dfPixel));
        iLine = static_cast<int>(floor(dfLine));

        /* --------------------------------------------------------------------
         */
        /*      Prepare report. */
        /* --------------------------------------------------------------------
         */
        CPLString osXmlLine;

        if (bAsXML)
        {
            osXmlLine.Printf("<Report pixel=\"%d\" line=\"%d\">", iPixel,
                             iLine);
            osXML += osXmlLine;
            if (!osExtraContent.empty())
            {
                char *pszEscaped =
                    CPLEscapeString(osExtraContent.c_str(), -1, CPLES_XML);
                osXML += CPLString().Printf("  <ExtraInput>%s</ExtraInput>",
                                            pszEscaped);
                CPLFree(pszEscaped);
            }
        }
        else if (!bQuiet)
        {
            printf("Report:\n");
            CPLString osPixel, osLine;
            if (eInterpolation == GRIORA_NearestNeighbour)
            {
                osPixel.Printf("%d", iPixel);
                osLine.Printf("%d", iLine);
            }
            else
            {
                osPixel.Printf("%.15g", dfPixel);
                osLine.Printf("%.15g", dfLine);
            }
            printf("  Location: (%sP,%sL)\n", osPixel.c_str(), osLine.c_str());
            if (!osExtraContent.empty())
            {
                printf("  Extra input: %s\n", osExtraContent.c_str());
            }
        }
        else if (bEcho)
        {
            printf("%.15g%s%.15g%s", dfXIn, osFieldSep.c_str(), dfYIn,
                   osFieldSep.c_str());
        }

        bool bPixelReport = true;

        if (iPixel < 0 || iLine < 0 ||
            dfPixel > static_cast<double>(GDALGetRasterXSize(hSrcDS) + 1e-5) ||
            dfLine > static_cast<double>(GDALGetRasterYSize(hSrcDS) + 1e-5))
        {
            if (bAsXML)
                osXML += "<Alert>Location is off this file! No further details "
                         "to report.</Alert>";
            else if (bValOnly)
            {
                for (int i = 1; i < static_cast<int>(anBandList.size()); i++)
                {
                    printf("%s", osFieldSep.c_str());
                }
            }
            else if (!bQuiet)
                printf("\nLocation is off this file! No further details to "
                       "report.\n");
            bPixelReport = false;
            nRetCode = 1;
        }

        /* --------------------------------------------------------------------
         */
        /*      Process each band. */
        /* --------------------------------------------------------------------
         */
        for (int i = 0; bPixelReport && i < static_cast<int>(anBandList.size());
             i++)
        {
            GDALRasterBandH hBand = GDALGetRasterBand(hSrcDS, anBandList[i]);

            int iPixelToQuery = iPixel;
            int iLineToQuery = iLine;

            double dfPixelToQuery = dfPixel;
            double dfLineToQuery = dfLine;

            if (nOverview >= 0 && hBand != nullptr)
            {
                GDALRasterBandH hOvrBand = GDALGetOverview(hBand, nOverview);
                if (hOvrBand != nullptr)
                {
                    int nOvrXSize = GDALGetRasterBandXSize(hOvrBand);
                    int nOvrYSize = GDALGetRasterBandYSize(hOvrBand);
                    iPixelToQuery = static_cast<int>(
                        0.5 +
                        1.0 * iPixel / GDALGetRasterXSize(hSrcDS) * nOvrXSize);
                    iLineToQuery = static_cast<int>(
                        0.5 +
                        1.0 * iLine / GDALGetRasterYSize(hSrcDS) * nOvrYSize);
                    if (iPixelToQuery >= nOvrXSize)
                        iPixelToQuery = nOvrXSize - 1;
                    if (iLineToQuery >= nOvrYSize)
                        iLineToQuery = nOvrYSize - 1;
                    dfPixelToQuery =
                        dfPixel / GDALGetRasterXSize(hSrcDS) * nOvrXSize;
                    dfLineToQuery =
                        dfLine / GDALGetRasterYSize(hSrcDS) * nOvrYSize;
                }
                else
                {
                    CPLError(CE_Failure, CPLE_AppDefined,
                             "Cannot get overview %d of band %d", nOverview + 1,
                             anBandList[i]);
                }
                hBand = hOvrBand;
            }

            if (hBand == nullptr)
                continue;

            if (bAsXML)
            {
                osXmlLine.Printf("<BandReport band=\"%d\">", anBandList[i]);
                osXML += osXmlLine;
            }
            else if (!bQuiet)
            {
                printf("  Band %d:\n", anBandList[i]);
            }

            /* --------------------------------------------------------------------
             */
            /*      Request location info for this location.  It is possible */
            /*      only the VRT driver actually supports this. */
            /* --------------------------------------------------------------------
             */
            CPLString osItem;

            osItem.Printf("Pixel_%d_%d", iPixelToQuery, iLineToQuery);

            const char *pszLI =
                GDALGetMetadataItem(hBand, osItem, "LocationInfo");

            if (pszLI != nullptr)
            {
                if (bAsXML)
                    osXML += pszLI;
                else if (!bQuiet)
                    printf("    %s\n", pszLI);
                else if (bLIFOnly)
                {
                    /* Extract all files, if any. */

                    CPLXMLNode *psRoot = CPLParseXMLString(pszLI);

                    if (psRoot != nullptr && psRoot->psChild != nullptr &&
                        psRoot->eType == CXT_Element &&
                        EQUAL(psRoot->pszValue, "LocationInfo"))
                    {
                        for (CPLXMLNode *psNode = psRoot->psChild;
                             psNode != nullptr; psNode = psNode->psNext)
                        {
                            if (psNode->eType == CXT_Element &&
                                EQUAL(psNode->pszValue, "File") &&
                                psNode->psChild != nullptr)
                            {
                                char *pszUnescaped =
                                    CPLUnescapeString(psNode->psChild->pszValue,
                                                      nullptr, CPLES_XML);
                                printf("%s\n", pszUnescaped);
                                CPLFree(pszUnescaped);
                            }
                        }
                    }
                    CPLDestroyXMLNode(psRoot);
                }
            }

            /* --------------------------------------------------------------------
             */
            /*      Report the pixel value of this band. */
            /* --------------------------------------------------------------------
             */
            double adfPixel[2] = {0, 0};
            const bool bIsComplex = CPL_TO_BOOL(
                GDALDataTypeIsComplex(GDALGetRasterDataType(hBand)));

            CPLErrorReset();
            CPLErr err = GDALRasterInterpolateAtPoint(
                hBand, dfPixelToQuery, dfLineToQuery, eInterpolation,
                &adfPixel[0], &adfPixel[1]);

            // GDALRasterInterpolateAtPoint() returns false on nodata
            bool bIsNoData = false;
            if (err != CE_None && CPLGetLastErrorType() == CE_None)
            {
                int bHasNoData = FALSE;
                const double dfNoData =
                    GDALGetRasterNoDataValue(hBand, &bHasNoData);
                if (bHasNoData)
                {
                    bIsNoData = true;
                    adfPixel[0] = adfPixel[1] = dfNoData;
                    err = CE_None;
                }
            }
            if (err == CE_None)
            {
                CPLString osValue;

                if (bIsComplex)
                    osValue.Printf("%.15g+%.15gi", adfPixel[0], adfPixel[1]);
                else
                    osValue.Printf("%.15g", adfPixel[0]);

                if (bAsXML)
                {
                    osXML += "<Value>";
                    osXML += osValue;
                    osXML += "</Value>";
                }
                else if (!bQuiet)
                    printf("    Value: %s\n", osValue.c_str());
                else if (bValOnly)
                {
                    if (i > 0)
                        printf("%s", osFieldSep.c_str());
                    printf("%s", osValue.c_str());
                }

                // Report unscaled if we have scale/offset values.
                int bSuccess;

                double dfOffset = GDALGetRasterOffset(hBand, &bSuccess);
                // TODO: Should we turn on checking of bSuccess?
                // Alternatively, delete these checks and put a comment as to
                // why checking bSuccess does not matter.
#if 0
                if (bSuccess == FALSE)
                {
                    CPLError( CE_Debug, CPLE_AppDefined,
                              "Unable to get raster offset." );
                }
#endif
                double dfScale = GDALGetRasterScale(hBand, &bSuccess);
#if 0
                if (bSuccess == FALSE)
                {
                    CPLError( CE_Debug, CPLE_AppDefined,
                              "Unable to get raster scale." );
                }
#endif
                if (!bIsNoData && (dfOffset != 0.0 || dfScale != 1.0))
                {
                    adfPixel[0] = adfPixel[0] * dfScale + dfOffset;

                    if (bIsComplex)
                    {
                        adfPixel[1] = adfPixel[1] * dfScale + dfOffset;
                        osValue.Printf("%.15g+%.15gi", adfPixel[0],
                                       adfPixel[1]);
                    }
                    else
                        osValue.Printf("%.15g", adfPixel[0]);

                    if (bAsXML)
                    {
                        osXML += "<DescaledValue>";
                        osXML += osValue;
                        osXML += "</DescaledValue>";
                    }
                    else if (!bQuiet)
                        printf("    Descaled Value: %s\n", osValue.c_str());
                }
            }
            else
            {
                nRetCode = 1;
                if (bAsXML)
                {
                    osXML += "<IOError/>";
                }
                else if (bValOnly)
                {
                    if (i > 0)
                        printf("%s", osFieldSep.c_str());
                }
            }
            if (bAsXML)
                osXML += "</BandReport>";
        }

        osXML += "</Report>";

        if (bValOnly)
        {
            if (!osExtraContent.empty() && osFieldSep != "\n")
                printf("%s%s", osFieldSep.c_str(), osExtraContent.c_str());
            printf("\n");
        }

        if (bIsXYSpecifiedAsArgument)
            break;

        osExtraContent.clear();
        if (fgets(szLine, sizeof(szLine) - 1, stdin))
        {
            const CPLStringList aosTokens(CSLTokenizeString(szLine));
            const int nCount = aosTokens.size();

            ++nLine;
            if (nCount < 2)
            {
                fprintf(stderr, "Not enough values at line %d\n", nLine);
                continue;
            }
            else
            {
                dfGeoX = CPLAtof(aosTokens[0]);
                dfGeoY = CPLAtof(aosTokens[1]);
                if (!bIgnoreExtraInput)
                {
                    for (int i = 2; i < nCount; ++i)
                    {
                        if (!osExtraContent.empty())
                            osExtraContent += ' ';
                        osExtraContent += aosTokens[i];
                    }
                    while (!osExtraContent.empty() &&
                           isspace(static_cast<int>(osExtraContent.back())))
                    {
                        osExtraContent.pop_back();
                    }
                }
            }
        }
        else
        {
            break;
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Finalize xml report and print.                                  */
    /* -------------------------------------------------------------------- */
    if (bAsXML)
    {
        CPLXMLNode *psRoot = CPLParseXMLString(osXML);
        char *pszFormattedXML = CPLSerializeXMLTree(psRoot);
        CPLDestroyXMLNode(psRoot);

        printf("%s", pszFormattedXML);
        CPLFree(pszFormattedXML);
    }

    /* -------------------------------------------------------------------- */
    /*      Cleanup                                                         */
    /* -------------------------------------------------------------------- */
    if (hCT)
    {
        OSRDestroySpatialReference(hSrcSRS);
        OCTDestroyCoordinateTransformation(hCT);
    }

    GDALClose(hSrcDS);

    GDALDumpOpenDatasets(stderr);
    GDALDestroyDriverManager();

    return nRetCode;
}

MAIN_END