/******************************************************************************
 *
 * Project:  Vertical Datum Transformation
 * Purpose:  Implementation of NOAA .gtx vertical datum shift file format.
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 2010, Frank Warmerdam
 * Copyright (c) 2010-2013, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_string.h"
#include "gdal_frmts.h"
#include "gdal_priv.h"
#include "ogr_srs_api.h"
#include "rawdataset.h"

#include <algorithm>
#include <limits>

/**

NOAA .GTX Vertical Datum Grid Shift Format

All values are bigendian

Header
------

float64  latitude_of_origin
float64  longitude_of_origin (0-360)
float64  cell size (x?y?)
float64  cell size (x?y?)
int32    length in pixels
int32    width in pixels

Data
----

float32  * width in pixels * length in pixels

Values are an offset in meters between two vertical datums.

**/

/************************************************************************/
/* ==================================================================== */
/*                              GTXDataset                              */
/* ==================================================================== */
/************************************************************************/

class GTXDataset final : public RawDataset
{
    VSILFILE *fpImage = nullptr;  // image data file.

    OGRSpatialReference m_oSRS{};
    GDALGeoTransform m_gt{};

    CPL_DISALLOW_COPY_ASSIGN(GTXDataset)

    CPLErr Close() override;

  public:
    GTXDataset()
    {
        m_oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
        m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    }

    ~GTXDataset() override;

    CPLErr GetGeoTransform(GDALGeoTransform &gt) const override;
    CPLErr SetGeoTransform(const GDALGeoTransform &gt) override;

    const OGRSpatialReference *GetSpatialRef() const override
    {
        return &m_oSRS;
    }

    static GDALDataset *Open(GDALOpenInfo *);
    static int Identify(GDALOpenInfo *);
    static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
                               int nBands, GDALDataType eType,
                               char **papszOptions);
};

/************************************************************************/
/* ==================================================================== */
/*                           GTXRasterBand                              */
/* ==================================================================== */
/************************************************************************/

class GTXRasterBand final : public RawRasterBand
{
    CPL_DISALLOW_COPY_ASSIGN(GTXRasterBand)

  public:
    GTXRasterBand(GDALDataset *poDS, int nBand, VSILFILE *fpRaw,
                  vsi_l_offset nImgOffset, int nPixelOffset, int nLineOffset,
                  GDALDataType eDataType, int bNativeOrder);

    ~GTXRasterBand() override;

    double GetNoDataValue(int *pbSuccess = nullptr) override;
};

/************************************************************************/
/*                            GTXRasterBand()                           */
/************************************************************************/

GTXRasterBand::GTXRasterBand(GDALDataset *poDSIn, int nBandIn,
                             VSILFILE *fpRawIn, vsi_l_offset nImgOffsetIn,
                             int nPixelOffsetIn, int nLineOffsetIn,
                             GDALDataType eDataTypeIn, int bNativeOrderIn)
    : RawRasterBand(poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
                    nLineOffsetIn, eDataTypeIn, bNativeOrderIn,
                    RawRasterBand::OwnFP::NO)
{
}

/************************************************************************/
/*                           ~GTXRasterBand()                           */
/************************************************************************/

GTXRasterBand::~GTXRasterBand()
{
}

/************************************************************************/
/*                           GetNoDataValue()                           */
/************************************************************************/

double GTXRasterBand::GetNoDataValue(int *pbSuccess)
{
    if (pbSuccess)
        *pbSuccess = TRUE;
    int bSuccess = FALSE;
    double dfNoData = GDALPamRasterBand::GetNoDataValue(&bSuccess);
    if (bSuccess)
    {
        return dfNoData;
    }
    return -88.8888;
}

/************************************************************************/
/* ==================================================================== */
/*                              GTXDataset                              */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                            ~GTXDataset()                             */
/************************************************************************/

GTXDataset::~GTXDataset()

{
    GTXDataset::Close();
}

/************************************************************************/
/*                              Close()                                 */
/************************************************************************/

CPLErr GTXDataset::Close()
{
    CPLErr eErr = CE_None;
    if (nOpenFlags != OPEN_FLAGS_CLOSED)
    {
        if (GTXDataset::FlushCache(true) != CE_None)
            eErr = CE_Failure;

        if (fpImage)
        {
            if (VSIFCloseL(fpImage) != 0)
            {
                CPLError(CE_Failure, CPLE_FileIO, "I/O error");
                eErr = CE_Failure;
            }
        }

        if (GDALPamDataset::Close() != CE_None)
            eErr = CE_Failure;
    }
    return eErr;
}

/************************************************************************/
/*                              Identify()                              */
/************************************************************************/

int GTXDataset::Identify(GDALOpenInfo *poOpenInfo)

