/******************************************************************************
 *
 * Project:  GDAL
 * Purpose:  Implementation of NOAA .b format used for GEOCON / NADCON5 grids
 * Author:   Even Rouault <even dot rouault at spatialys.com>
 *
 ******************************************************************************
 * Copyright (c) 2022, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

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

#include <limits>

// Specification of the format is at "paragraph 10.2 ".b" grids (GEOCON and
// NADCON 5.0)" of "NOAA Technical Report NOS NGS 63" at
// https://geodesy.noaa.gov/library/pdfs/NOAA_TR_NOS_NGS_0063.pdf

constexpr int HEADER_SIZE = 52;
constexpr int FORTRAN_HEADER_SIZE = 4;
constexpr int FORTRAN_TRAILER_SIZE = 4;

/************************************************************************/
/* ==================================================================== */
/*                          NOAA_B_Dataset                              */
/* ==================================================================== */
/************************************************************************/

class NOAA_B_Dataset final : public RawDataset
{
    OGRSpatialReference m_oSRS{};
    GDALGeoTransform m_gt{};

    CPL_DISALLOW_COPY_ASSIGN(NOAA_B_Dataset)

    static int IdentifyEx(GDALOpenInfo *poOpenInfo, bool &bBigEndianOut);

    CPLErr Close() override
    {
        return GDALPamDataset::Close();
    }

  public:
    NOAA_B_Dataset()
    {
        m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    }

    CPLErr GetGeoTransform(GDALGeoTransform &gt) const override;

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

    static GDALDataset *Open(GDALOpenInfo *);
    static int Identify(GDALOpenInfo *);
};

/************************************************************************/
/* ==================================================================== */
/*                          NOAA_B_Dataset                              */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                           GetHeaderValues()                          */
/************************************************************************/

static void GetHeaderValues(const GDALOpenInfo *poOpenInfo, double &dfSWLat,
                            double &dfSWLon, double &dfDeltaLat,
                            double &dfDeltaLon, int32_t &nRows, int32_t &nCols,
                            int32_t &iKind, bool bBigEndian)
{
    const auto ReadFloat64 = [bBigEndian](const GByte *&ptr)
    {
        double v;
        memcpy(&v, ptr, sizeof(v));
        ptr += sizeof(v);
        if (bBigEndian)
            CPL_MSBPTR64(&v);
        else
            CPL_LSBPTR64(&v);
        return v;
    };

    const auto ReadInt32 = [bBigEndian](const GByte *&ptr)
    {
        int32_t v;
        memcpy(&v, ptr, sizeof(v));
        ptr += sizeof(v);
        if (bBigEndian)
            CPL_MSBPTR32(&v);
        else
            CPL_LSBPTR32(&v);
        return v;
    };

    const GByte *ptr = poOpenInfo->pabyHeader + FORTRAN_HEADER_SIZE;

    dfSWLat = ReadFloat64(ptr);
    dfSWLon = ReadFloat64(ptr);
    dfDeltaLat = ReadFloat64(ptr);
    dfDeltaLon = ReadFloat64(ptr);

    nRows = ReadInt32(ptr);
    nCols = ReadInt32(ptr);
    iKind = ReadInt32(ptr);
}

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

int NOAA_B_Dataset::IdentifyEx(GDALOpenInfo *poOpenInfo, bool &bBigEndianOut)

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

#if !defined(FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION)
    if (!poOpenInfo->IsExtensionEqualToCI("b"))
        return FALSE;