{
    if (poOpenInfo->nHeaderBytes < 40)
        return FALSE;

    if (!poOpenInfo->IsExtensionEqualToCI("gtx"))
        return FALSE;

    return TRUE;
}

/************************************************************************/
/*                                Open()                                */
/************************************************************************/

GDALDataset *GTXDataset::Open(GDALOpenInfo *poOpenInfo)

{
    if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
        return nullptr;

    /* -------------------------------------------------------------------- */
    /*      Create a corresponding GDALDataset.                             */
    /* -------------------------------------------------------------------- */
    auto poDS = std::make_unique<GTXDataset>();

    poDS->eAccess = poOpenInfo->eAccess;
    std::swap(poDS->fpImage, poOpenInfo->fpL);

    /* -------------------------------------------------------------------- */
    /*      Read the header.                                                */
    /* -------------------------------------------------------------------- */
    double gt[6] = {0};
    CPL_IGNORE_RET_VAL(VSIFReadL(&gt[3], 8, 1, poDS->fpImage));
    CPL_IGNORE_RET_VAL(VSIFReadL(&gt[0], 8, 1, poDS->fpImage));
    CPL_IGNORE_RET_VAL(VSIFReadL(&gt[5], 8, 1, poDS->fpImage));
    CPL_IGNORE_RET_VAL(VSIFReadL(&gt[1], 8, 1, poDS->fpImage));

    CPL_IGNORE_RET_VAL(VSIFReadL(&(poDS->nRasterYSize), 4, 1, poDS->fpImage));
    CPL_IGNORE_RET_VAL(VSIFReadL(&(poDS->nRasterXSize), 4, 1, poDS->fpImage));

    CPL_MSBPTR32(&(poDS->nRasterYSize));
    CPL_MSBPTR32(&(poDS->nRasterXSize));

    CPL_MSBPTR64(&gt[0]);
    CPL_MSBPTR64(&gt[1]);
    CPL_MSBPTR64(&gt[3]);
    CPL_MSBPTR64(&gt[5]);

    poDS->m_gt = GDALGeoTransform(gt);
    poDS->m_gt[3] +=
        poDS->m_gt[5] * (static_cast<double>(poDS->nRasterYSize) - 1);

    poDS->m_gt[0] -= poDS->m_gt[1] * 0.5;
    poDS->m_gt[3] += poDS->m_gt[5] * 0.5;

    poDS->m_gt[5] *= -1;

    if (CPLFetchBool(poOpenInfo->papszOpenOptions,
                     "SHIFT_ORIGIN_IN_MINUS_180_PLUS_180", false))
    {
        if (poDS->m_gt[0] < -180.0 - poDS->m_gt[1])
            poDS->m_gt[0] += 360.0;
        else if (poDS->m_gt[0] > 180.0)
            poDS->m_gt[0] -= 360.0;
    }

    if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
        static_cast<vsi_l_offset>(poDS->nRasterXSize) * poDS->nRasterYSize >
            std::numeric_limits<vsi_l_offset>::max() / sizeof(double))
    {
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*      Guess the data type. Since October 1, 2009, it should be        */
    /*      Float32. Before it was double.                                  */
    /* -------------------------------------------------------------------- */
    CPL_IGNORE_RET_VAL(VSIFSeekL(poDS->fpImage, 0, SEEK_END));
    const vsi_l_offset nSize = VSIFTellL(poDS->fpImage);

    GDALDataType eDT = GDT_Float32;
    if (nSize - 40 == sizeof(double) *
                          static_cast<vsi_l_offset>(poDS->nRasterXSize) *
                          poDS->nRasterYSize)
        eDT = GDT_Float64;
    const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
    if (nDTSize <= 0 || poDS->nRasterXSize > INT_MAX / nDTSize)
    {
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*      Create band information object.                                 */
    /* -------------------------------------------------------------------- */
    auto poBand = std::make_unique<GTXRasterBand>(
        poDS.get(), 1, poDS->fpImage,
        static_cast<vsi_l_offset>(poDS->nRasterYSize - 1) * poDS->nRasterXSize *
                nDTSize +
            40,
        nDTSize, poDS->nRasterXSize * -nDTSize, eDT, !CPL_IS_LSB);
    if (!poBand->IsValid())
        return nullptr;
    poDS->SetBand(1, std::move(poBand));

    /* -------------------------------------------------------------------- */
    /*      Initialize any PAM information.                                 */
    /* -------------------------------------------------------------------- */
    poDS->SetDescription(poOpenInfo->pszFilename);
    poDS->TryLoadXML();

    /* -------------------------------------------------------------------- */
    /*      Check for overviews.                                            */
    /* -------------------------------------------------------------------- */
    poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);

    return poDS.release();
}

/************************************************************************/
/*                          GetGeoTransform()                           */
/************************************************************************/

CPLErr GTXDataset::GetGeoTransform(GDALGeoTransform &gt) const

{
    gt = m_gt;
    return CE_None;
}

/************************************************************************/
/*                          SetGeoTransform()                           */
/************************************************************************/

CPLErr GTXDataset::SetGeoTransform(const GDALGeoTransform &gt)
{
    if (gt[2] != 0.0 || gt[4] != 0.0)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Attempt to write skewed or rotated geotransform to gtx.");
        return CE_Failure;
    }

    m_gt = gt;

    const double dfXOrigin = m_gt[0] + 0.5 * m_gt[1];
    const double dfYOrigin = m_gt[3] + (nRasterYSize - 0.5) * m_gt[5];
    const double dfWidth = m_gt[1];
    const double dfHeight = -m_gt[5];

    unsigned char header[32] = {'\0'};
    memcpy(header + 0, &dfYOrigin, 8);
    CPL_MSBPTR64(header + 0);

    memcpy(header + 8, &dfXOrigin, 8);
    CPL_MSBPTR64(header + 8);

    memcpy(header + 16, &dfHeight, 8);
    CPL_MSBPTR64(header + 16);

    memcpy(header + 24, &dfWidth, 8);
    CPL_MSBPTR64(header + 24);

    if (VSIFSeekL(fpImage, SEEK_SET, 0) != 0 ||
        VSIFWriteL(header, 32, 1, fpImage) != 1)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Attempt to write geotransform header to GTX failed.");
        return CE_Failure;
    }

    return CE_None;
}

/************************************************************************/
/*                               Create()                               */
/************************************************************************/

GDALDataset *GTXDataset::Create(const char *pszFilename, int nXSize, int nYSize,
                                int /* nBands */, GDALDataType eType,
                                char ** /* papszOptions */)
{
    if (eType != GDT_Float32)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Attempt to create gtx file with unsupported data type '%s'.",
                 GDALGetDataTypeName(eType));
        return nullptr;
    }

    if (!EQUAL(CPLGetExtensionSafe(pszFilename).c_str(), "gtx"))
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Attempt to create gtx file with extension other than gtx.");
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*      Try to create the file.                                         */
    /* -------------------------------------------------------------------- */
    VSILFILE *fp = VSIFOpenL(pszFilename, "wb");
    if (fp == nullptr)
    {
        CPLError(CE_Failure, CPLE_OpenFailed,
                 "Attempt to create file `%s' failed.\n", pszFilename);
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*      Write out the header with stub georeferencing.                  */
    /* -------------------------------------------------------------------- */

    unsigned char header[40] = {'\0'};
    double dfYOrigin = 0.0;
    memcpy(header + 0, &dfYOrigin, 8);
    CPL_MSBPTR64(header + 0);

    double dfXOrigin = 0.0;
    memcpy(header + 8, &dfXOrigin, 8);
    CPL_MSBPTR64(header + 8);

    double dfYSize = 0.01;
    memcpy(header + 16, &dfYSize, 8);
    CPL_MSBPTR64(header + 16);

    double dfXSize = 0.01;
    memcpy(header + 24, &dfXSize, 8);
    CPL_MSBPTR64(header + 24);

    GInt32 nYSize32 = nYSize;
    memcpy(header + 32, &nYSize32, 4);
    CPL_MSBPTR32(header + 32);

    GInt32 nXSize32 = nXSize;
    memcpy(header + 36, &nXSize32, 4);
    CPL_MSBPTR32(header + 36);

    CPL_IGNORE_RET_VAL(VSIFWriteL(header, 40, 1, fp));
    CPL_IGNORE_RET_VAL(VSIFCloseL(fp));

    return GDALDataset::FromHandle(GDALOpen(pszFilename, GA_Update));
}

/************************************************************************/
/*                          GDALRegister_GTX()                          */
/************************************************************************/

void GDALRegister_GTX()

{
    if (GDALGetDriverByName("GTX") != nullptr)
        return;

    GDALDriver *poDriver = new GDALDriver();

    poDriver->SetDescription("GTX");
    poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "NOAA Vertical Datum .GTX");
    poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "gtx");
    poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    // poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
    //                            "frmt_various.html#GTX" );
    poDriver->SetMetadataItem(
        GDAL_DMD_OPENOPTIONLIST,
        "<OpenOptionList>"
        "   <Option name='SHIFT_ORIGIN_IN_MINUS_180_PLUS_180' type='boolean' "
        "description='Whether to apply a +/-360 deg shift to the longitude of "
        "the top left corner so that it is in the [-180,180] range' "
        "default='NO'/>"
        "</OpenOptionList>");

    poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Float32");

    poDriver->pfnOpen = GTXDataset::Open;
    poDriver->pfnIdentify = GTXDataset::Identify;
    poDriver->pfnCreate = GTXDataset::Create;

    GetGDALDriverManager()->RegisterDriver(poDriver);
}