#endif

    // Sanity checks on header
    double dfSWLat;
    double dfSWLon;
    double dfDeltaLat;
    double dfDeltaLon;
    int32_t nRows;
    int32_t nCols;
    int32_t iKind;

    // Fun... nadcon5 files are encoded in big-endian, but vertcon3 files...
    // in little-endian. We could probably figure that out directly from the
    // 4 bytes which are 0x00 0x00 0x00 0x2C for nadcon5, and the reverse for
    // vertcon3, but the semantics of those 4 bytes is undocumented.
    // So try both possibilities and rely on sanity checks.
    for (int i = 0; i < 2; ++i)
    {
        const bool bBigEndian = i == 0 ? true : false;
        GetHeaderValues(poOpenInfo, dfSWLat, dfSWLon, dfDeltaLat, dfDeltaLon,
                        nRows, nCols, iKind, bBigEndian);
        if (!(fabs(dfSWLat) <= 90))
            continue;
        if (!(fabs(dfSWLon) <=
              360))  // NADCON5 grids typically have SWLon > 180
            continue;
        if (!(dfDeltaLat > 0 && dfDeltaLat <= 1))
            continue;
        if (!(dfDeltaLon > 0 && dfDeltaLon <= 1))
            continue;
        if (!(nRows > 0 && dfSWLat + (nRows - 1) * dfDeltaLat <= 90))
            continue;
        if (!(nCols > 0 && (nCols - 1) * dfDeltaLon <= 360))
            continue;
        if (!(iKind >= -1 && iKind <= 2))
            continue;

        bBigEndianOut = bBigEndian;
        return TRUE;
    }
    return FALSE;
}

int NOAA_B_Dataset::Identify(GDALOpenInfo *poOpenInfo)

{
    bool bBigEndian = false;
    return IdentifyEx(poOpenInfo, bBigEndian);
}

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

CPLErr NOAA_B_Dataset::GetGeoTransform(GDALGeoTransform &gt) const

{
    gt = m_gt;
    return CE_None;
}

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

GDALDataset *NOAA_B_Dataset::Open(GDALOpenInfo *poOpenInfo)

{
    bool bBigEndian = false;
    if (!IdentifyEx(poOpenInfo, bBigEndian) || poOpenInfo->fpL == nullptr ||
        poOpenInfo->eAccess == GA_Update)
    {
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*      Read the header.                                                */
    /* -------------------------------------------------------------------- */
    double dfSWLat;
    double dfSWLon;
    double dfDeltaLat;
    double dfDeltaLon;
    int32_t nRows;
    int32_t nCols;
    int32_t iKind;
    GetHeaderValues(poOpenInfo, dfSWLat, dfSWLon, dfDeltaLat, dfDeltaLon, nRows,
                    nCols, iKind, bBigEndian);

    if (iKind == -1)
    {
        CPLError(CE_Failure, CPLE_NotSupported,
                 "KIND = -1 in NOAA .b dataset not supported");
        return nullptr;
    }

    const GDALDataType eDT =
        // iKind == -1 ? GDT_Int16 :
        iKind == 0   ? GDT_Int32
        : iKind == 1 ? GDT_Float32
                     : GDT_Int16;
    const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
    if (!GDALCheckDatasetDimensions(nCols, nRows) ||
        (nDTSize > 0 && static_cast<vsi_l_offset>(nCols) * nRows >
                            std::numeric_limits<vsi_l_offset>::max() / nDTSize))
    {
        return nullptr;
    }
    if (nDTSize > 0 && nCols > (std::numeric_limits<int>::max() -
                                FORTRAN_HEADER_SIZE - FORTRAN_TRAILER_SIZE) /
                                   nDTSize)
    {
        return nullptr;
    }
    const int nLineSize =
        FORTRAN_HEADER_SIZE + nCols * nDTSize + FORTRAN_TRAILER_SIZE;

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

    poDS->nRasterXSize = nCols;
    poDS->nRasterYSize = nRows;

    // Adjust longitude > 180 to [-180, 180] range
    if (dfSWLon > 180)
        dfSWLon -= 360;

    // Convert from south-west center-of-pixel convention to
    // north-east pixel-corner convention
    poDS->m_gt[0] = dfSWLon - dfDeltaLon / 2;
    poDS->m_gt[1] = dfDeltaLon;
    poDS->m_gt[2] = 0.0;
    poDS->m_gt[3] = dfSWLat + (nRows - 1) * dfDeltaLat + dfDeltaLat / 2;
    poDS->m_gt[4] = 0.0;
    poDS->m_gt[5] = -dfDeltaLat;

    /* -------------------------------------------------------------------- */
    /*      Create band information object.                                 */
    /* -------------------------------------------------------------------- */

    // Borrow file handle
    VSILFILE *fpImage = poOpenInfo->fpL;
    poOpenInfo->fpL = nullptr;

    // Records are presented from the southern-most to the northern-most
    auto poBand = RawRasterBand::Create(
        poDS.get(), 1, fpImage,
        // skip to beginning of northern-most line
        HEADER_SIZE +
            static_cast<vsi_l_offset>(poDS->nRasterYSize - 1) * nLineSize +
            FORTRAN_HEADER_SIZE,
        nDTSize, -nLineSize, eDT,
        bBigEndian ? RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN
                   : RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
        RawRasterBand::OwnFP::YES);
    if (!poBand)
        return nullptr;
    poDS->SetBand(1, std::move(poBand));

    /* -------------------------------------------------------------------- */
    /*      Guess CRS from filename.                                        */
    /* -------------------------------------------------------------------- */
    const std::string osFilename(CPLGetFilename(poOpenInfo->pszFilename));

    static const struct
    {
        const char *pszPrefix;
        int nEPSGCode;
    }
    // Cf https://geodesy.noaa.gov/pub/nadcon5/20160901release/Builds/
    asFilenameToCRS[] = {
        {"nadcon5.nad27.", 4267},       // NAD27
        {"nadcon5.pr40.", 4139},        // Puerto Rico (1940)
        {"nadcon5.ohd.", 4135},         // Old Hawaian
        {"nadcon5.sl1952.", 4136},      // Saint Lawrence Island (1952)
        {"nadcon5.sp1952.", 4137},      // Saint Paul Island (1952)
        {"nadcon5.sg1952.", 4138},      // Saint George Island (1952)
        {"nadcon5.as62.", 4169},        // American Samoa 1962
        {"nadcon5.gu63.", 4675},        // Guam 1963
        {"nadcon5.nad83_1986.", 4269},  // NAD83
        {"nadcon5.nad83_harn.", 4152},  // NAD83(HARN)
        {"nadcon5.nad83_1992.",
         4152},  // NAD83(1992) for Alaska is NAD83(HARN) in EPSG
        {"nadcon5.nad83_1993.",
         4152},  // NAD83(1993) for American Samoa, PRVI, Guam and Hawaii is
                 // NAD83(HARN) in EPSG
        {"nadcon5.nad83_1997.", 8545},  // NAD83(HARN Corrected)
        {"nadcon5.nad83_fbn.", 8860},   // NAD83(FBN)
        {"nadcon5.nad83_2002.",
         8860},  // NAD83(2002) for Alaska, PRVI and Guam is NAD83(FBN) in EPSG
        {"nadcon5.nad83_2007.", 4759},  // NAD83(NSRS2007)
    };

    for (const auto &sPair : asFilenameToCRS)
    {
        if (STARTS_WITH_CI(osFilename.c_str(), sPair.pszPrefix))
        {
            poDS->m_oSRS.importFromEPSG(sPair.nEPSGCode);
            break;
        }
    }
    if (poDS->m_oSRS.IsEmpty())
    {
        poDS->m_oSRS.importFromWkt(
            "GEOGCRS[\"Unspecified geographic CRS\",DATUM[\"Unspecified datum "
            "based on GRS80 ellipsoid\",ELLIPSOID[\"GRS "
            "1980\",6378137,298.257222101]],CS[ellipsoidal,2],AXIS[\"geodetic "
            "latitude (Lat)\",north,ANGLEUNIT[\"degree\",0.0174532925199433]], "
            "       AXIS[\"geodetic longitude "
            "(Lon)\",east,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433]]]");
    }

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

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

    return poDS.release();
}

/************************************************************************/
/*                        GDALRegister_NOAA_B()                         */
/************************************************************************/

void GDALRegister_NOAA_B()
{
    if (GDALGetDriverByName("NOAA_B") != nullptr)
        return;

    GDALDriver *poDriver = new GDALDriver();

    poDriver->SetDescription("NOAA_B");
    poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
                              "NOAA GEOCON/NADCON5 .b format");
    poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "b");
    poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/noaa_b.html");

    poDriver->pfnIdentify = NOAA_B_Dataset::Identify;
    poDriver->pfnOpen = NOAA_B_Dataset::Open;

    GetGDALDriverManager()->RegisterDriver(poDriver);
}