/******************************************************************************
 *
 * Project:  ESRI .hdr Driver
 * Purpose:  Implementation of EHdrDataset
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 1999, Frank Warmerdam <warmerdam@pobox.com>
 * Copyright (c) 2007-2013, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_port.h"
#include "ehdrdataset.h"
#include "gdal_priv.h"
#include "rawdataset.h"

#include <algorithm>
#include <cctype>
#include <cerrno>
#include <climits>
#include <cmath>
#include <cstddef>
#include <cstdio>
#include <cstdlib>
#include <cstring>

#include <limits>

#include "cpl_conv.h"
#include "cpl_error.h"
#include "cpl_progress.h"
#include "cpl_string.h"
#include "cpl_vsi.h"
#include "gdal.h"
#include "gdal_frmts.h"
#include "gdal_pam.h"
#include "gdal_priv.h"
#include "ogr_core.h"
#include "ogr_spatialref.h"

constexpr int HAS_MIN_FLAG = 0x1;
constexpr int HAS_MAX_FLAG = 0x2;
constexpr int HAS_MEAN_FLAG = 0x4;
constexpr int HAS_STDDEV_FLAG = 0x8;
constexpr int HAS_ALL_FLAGS =
    HAS_MIN_FLAG | HAS_MAX_FLAG | HAS_MEAN_FLAG | HAS_STDDEV_FLAG;

/************************************************************************/
/*                           EHdrRasterBand()                           */
/************************************************************************/

EHdrRasterBand::EHdrRasterBand(GDALDataset *poDSIn, int nBandIn,
                               VSILFILE *fpRawIn, vsi_l_offset nImgOffsetIn,
                               int nPixelOffsetIn, int nLineOffsetIn,
                               GDALDataType eDataTypeIn,
                               RawRasterBand::ByteOrder eByteOrderIn,
                               int nBitsIn)
    : RawRasterBand(poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
                    nLineOffsetIn, eDataTypeIn, eByteOrderIn,
                    RawRasterBand::OwnFP::NO),
      nBits(nBitsIn), nStartBit(0), nPixelOffsetBits(0), nLineOffsetBits(0),
      bNoDataSet(FALSE), dfNoData(0.0), dfMin(0.0), dfMax(0.0), dfMean(0.0),
      dfStdDev(0.0), minmaxmeanstddev(0)
{
    m_bValid = RawRasterBand::IsValid();

    EHdrDataset *poEDS = cpl::down_cast<EHdrDataset *>(poDS);

    if (nBits < 8)
    {
        const int nSkipBytes = atoi(poEDS->GetKeyValue("SKIPBYTES"));
        if (nSkipBytes < 0 || nSkipBytes > std::numeric_limits<int>::max() / 8)
        {
            m_bValid = false;
            CPLError(CE_Failure, CPLE_AppDefined, "Invalid SKIPBYTES: %d",
                     nSkipBytes);
            nStartBit = 0;
        }
        else
        {
            nStartBit = static_cast<vsi_l_offset>(nSkipBytes) * 8;
        }
        if (nBand >= 2)
        {
            GIntBig nBandRowBytes =
                CPLAtoGIntBig(poEDS->GetKeyValue("BANDROWBYTES"));
            if (nBandRowBytes < 0)
            {
                m_bValid = false;
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Invalid BANDROWBYTES: " CPL_FRMT_GIB, nBandRowBytes);
                nBandRowBytes = 0;
            }
            vsi_l_offset nRowBytes = 0;
            if (nBandRowBytes == 0)
                nRowBytes =
                    (static_cast<vsi_l_offset>(nBits) * poDS->GetRasterXSize() +
                     7) /
                    8;
            else
                nRowBytes = static_cast<vsi_l_offset>(nBandRowBytes);

            nStartBit += nRowBytes * (nBand - 1) * 8;
        }

        nPixelOffsetBits = nBits;
        GIntBig nTotalRowBytes =
            CPLAtoGIntBig(poEDS->GetKeyValue("TOTALROWBYTES"));
        if (nTotalRowBytes < 0 ||
            nTotalRowBytes > GINTBIG_MAX / 8 / poDS->GetRasterYSize())
        {
            m_bValid = false;
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Invalid TOTALROWBYTES: " CPL_FRMT_GIB, nTotalRowBytes);
            nTotalRowBytes = 0;
        }
        if (nTotalRowBytes > 0)
            nLineOffsetBits = static_cast<vsi_l_offset>(nTotalRowBytes * 8);
        else
            nLineOffsetBits = static_cast<vsi_l_offset>(nPixelOffsetBits) *
                              poDS->GetRasterXSize();

        nBlockXSize = poDS->GetRasterXSize();
        nBlockYSize = 1;

        SetMetadataItem("NBITS", CPLString().Printf("%d", nBits),
                        "IMAGE_STRUCTURE");
    }
}

/************************************************************************/
/*                             IReadBlock()                             */
/************************************************************************/

CPLErr EHdrRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)

{
    if (nBits >= 8)
        return RawRasterBand::IReadBlock(nBlockXOff, nBlockYOff, pImage);

    // Establish desired position.
    const vsi_l_offset nLineStart =
        (nStartBit + nLineOffsetBits * nBlockYOff) / 8;
    int iBitOffset =
        static_cast<int>((nStartBit + nLineOffsetBits * nBlockYOff) % 8);
    const vsi_l_offset nLineEnd =
        (nStartBit + nLineOffsetBits * nBlockYOff +
         static_cast<vsi_l_offset>(nPixelOffsetBits) * nBlockXSize - 1) /
        8;
    const vsi_l_offset nLineBytesBig = nLineEnd - nLineStart + 1;
    if (nLineBytesBig >
        static_cast<vsi_l_offset>(std::numeric_limits<int>::max()))
        return CE_Failure;
    const unsigned int nLineBytes = static_cast<unsigned int>(nLineBytesBig);

    // Read data into buffer.
    GByte *pabyBuffer = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nLineBytes));
    if (pabyBuffer == nullptr)
        return CE_Failure;

    if (VSIFSeekL(GetFPL(), nLineStart, SEEK_SET) != 0 ||
        VSIFReadL(pabyBuffer, 1, nLineBytes, GetFPL()) != nLineBytes)
    {
        CPLError(CE_Failure, CPLE_FileIO,
                 "Failed to read %u bytes at offset %lu.\n%s", nLineBytes,
                 static_cast<unsigned long>(nLineStart), VSIStrerror(errno));
        CPLFree(pabyBuffer);
        return CE_Failure;
    }

    // Copy data, promoting to 8bit.
    for (int iX = 0, iPixel = 0; iX < nBlockXSize; iX++)
    {
        int nOutWord = 0;

        for (int iBit = 0; iBit < nBits; iBit++)
        {
            if (pabyBuffer[iBitOffset >> 3] & (0x80 >> (iBitOffset & 7)))
                nOutWord |= (1 << (nBits - 1 - iBit));
            iBitOffset++;
        }

        iBitOffset = iBitOffset + nPixelOffsetBits - nBits;

        reinterpret_cast<GByte *>(pImage)[iPixel++] =
            static_cast<GByte>(nOutWord);
    }

    CPLFree(pabyBuffer);

    return CE_None;
}

/************************************************************************/
/*                            IWriteBlock()                             */
/************************************************************************/

CPLErr EHdrRasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff, void *pImage)

{
    if (nBits >= 8)
        return RawRasterBand::IWriteBlock(nBlockXOff, nBlockYOff, pImage);

    // Establish desired position.
    const vsi_l_offset nLineStart =
        (nStartBit + nLineOffsetBits * nBlockYOff) / 8;
    int iBitOffset =
        static_cast<int>((nStartBit + nLineOffsetBits * nBlockYOff) % 8);
    const vsi_l_offset nLineEnd =
        (nStartBit + nLineOffsetBits * nBlockYOff +
         static_cast<vsi_l_offset>(nPixelOffsetBits) * nBlockXSize - 1) /
        8;
    const vsi_l_offset nLineBytesBig = nLineEnd - nLineStart + 1;
    if (nLineBytesBig >
        static_cast<vsi_l_offset>(std::numeric_limits<int>::max()))
        return CE_Failure;
    const unsigned int nLineBytes = static_cast<unsigned int>(nLineBytesBig);

    // Read data into buffer.
    GByte *pabyBuffer = static_cast<GByte *>(VSI_CALLOC_VERBOSE(nLineBytes, 1));
    if (pabyBuffer == nullptr)
        return CE_Failure;

    if (VSIFSeekL(GetFPL(), nLineStart, SEEK_SET) != 0)
    {
        CPLError(CE_Failure, CPLE_FileIO,
                 "Failed to read %u bytes at offset %lu.\n%s", nLineBytes,
                 static_cast<unsigned long>(nLineStart), VSIStrerror(errno));
        CPLFree(pabyBuffer);
        return CE_Failure;
    }

    CPL_IGNORE_RET_VAL(VSIFReadL(pabyBuffer, nLineBytes, 1, GetFPL()));

    // Copy data, promoting to 8bit.
    for (int iX = 0, iPixel = 0; iX < nBlockXSize; iX++)
    {
        const int nOutWord = reinterpret_cast<GByte *>(pImage)[iPixel++];

        for (int iBit = 0; iBit < nBits; iBit++)
        {
            if (nOutWord & (1 << (nBits - 1 - iBit)))
                pabyBuffer[iBitOffset >> 3] |= (0x80 >> (iBitOffset & 7));
            else
                pabyBuffer[iBitOffset >> 3] &= ~((0x80 >> (iBitOffset & 7)));

            iBitOffset++;
        }

        iBitOffset = iBitOffset + nPixelOffsetBits - nBits;
    }

    // Write the data back out.
    if (VSIFSeekL(GetFPL(), nLineStart, SEEK_SET) != 0 ||
        VSIFWriteL(pabyBuffer, 1, nLineBytes, GetFPL()) != nLineBytes)
    {
        CPLError(CE_Failure, CPLE_FileIO,
                 "Failed to write %u bytes at offset %lu.\n%s", nLineBytes,
                 static_cast<unsigned long>(nLineStart), VSIStrerror(errno));
        return CE_Failure;
    }

    CPLFree(pabyBuffer);

    return CE_None;
}

/************************************************************************/
/*                             IRasterIO()                              */
/************************************************************************/

CPLErr EHdrRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
                                 int nXSize, int nYSize, void *pData,
                                 int nBufXSize, int nBufYSize,
                                 GDALDataType eBufType, GSpacing nPixelSpace,
                                 GSpacing nLineSpace,
                                 GDALRasterIOExtraArg *psExtraArg)

{
    // Defer to RawRasterBand
    if (nBits >= 8)
        return RawRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
                                        pData, nBufXSize, nBufYSize, eBufType,
                                        nPixelSpace, nLineSpace, psExtraArg);

    // Force use of IReadBlock() and IWriteBlock()
    return GDALRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
                                     pData, nBufXSize, nBufYSize, eBufType,
                                     nPixelSpace, nLineSpace, psExtraArg);
}

/************************************************************************/
/*                              OSR_GDS()                               */
/************************************************************************/

static const char *OSR_GDS(char *pszResult, int nResultLen, char **papszNV,
                           const char *pszField, const char *pszDefaultValue)

{
    if (papszNV == nullptr || papszNV[0] == nullptr)
        return pszDefaultValue;

    int iLine = 0;  // Used after for.
    for (; papszNV[iLine] != nullptr &&
           !EQUALN(papszNV[iLine], pszField, strlen(pszField));
         iLine++)
    {
    }

    if (papszNV[iLine] == nullptr)
        return pszDefaultValue;

    char **papszTokens = CSLTokenizeString(papszNV[iLine]);

    if (CSLCount(papszTokens) > 1)
        strncpy(pszResult, papszTokens[1], nResultLen - 1);
    else
        strncpy(pszResult, pszDefaultValue, nResultLen - 1);
    pszResult[nResultLen - 1] = '\0';

    CSLDestroy(papszTokens);
    return pszResult;
}

/************************************************************************/
/* ==================================================================== */
/*                            EHdrDataset                               */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                            EHdrDataset()                             */
/************************************************************************/

EHdrDataset::EHdrDataset()
    : fpImage(nullptr), osHeaderExt("hdr"), bGotTransform(false),
      bHDRDirty(false), papszHDR(nullptr), bCLRDirty(false)
{
    m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
}

/************************************************************************/
/*                            ~EHdrDataset()                            */
/************************************************************************/

EHdrDataset::~EHdrDataset()

{
    EHdrDataset::Close();
}

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

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

        if (nBands > 0 && GetAccess() == GA_Update)
        {
            int bNoDataSet;
            RawRasterBand *poBand =
                cpl::down_cast<RawRasterBand *>(GetRasterBand(1));

            const double dfNoData = poBand->GetNoDataValue(&bNoDataSet);
            if (bNoDataSet)
            {
                ResetKeyValue("NODATA", CPLString().Printf("%.8g", dfNoData));
            }

            if (bCLRDirty)
                RewriteCLR(poBand);

            if (bHDRDirty)
            {
                if (RewriteHDR() != CE_None)
                    eErr = CE_Failure;
            }
        }

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

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

/************************************************************************/
/*                            GetKeyValue()                             */
/************************************************************************/

const char *EHdrDataset::GetKeyValue(const char *pszKey, const char *pszDefault)

{
    for (int i = 0; papszHDR[i] != nullptr; i++)
    {
        if (EQUALN(pszKey, papszHDR[i], strlen(pszKey)) &&
            isspace(static_cast<unsigned char>(papszHDR[i][strlen(pszKey)])))
        {
            const char *pszValue = papszHDR[i] + strlen(pszKey);
            while (isspace(static_cast<unsigned char>(*pszValue)))
                pszValue++;

            return pszValue;
        }
    }

    return pszDefault;
}

/************************************************************************/
/*                           ResetKeyValue()                            */
/*                                                                      */
/*      Replace or add the keyword with the indicated value in the      */
/*      papszHDR list.                                                  */
/************************************************************************/

void EHdrDataset::ResetKeyValue(const char *pszKey, const char *pszValue)

{
    if (strlen(pszValue) > 65)
    {
        CPLAssert(strlen(pszValue) <= 65);
        return;
    }

    char szNewLine[82] = {'\0'};
    snprintf(szNewLine, sizeof(szNewLine), "%-15s%s", pszKey, pszValue);

    for (int i = CSLCount(papszHDR) - 1; i >= 0; i--)
    {
        if (EQUALN(papszHDR[i], szNewLine, strlen(pszKey) + 1))
        {
            if (strcmp(papszHDR[i], szNewLine) != 0)
            {
                CPLFree(papszHDR[i]);
                papszHDR[i] = CPLStrdup(szNewLine);
                bHDRDirty = true;
            }
            return;
        }
    }

    bHDRDirty = true;
    papszHDR = CSLAddString(papszHDR, szNewLine);
}

/************************************************************************/
/*                           RewriteCLR()                               */
/************************************************************************/

void EHdrDataset::RewriteCLR(GDALRasterBand *poBand) const

{
    CPLString osCLRFilename = CPLResetExtensionSafe(GetDescription(), "clr");
    GDALColorTable *poTable = poBand->GetColorTable();
    GDALRasterAttributeTable *poRAT = poBand->GetDefaultRAT();
    if (poTable || poRAT)
    {
        VSILFILE *fp = VSIFOpenL(osCLRFilename, "wt");
        if (fp != nullptr)
        {
            // Write RAT in priority if both are defined
            if (poRAT)
            {
                for (int iEntry = 0; iEntry < poRAT->GetRowCount(); iEntry++)
                {
                    CPLString oLine;
                    oLine.Printf("%3d %3d %3d %3d\n",
                                 poRAT->GetValueAsInt(iEntry, 0),
                                 poRAT->GetValueAsInt(iEntry, 1),
                                 poRAT->GetValueAsInt(iEntry, 2),
                                 poRAT->GetValueAsInt(iEntry, 3));
                    if (VSIFWriteL(reinterpret_cast<void *>(
                                       const_cast<char *>(oLine.c_str())),
                                   strlen(oLine), 1, fp) != 1)
                    {
                        CPLError(CE_Failure, CPLE_FileIO,
                                 "Error while write color table");
                        CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
                        return;
                    }
                }
            }
            else
            {
                for (int iColor = 0; iColor < poTable->GetColorEntryCount();
                     iColor++)
                {
                    GDALColorEntry sEntry;
                    poTable->GetColorEntryAsRGB(iColor, &sEntry);

                    // I wish we had a way to mark transparency.
                    CPLString oLine;
                    oLine.Printf("%3d %3d %3d %3d\n", iColor, sEntry.c1,
                                 sEntry.c2, sEntry.c3);
                    if (VSIFWriteL(reinterpret_cast<void *>(
                                       const_cast<char *>(oLine.c_str())),
                                   strlen(oLine), 1, fp) != 1)
                    {
                        CPLError(CE_Failure, CPLE_FileIO,
                                 "Error while write color table");
                        CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
                        return;
                    }
                }
            }
            if (VSIFCloseL(fp) != 0)
            {
                CPLError(CE_Failure, CPLE_FileIO,
                         "Error while write color table");
            }
        }
        else
        {
            CPLError(CE_Failure, CPLE_OpenFailed,
                     "Unable to create color file %s.", osCLRFilename.c_str());
        }
    }
    else
    {
        VSIUnlink(osCLRFilename);
    }
}

/************************************************************************/
/*                           SetSpatialRef()                            */
/************************************************************************/

CPLErr EHdrDataset::SetSpatialRef(const OGRSpatialReference *poSRS)

{
    // Reset coordinate system on the dataset.
    m_oSRS.Clear();
    if (poSRS == nullptr)
        return CE_None;

    m_oSRS = *poSRS;
    // Convert to ESRI WKT.
    char *pszESRI_SRS = nullptr;
    const char *const apszOptions[] = {"FORMAT=WKT1_ESRI", nullptr};
    m_oSRS.exportToWkt(&pszESRI_SRS, apszOptions);

    if (pszESRI_SRS)
    {
        // Write to .prj file.
        CPLString osPrjFilename =
            CPLResetExtensionSafe(GetDescription(), "prj");
        VSILFILE *fp = VSIFOpenL(osPrjFilename.c_str(), "wt");
        if (fp != nullptr)
        {
            size_t nCount = VSIFWriteL(pszESRI_SRS, strlen(pszESRI_SRS), 1, fp);
            nCount += VSIFWriteL("\n", 1, 1, fp);
            if (VSIFCloseL(fp) != 0 || nCount != 2)
            {
                CPLFree(pszESRI_SRS);
                return CE_Failure;
            }
        }

        CPLFree(pszESRI_SRS);
    }

    return CE_None;
}

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

CPLErr EHdrDataset::GetGeoTransform(GDALGeoTransform &gt) const

{
    if (bGotTransform)
    {
        gt = m_gt;
        return CE_None;
    }

    return GDALPamDataset::GetGeoTransform(gt);
}

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

CPLErr EHdrDataset::SetGeoTransform(const GDALGeoTransform &gt)

{
    // We only support non-rotated images with info in the .HDR file.
    if (gt[2] != 0.0 || gt[4] != 0.0)
    {
        return GDALPamDataset::SetGeoTransform(gt);
    }

    // Record new geotransform.
    bGotTransform = true;
    m_gt = gt;

    // Strip out all old geotransform keywords from HDR records.
    for (int i = CSLCount(papszHDR) - 1; i >= 0; i--)
    {
        if (STARTS_WITH_CI(papszHDR[i], "ul") ||
            STARTS_WITH_CI(papszHDR[i] + 1, "ll") ||
            STARTS_WITH_CI(papszHDR[i], "cell") ||
            STARTS_WITH_CI(papszHDR[i] + 1, "dim"))
        {
            papszHDR = CSLRemoveStrings(papszHDR, i, 1, nullptr);
        }
    }

    // Set the transformation information.
    CPLString oValue;

    oValue.Printf("%.15g", m_gt[0] + m_gt[1] * 0.5);
    ResetKeyValue("ULXMAP", oValue);

    oValue.Printf("%.15g", m_gt[3] + m_gt[5] * 0.5);
    ResetKeyValue("ULYMAP", oValue);

    oValue.Printf("%.15g", m_gt[1]);
    ResetKeyValue("XDIM", oValue);

    oValue.Printf("%.15g", fabs(m_gt[5]));
    ResetKeyValue("YDIM", oValue);

    return CE_None;
}

/************************************************************************/
/*                             RewriteHDR()                             */
/************************************************************************/

CPLErr EHdrDataset::RewriteHDR()

{
    const CPLString osPath = CPLGetPathSafe(GetDescription());
    const CPLString osName = CPLGetBasenameSafe(GetDescription());
    CPLString osHDRFilename =
        CPLFormCIFilenameSafe(osPath, osName, osHeaderExt);

    // Write .hdr file.
    VSILFILE *fp = VSIFOpenL(osHDRFilename, "wt");

    if (fp == nullptr)
    {
        CPLError(CE_Failure, CPLE_OpenFailed, "Failed to rewrite .hdr file %s.",
                 osHDRFilename.c_str());
        return CE_Failure;
    }

    for (int i = 0; papszHDR[i] != nullptr; i++)
    {
        size_t nCount = VSIFWriteL(papszHDR[i], strlen(papszHDR[i]), 1, fp);
        nCount += VSIFWriteL("\n", 1, 1, fp);
        if (nCount != 2)
        {
            CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
            return CE_Failure;
        }
    }

    bHDRDirty = false;

    if (VSIFCloseL(fp) != 0)
        return CE_Failure;

    return CE_None;
}

/************************************************************************/
/*                             RewriteSTX()                             */
/************************************************************************/

CPLErr EHdrDataset::RewriteSTX() const
{
    const CPLString osPath = CPLGetPathSafe(GetDescription());
    const CPLString osName = CPLGetBasenameSafe(GetDescription());
    const CPLString osSTXFilename =
        CPLFormCIFilenameSafe(osPath, osName, "stx");

    VSILFILE *fp = VSIFOpenL(osSTXFilename, "wt");
    if (fp == nullptr)
    {
        CPLDebug("EHDR", "Failed to rewrite .stx file %s.",
                 osSTXFilename.c_str());
        return CE_Failure;
    }

    bool bOK = true;
    for (int i = 0; bOK && i < nBands; ++i)
    {
        EHdrRasterBand *poBand =
            reinterpret_cast<EHdrRasterBand *>(papoBands[i]);
        bOK &= VSIFPrintfL(fp, "%d %.10f %.10f ", i + 1, poBand->dfMin,
                           poBand->dfMax) >= 0;
        if (poBand->minmaxmeanstddev & HAS_MEAN_FLAG)
            bOK &= VSIFPrintfL(fp, "%.10f ", poBand->dfMean) >= 0;
        else
            bOK &= VSIFPrintfL(fp, "# ") >= 0;

        if (poBand->minmaxmeanstddev & HAS_STDDEV_FLAG)
            bOK &= VSIFPrintfL(fp, "%.10f\n", poBand->dfStdDev) >= 0;
        else
            bOK &= VSIFPrintfL(fp, "#\n") >= 0;
    }

    if (VSIFCloseL(fp) != 0)
        bOK = false;

    return bOK ? CE_None : CE_Failure;
}

/************************************************************************/
/*                              ReadSTX()                               */
/************************************************************************/

CPLErr EHdrDataset::ReadSTX() const
{
    const CPLString osPath = CPLGetPathSafe(GetDescription());
    const CPLString osName = CPLGetBasenameSafe(GetDescription());
    const CPLString osSTXFilename =
        CPLFormCIFilenameSafe(osPath, osName, "stx");

    VSILFILE *fp = VSIFOpenL(osSTXFilename, "rt");
    if (fp == nullptr)
        return CE_None;

    const char *pszLine = nullptr;
    while ((pszLine = CPLReadLineL(fp)) != nullptr)
    {
        char **papszTokens =
            CSLTokenizeStringComplex(pszLine, " \t", TRUE, FALSE);
        const int nTokens = CSLCount(papszTokens);
        if (nTokens >= 5)
        {
            const int i = atoi(papszTokens[0]);
            if (i > 0 && i <= nBands)
            {
                EHdrRasterBand *poBand =
                    reinterpret_cast<EHdrRasterBand *>(papoBands[i - 1]);
                poBand->dfMin = CPLAtof(papszTokens[1]);
                poBand->dfMax = CPLAtof(papszTokens[2]);

                int bNoDataSet = FALSE;
                const double dfNoData = poBand->GetNoDataValue(&bNoDataSet);
                if (bNoDataSet && dfNoData == poBand->dfMin)
                {
                    // Triggered by
                    // /vsicurl/http://eros.usgs.gov/archive/nslrsda/GeoTowns/HongKong/srtm/n22e113.zip/n22e113.bil
                    CPLDebug(
                        "EHDr",
                        "Ignoring .stx file where min == nodata. "
                        "The nodata value should not be taken into account "
                        "in minimum value computation.");
                    CSLDestroy(papszTokens);
                    papszTokens = nullptr;
                    break;
                }

                poBand->minmaxmeanstddev = HAS_MIN_FLAG | HAS_MAX_FLAG;
                // Reads optional mean and stddev.
                if (!EQUAL(papszTokens[3], "#"))
                {
                    poBand->dfMean = CPLAtof(papszTokens[3]);
                    poBand->minmaxmeanstddev |= HAS_MEAN_FLAG;
                }
                if (!EQUAL(papszTokens[4], "#"))
                {
                    poBand->dfStdDev = CPLAtof(papszTokens[4]);
                    poBand->minmaxmeanstddev |= HAS_STDDEV_FLAG;
                }

                if (nTokens >= 6 && !EQUAL(papszTokens[5], "#"))
                    poBand->SetMetadataItem("STRETCHMIN", papszTokens[5],
                                            "RENDERING_HINTS");

                if (nTokens >= 7 && !EQUAL(papszTokens[6], "#"))
                    poBand->SetMetadataItem("STRETCHMAX", papszTokens[6],
                                            "RENDERING_HINTS");
            }
        }

        CSLDestroy(papszTokens);
    }

    CPL_IGNORE_RET_VAL(VSIFCloseL(fp));

    return CE_None;
}

/************************************************************************/
/*                      GetImageRepFilename()                           */
/************************************************************************/

// Check for IMAGE.REP (Spatiocarte Defense 1.0) or name_of_image.rep
// if it is a GIS-GeoSPOT image.
// For the specification of SPDF (in French), see
//   http://eden.ign.fr/download/pub/doc/emabgi/spdf10.pdf/download

CPLString EHdrDataset::GetImageRepFilename(const char *pszFilename)
{

    const CPLString osPath = CPLGetPathSafe(pszFilename);
    const CPLString osName = CPLGetBasenameSafe(pszFilename);
    CPLString osREPFilename = CPLFormCIFilenameSafe(osPath, osName, "rep");

    VSIStatBufL sStatBuf;
    if (VSIStatExL(osREPFilename.c_str(), &sStatBuf, VSI_STAT_EXISTS_FLAG) == 0)
        return osREPFilename;

    if (EQUAL(CPLGetFilename(pszFilename), "imspatio.bil") ||
        EQUAL(CPLGetFilename(pszFilename), "haspatio.bil"))
    {
        CPLString osImageRepFilename(
            CPLFormCIFilenameSafe(osPath, "image", "rep"));
        if (VSIStatExL(osImageRepFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG) ==
            0)
            return osImageRepFilename;

        // Try in the upper directories if not found in the BIL image directory.
        CPLString dirName(CPLGetDirnameSafe(osPath));
        if (CPLIsFilenameRelative(osPath.c_str()))
        {
            char *cwd = CPLGetCurrentDir();
            if (cwd)
            {
                dirName = CPLFormFilenameSafe(cwd, dirName.c_str(), nullptr);
                CPLFree(cwd);
            }
        }
        while (dirName[0] != 0 && EQUAL(dirName, ".") == FALSE &&
               EQUAL(dirName, "/") == FALSE)
        {
            osImageRepFilename =
                CPLFormCIFilenameSafe(dirName.c_str(), "image", "rep");
            if (VSIStatExL(osImageRepFilename.c_str(), &sStatBuf,
                           VSI_STAT_EXISTS_FLAG) == 0)
                return osImageRepFilename;

            // Don't try to recurse above the 'image' subdirectory.
            if (EQUAL(dirName, "image"))
            {
                break;
            }
            dirName = CPLString(CPLGetDirnameSafe(dirName));
        }
    }
    return CPLString();
}

/************************************************************************/
/*                            GetFileList()                             */
/************************************************************************/

char **EHdrDataset::GetFileList()

{
    const CPLString osPath = CPLGetPathSafe(GetDescription());
    const CPLString osName = CPLGetBasenameSafe(GetDescription());

    // Main data file, etc.
    char **papszFileList = GDALPamDataset::GetFileList();

    // Header file.
    CPLString osFilename = CPLFormCIFilenameSafe(osPath, osName, osHeaderExt);
    papszFileList = CSLAddString(papszFileList, osFilename);

    // Statistics file
    osFilename = CPLFormCIFilenameSafe(osPath, osName, "stx");
    VSIStatBufL sStatBuf;
    if (VSIStatExL(osFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG) == 0)
        papszFileList = CSLAddString(papszFileList, osFilename);

    // color table file.
    osFilename = CPLFormCIFilenameSafe(osPath, osName, "clr");
    if (VSIStatExL(osFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG) == 0)
        papszFileList = CSLAddString(papszFileList, osFilename);

    // projections file.
    osFilename = CPLFormCIFilenameSafe(osPath, osName, "prj");
    if (VSIStatExL(osFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG) == 0)
        papszFileList = CSLAddString(papszFileList, osFilename);

    const CPLString imageRepFilename = GetImageRepFilename(GetDescription());
    if (!imageRepFilename.empty())
        papszFileList = CSLAddString(papszFileList, imageRepFilename.c_str());

    return papszFileList;
}

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

GDALDataset *EHdrDataset::Open(GDALOpenInfo *poOpenInfo)

{
    return Open(poOpenInfo, true);
}

GDALDataset *EHdrDataset::Open(GDALOpenInfo *poOpenInfo, bool bFileSizeCheck)

{
    // Assume the caller is pointing to the binary (i.e. .bil) file.
    if (poOpenInfo->nHeaderBytes < 2 || poOpenInfo->fpL == nullptr ||
        (!poOpenInfo->IsSingleAllowedDriver("EHdr") &&
         poOpenInfo->IsExtensionEqualToCI("zarr")))
    {
        return nullptr;
    }

    // Tear apart the filename to form a .HDR filename.
    const CPLString osPath = CPLGetPathSafe(poOpenInfo->pszFilename);
    const CPLString osName = CPLGetBasenameSafe(poOpenInfo->pszFilename);

    const char *pszHeaderExt = "hdr";
    if (poOpenInfo->IsExtensionEqualToCI("SRC") && osName.size() == 7 &&
        (osName[0] == 'e' || osName[0] == 'E' || osName[0] == 'w' ||
         osName[0] == 'W') &&
        (osName[4] == 'n' || osName[4] == 'N' || osName[4] == 's' ||
         osName[4] == 'S'))
    {
        // It is a GTOPO30 or SRTM30 source file, whose header extension is .sch
        // see http://dds.cr.usgs.gov/srtm/version1/SRTM30/GTOPO30_Documentation
        pszHeaderExt = "sch";
    }

    CSLConstList papszSiblingFiles = poOpenInfo->GetSiblingFiles();
    CPLString osHDRFilename;
    if (papszSiblingFiles)
    {
        const int iFile = CSLFindString(
            papszSiblingFiles,
            CPLFormFilenameSafe(nullptr, osName, pszHeaderExt).c_str());
        if (iFile < 0)  // Return if there is no corresponding .hdr file.
            return nullptr;

        osHDRFilename =
            CPLFormFilenameSafe(osPath, papszSiblingFiles[iFile], nullptr);
    }
    else
    {
        osHDRFilename = CPLFormCIFilenameSafe(osPath, osName, pszHeaderExt);
    }

    const bool bSelectedHDR = EQUAL(osHDRFilename, poOpenInfo->pszFilename);

    // Do we have a .hdr file?
    VSILFILE *fp = VSIFOpenL(osHDRFilename, "r");
    if (fp == nullptr)
    {
        return nullptr;
    }

    // Is this file an ESRI header file?  Read a few lines of text
    // searching for something starting with nrows or ncols.
    int nRows = -1;
    int nCols = -1;
    int l_nBands = 1;
    int nSkipBytes = 0;
    double dfULXMap = 0.5;
    double dfULYMap = 0.5;
    double dfYLLCorner = -123.456;
    int bCenter = TRUE;
    double dfXDim = 1.0;
    double dfYDim = 1.0;
    double dfNoData = 0.0;
    int nLineCount = 0;
    int bNoDataSet = FALSE;
    GDALDataType eDataType = GDT_Byte;
    int nBits = -1;
    char chByteOrder = 'M';
    char chPixelType = 'N';  // Not defined.
    char szLayout[10] = "BIL";
    char **papszHDR = nullptr;
    int bHasInternalProjection = FALSE;
    int bHasMin = FALSE;
    int bHasMax = FALSE;
    double dfMin = 0;
    double dfMax = 0;

    const char *pszLine = nullptr;
    while ((pszLine = CPLReadLineL(fp)) != nullptr)
    {
        nLineCount++;

        if (nLineCount > 50 || strlen(pszLine) > 1000)
            break;

        papszHDR = CSLAddString(papszHDR, pszLine);

        char **papszTokens =
            CSLTokenizeStringComplex(pszLine, " \t", TRUE, FALSE);
        if (CSLCount(papszTokens) < 2)
        {
            CSLDestroy(papszTokens);
            continue;
        }

        if (EQUAL(papszTokens[0], "ncols"))
        {
            nCols = atoi(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "nrows"))
        {
            nRows = atoi(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "skipbytes"))
        {
            nSkipBytes = atoi(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "ulxmap") ||
                 EQUAL(papszTokens[0], "xllcorner") ||
                 EQUAL(papszTokens[0], "xllcenter"))
        {
            dfULXMap = CPLAtofM(papszTokens[1]);
            if (EQUAL(papszTokens[0], "xllcorner"))
                bCenter = FALSE;
        }
        else if (EQUAL(papszTokens[0], "ulymap"))
        {
            dfULYMap = CPLAtofM(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "yllcorner") ||
                 EQUAL(papszTokens[0], "yllcenter"))
        {
            dfYLLCorner = CPLAtofM(papszTokens[1]);
            if (EQUAL(papszTokens[0], "yllcorner"))
                bCenter = FALSE;
        }
        else if (EQUAL(papszTokens[0], "xdim"))
        {
            dfXDim = CPLAtofM(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "ydim"))
        {
            dfYDim = CPLAtofM(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "cellsize"))
        {
            dfXDim = CPLAtofM(papszTokens[1]);
            dfYDim = dfXDim;
        }
        else if (EQUAL(papszTokens[0], "nbands"))
        {
            l_nBands = atoi(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "layout"))
        {
            snprintf(szLayout, sizeof(szLayout), "%s", papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "NODATA_value") ||
                 EQUAL(papszTokens[0], "NODATA"))
        {
            dfNoData = CPLAtofM(papszTokens[1]);
            bNoDataSet = TRUE;
        }
        else if (EQUAL(papszTokens[0], "NBITS"))
        {
            nBits = atoi(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "PIXELTYPE"))
        {
            chPixelType = static_cast<char>(
                toupper(static_cast<unsigned char>(papszTokens[1][0])));
        }
        else if (EQUAL(papszTokens[0], "byteorder"))
        {
            chByteOrder = static_cast<char>(
                toupper(static_cast<unsigned char>(papszTokens[1][0])));
        }

        // http://www.worldclim.org/futdown.htm have the projection extensions
        else if (EQUAL(papszTokens[0], "Projection"))
        {
            bHasInternalProjection = TRUE;
        }
        else if (EQUAL(papszTokens[0], "MinValue") ||
                 EQUAL(papszTokens[0], "MIN_VALUE"))
        {
            dfMin = CPLAtofM(papszTokens[1]);
            bHasMin = TRUE;
        }
        else if (EQUAL(papszTokens[0], "MaxValue") ||
                 EQUAL(papszTokens[0], "MAX_VALUE"))
        {
            dfMax = CPLAtofM(papszTokens[1]);
            bHasMax = TRUE;
        }

        CSLDestroy(papszTokens);
    }

    CPL_IGNORE_RET_VAL(VSIFCloseL(fp));

    // Did we get the required keywords?  If not, return with this never having
    // been considered to be a match. This isn't an error!
    if (nRows == -1 || nCols == -1)
    {
        CSLDestroy(papszHDR);
        return nullptr;
    }

    if (!GDALCheckDatasetDimensions(nCols, nRows) ||
        !GDALCheckBandCount(l_nBands, FALSE))
    {
        CSLDestroy(papszHDR);
        return nullptr;
    }

    // Has the caller selected the .hdr file to open?
    if (bSelectedHDR)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "The selected file is an ESRI BIL header file, but to "
                 "open ESRI BIL datasets, the data file should be selected "
                 "instead of the .hdr file.  Please try again selecting "
                 "the data file (often with the extension .bil) corresponding "
                 "to the header file: %s",
                 poOpenInfo->pszFilename);
        CSLDestroy(papszHDR);
        return nullptr;
    }

    // If we aren't sure of the file type, check the data file size.  If it is 4
    // bytes or more per pixel then we assume it is floating point data.
    if (nBits == -1 && chPixelType == 'N')
    {
        VSIStatBufL sStatBuf;
        if (VSIStatL(poOpenInfo->pszFilename, &sStatBuf) == 0)
        {
            const vsi_l_offset nBytes =
                sStatBuf.st_size / nCols / nRows / l_nBands;
            // Exit now if nBytes value does not make sense to avoid later overflow in below multiplication
            if (nBytes > 8)
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "EHdr driver cannot infer NBITS value");
                CSLDestroy(papszHDR);
                return nullptr;
            }
            if (nBytes > 0 && nBytes != 3)
                nBits = static_cast<int>(nBytes * 8);

            if (nBytes == 4)
                chPixelType = 'F';
        }
    }

    // If the extension is FLT it is likely a floating point file.
    if (chPixelType == 'N')
    {
        if (poOpenInfo->IsExtensionEqualToCI("FLT"))
            chPixelType = 'F';
    }

    // If we have a negative nodata value, assume that the
    // pixel type is signed. This is necessary for datasets from
    // http://www.worldclim.org/futdown.htm

    if (bNoDataSet && dfNoData < 0 && chPixelType == 'N')
    {
        chPixelType = 'S';
    }

    auto poDS = std::make_unique<EHdrDataset>();

    poDS->osHeaderExt = pszHeaderExt;

    poDS->nRasterXSize = nCols;
    poDS->nRasterYSize = nRows;
    poDS->papszHDR = papszHDR;
    std::swap(poDS->fpImage, poOpenInfo->fpL);
    poDS->eAccess = poOpenInfo->eAccess;

    // Figure out the data type.
    if (nBits == 16)
    {
        if (chPixelType == 'S')
            eDataType = GDT_Int16;
        else
            eDataType = GDT_UInt16;  // Default
    }
    else if (nBits == 32)
    {
        if (chPixelType == 'S')
            eDataType = GDT_Int32;
        else if (chPixelType == 'F')
            eDataType = GDT_Float32;
        else
            eDataType = GDT_UInt32;  // Default
    }
    else if (nBits >= 1 && nBits <= 8)
    {
        if (chPixelType == 'S')
            eDataType = GDT_Int8;
        else
            eDataType = GDT_Byte;
        nBits = 8;
    }
    else if (nBits == -1)
    {
        if (chPixelType == 'F')
        {
            eDataType = GDT_Float32;
            nBits = 32;
        }
        else
        {
            eDataType = GDT_Byte;
            nBits = 8;
        }
    }
    else
    {
        CPLError(CE_Failure, CPLE_NotSupported,
                 "EHdr driver does not support %d NBITS value.", nBits);
        return nullptr;
    }

    // Compute the line offset.
    const int nItemSize = GDALGetDataTypeSizeBytes(eDataType);
    CPLAssert(nItemSize != 0);
    CPLAssert(l_nBands != 0);

    int nPixelOffset = 0;
    int nLineOffset = 0;
    vsi_l_offset nBandOffset = 0;

    if (EQUAL(szLayout, "BIP"))
    {
        if (nCols > std::numeric_limits<int>::max() / (nItemSize * l_nBands))
        {
            CPLError(CE_Failure, CPLE_AppDefined, "Int overflow occurred.");
            return nullptr;
        }
        nPixelOffset = nItemSize * l_nBands;
        nLineOffset = nPixelOffset * nCols;
        nBandOffset = static_cast<vsi_l_offset>(nItemSize);
    }
    else if (EQUAL(szLayout, "BSQ"))
    {
        if (nCols > std::numeric_limits<int>::max() / nItemSize)
        {
            CPLError(CE_Failure, CPLE_AppDefined, "Int overflow occurred.");
            return nullptr;
        }
        nPixelOffset = nItemSize;
        nLineOffset = nPixelOffset * nCols;
        nBandOffset = static_cast<vsi_l_offset>(nLineOffset) * nRows;
    }
    else
    {
        // Assume BIL.
        if (nCols > std::numeric_limits<int>::max() / (nItemSize * l_nBands))
        {
            CPLError(CE_Failure, CPLE_AppDefined, "Int overflow occurred.");
            return nullptr;
        }
        nPixelOffset = nItemSize;
        nLineOffset = nItemSize * l_nBands * nCols;
        nBandOffset = static_cast<vsi_l_offset>(nItemSize) * nCols;
    }

    if (nBits >= 8 && bFileSizeCheck &&
        !RAWDatasetCheckMemoryUsage(
            poDS->nRasterXSize, poDS->nRasterYSize, l_nBands, nItemSize,
            nPixelOffset, nLineOffset, nSkipBytes, nBandOffset, poDS->fpImage))
    {
        return nullptr;
    }

    poDS->SetDescription(poOpenInfo->pszFilename);
    poDS->PamInitialize();

    // Create band information objects.
    for (int i = 0; i < l_nBands; i++)
    {
        auto poBand = std::make_unique<EHdrRasterBand>(
            poDS.get(), i + 1, poDS->fpImage, nSkipBytes + nBandOffset * i,
            nPixelOffset, nLineOffset, eDataType,
            chByteOrder == 'I' || chByteOrder == 'L'
                ? RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN
                : RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN,
            nBits);
        if (!poBand->IsValid())
            return nullptr;

        poBand->bNoDataSet = bNoDataSet;
        poBand->dfNoData = dfNoData;

        if (bHasMin && bHasMax)
        {
            poBand->dfMin = dfMin;
            poBand->dfMax = dfMax;
            poBand->minmaxmeanstddev = HAS_MIN_FLAG | HAS_MAX_FLAG;
        }

        poDS->SetBand(i + 1, std::move(poBand));
    }

    // If we didn't get bounds in the .hdr, look for a worldfile.
    if (dfYLLCorner != -123.456)
    {
        if (bCenter)
            dfULYMap = dfYLLCorner + (nRows - 1) * dfYDim;
        else
            dfULYMap = dfYLLCorner + nRows * dfYDim;
    }

    if (dfULXMap != 0.5 || dfULYMap != 0.5 || dfXDim != 1.0 || dfYDim != 1.0)
    {
        poDS->bGotTransform = true;

        if (bCenter)
        {
            poDS->m_gt[0] = dfULXMap - dfXDim * 0.5;
            poDS->m_gt[1] = dfXDim;
            poDS->m_gt[2] = 0.0;
            poDS->m_gt[3] = dfULYMap + dfYDim * 0.5;
            poDS->m_gt[4] = 0.0;
            poDS->m_gt[5] = -dfYDim;
        }
        else
        {
            poDS->m_gt[0] = dfULXMap;
            poDS->m_gt[1] = dfXDim;
            poDS->m_gt[2] = 0.0;
            poDS->m_gt[3] = dfULYMap;
            poDS->m_gt[4] = 0.0;
            poDS->m_gt[5] = -dfYDim;
        }
    }

    if (!poDS->bGotTransform)
        poDS->bGotTransform = CPL_TO_BOOL(GDALReadWorldFile(
            poOpenInfo->pszFilename, nullptr, poDS->m_gt.data()));

    if (!poDS->bGotTransform)
        poDS->bGotTransform = CPL_TO_BOOL(GDALReadWorldFile(
            poOpenInfo->pszFilename, "wld", poDS->m_gt.data()));

    // Check for a .prj file.
    std::string osPrjFilename = CPLFormCIFilenameSafe(osPath, osName, "prj");

    fp = VSIFOpenL(osPrjFilename.c_str(), "r");

    // .hdr files from http://www.worldclim.org/futdown.htm have the projection
    // info in the .hdr file itself.
    if (fp == nullptr && bHasInternalProjection)
    {
        osPrjFilename = std::move(osHDRFilename);
        fp = VSIFOpenL(osPrjFilename.c_str(), "r");
    }

    if (fp != nullptr)
    {
        CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
        fp = nullptr;

        char **papszLines = CSLLoad(osPrjFilename.c_str());

        if (poDS->m_oSRS.importFromESRI(papszLines) == OGRERR_NONE)
        {
            // If geographic values are in seconds, we must transform.
            // Is there a code for minutes too?
            char szResult[80] = {'\0'};
            if (poDS->m_oSRS.IsGeographic() &&
                EQUAL(OSR_GDS(szResult, sizeof(szResult), papszLines, "Units",
                              ""),
                      "DS"))
            {
                poDS->m_gt[0] /= 3600.0;
                poDS->m_gt[1] /= 3600.0;
                poDS->m_gt[2] /= 3600.0;
                poDS->m_gt[3] /= 3600.0;
                poDS->m_gt[4] /= 3600.0;
                poDS->m_gt[5] /= 3600.0;
            }
        }
        else
        {
            poDS->m_oSRS.Clear();
        }

        CSLDestroy(papszLines);
    }
    else
    {
        // Check for IMAGE.REP (Spatiocarte Defense 1.0) or name_of_image.rep
        // if it is a GIS-GeoSPOT image
        // For the specification of SPDF (in French), see
        //   http://eden.ign.fr/download/pub/doc/emabgi/spdf10.pdf/download
        const CPLString szImageRepFilename =
            GetImageRepFilename(poOpenInfo->pszFilename);
        if (!szImageRepFilename.empty())
        {
            fp = VSIFOpenL(szImageRepFilename.c_str(), "r");
        }
        if (fp != nullptr)
        {
            bool bUTM = false;
            bool bWGS84 = false;
            int bNorth = FALSE;
            bool bSouth = false;
            int utmZone = 0;

            while ((pszLine = CPLReadLineL(fp)) != nullptr)
            {
                if (STARTS_WITH(pszLine, "PROJ_ID") && strstr(pszLine, "UTM"))
                {
                    bUTM = true;
                }
                else if (STARTS_WITH(pszLine, "PROJ_ZONE"))
                {
                    const char *c = strchr(pszLine, '"');
                    if (c)
                    {
                        c++;
                        if (*c >= '0' && *c <= '9')
                        {
                            utmZone = atoi(c);
                            if (utmZone >= 1 && utmZone <= 60)
                            {
                                if (strstr(pszLine, "Nord") ||
                                    strstr(pszLine, "NORD"))
                                {
                                    bNorth = TRUE;
                                }
                                else if (strstr(pszLine, "Sud") ||
                                         strstr(pszLine, "SUD"))
                                {
                                    bSouth = true;
                                }
                            }
                        }
                    }
                }
                else if (STARTS_WITH(pszLine, "PROJ_CODE") &&
                         strstr(pszLine, "FR-MINDEF"))
                {
                    const char *c = strchr(pszLine, 'A');
                    if (c)
                    {
                        c++;
                        if (*c >= '0' && *c <= '9')
                        {
                            utmZone = atoi(c);
                            if (utmZone >= 1 && utmZone <= 60)
                            {
                                if (c[1] == 'N' ||
                                    (c[1] != '\0' && c[2] == 'N'))
                                {
                                    bNorth = TRUE;
                                }
                                else if (c[1] == 'S' ||
                                         (c[1] != '\0' && c[2] == 'S'))
                                {
                                    bSouth = true;
                                }
                            }
                        }
                    }
                }
                else if (STARTS_WITH(pszLine, "HORIZ_DATUM") &&
                         (strstr(pszLine, "WGS 84") ||
                          strstr(pszLine, "WGS84")))
                {
                    bWGS84 = true;
                }
                else if (STARTS_WITH(pszLine, "MAP_NUMBER"))
                {
                    const char *c = strchr(pszLine, '"');
                    if (c)
                    {
                        char *pszMapNumber = CPLStrdup(c + 1);
                        char *c2 = strchr(pszMapNumber, '"');
                        if (c2)
                            *c2 = 0;
                        poDS->SetMetadataItem("SPDF_MAP_NUMBER", pszMapNumber);
                        CPLFree(pszMapNumber);
                    }
                }
                else if (STARTS_WITH(pszLine, "PRODUCTION_DATE"))
                {
                    const char *c = pszLine + strlen("PRODUCTION_DATE");
                    while (*c == ' ')
                        c++;
                    if (*c)
                    {
                        poDS->SetMetadataItem("SPDF_PRODUCTION_DATE", c);
                    }
                }
            }

            CPL_IGNORE_RET_VAL(VSIFCloseL(fp));

            if (utmZone >= 1 && utmZone <= 60 && bUTM && bWGS84 &&
                (bNorth || bSouth))
            {
                poDS->m_oSRS.importFromEPSG((bNorth ? 32600 : 32700) + utmZone);
            }
            else
            {
                CPLError(CE_Warning, CPLE_NotSupported,
                         "Cannot retrieve projection from IMAGE.REP");
            }
        }
    }

    // Check for a color table.
    const std::string osCLRFilename =
        CPLFormCIFilenameSafe(osPath, osName, "clr");

    // Only read the .clr for byte, int16 or uint16 bands.
    if (nItemSize <= 2)
        fp = VSIFOpenL(osCLRFilename.c_str(), "r");
    else
        fp = nullptr;

    if (fp != nullptr)
    {
        std::shared_ptr<GDALRasterAttributeTable> poRat(
            new GDALDefaultRasterAttributeTable());
        poRat->CreateColumn("Value", GFT_Integer, GFU_Generic);
        poRat->CreateColumn("Red", GFT_Integer, GFU_Red);
        poRat->CreateColumn("Green", GFT_Integer, GFU_Green);
        poRat->CreateColumn("Blue", GFT_Integer, GFU_Blue);

        poDS->m_poColorTable.reset(new GDALColorTable());

        bool bHasFoundNonCTValues = false;
        int nRatRow = 0;

        while (true)
        {
            pszLine = CPLReadLineL(fp);
            if (!pszLine)
                break;

            if (*pszLine == '#' || *pszLine == '!')
                continue;

            char **papszValues =
                CSLTokenizeString2(pszLine, "\t ", CSLT_HONOURSTRINGS);

            if (CSLCount(papszValues) >= 4)
            {
                const int nIndex = atoi(papszValues[0]);
                poRat->SetValue(nRatRow, 0, nIndex);
                poRat->SetValue(nRatRow, 1, atoi(papszValues[1]));
                poRat->SetValue(nRatRow, 2, atoi(papszValues[2]));
                poRat->SetValue(nRatRow, 3, atoi(papszValues[3]));
                nRatRow++;

                if (nIndex >= 0 && nIndex < 65536)
                {
                    const GDALColorEntry oEntry = {
                        static_cast<short>(
                            std::clamp(atoi(papszValues[1]), 0, 255)),  // Red
                        static_cast<short>(
                            std::clamp(atoi(papszValues[2]), 0, 255)),  // Green
                        static_cast<short>(
                            std::clamp(atoi(papszValues[3]), 0, 255)),  // Blue
                        255};

                    poDS->m_poColorTable->SetColorEntry(nIndex, &oEntry);
                }
                else
                {
                    // Negative values are valid. At least we can find use of
                    // them here:
                    //   http://www.ngdc.noaa.gov/mgg/topo/elev/esri/clr/
                    // But, there's no way of representing them with GDAL color
                    // table model.
                    if (!bHasFoundNonCTValues)
                        CPLDebug("EHdr", "Ignoring color index : %d", nIndex);
                    bHasFoundNonCTValues = true;
                }
            }

            CSLDestroy(papszValues);
        }

        CPL_IGNORE_RET_VAL(VSIFCloseL(fp));

        if (bHasFoundNonCTValues)
        {
            poDS->m_poRAT.swap(poRat);
        }

        for (int i = 1; i <= poDS->nBands; i++)
        {
            EHdrRasterBand *poBand =
                cpl::down_cast<EHdrRasterBand *>(poDS->GetRasterBand(i));
            poBand->m_poColorTable = poDS->m_poColorTable;
            poBand->m_poRAT = poDS->m_poRAT;
            poBand->SetColorInterpretation(GCI_PaletteIndex);
        }

        poDS->bCLRDirty = false;
    }

    // Read statistics (.STX).
    poDS->ReadSTX();

    // Initialize any PAM information.
    poDS->TryLoadXML();

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

    return poDS.release();
}

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

GDALDataset *EHdrDataset::Create(const char *pszFilename, int nXSize,
                                 int nYSize, int nBandsIn, GDALDataType eType,
                                 char **papszParamList)

{
    // Verify input options.
    if (nBandsIn <= 0)
    {
        CPLError(CE_Failure, CPLE_NotSupported,
                 "EHdr driver does not support %d bands.", nBandsIn);
        return nullptr;
    }

    if (eType != GDT_Byte && eType != GDT_Int8 && eType != GDT_Float32 &&
        eType != GDT_UInt16 && eType != GDT_Int16 && eType != GDT_Int32 &&
        eType != GDT_UInt32)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Attempt to create ESRI .hdr labelled dataset with an illegal"
                 "data type (%s).",
                 GDALGetDataTypeName(eType));

        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.", pszFilename);
        return nullptr;
    }

    // Just write out a couple of bytes to establish the binary
    // file, and then close it.
    bool bOK = VSIFWriteL(reinterpret_cast<void *>(const_cast<char *>("\0\0")),
                          2, 1, fp) == 1;
    if (VSIFCloseL(fp) != 0)
        bOK = false;
    fp = nullptr;
    if (!bOK)
        return nullptr;

    // Create the hdr filename.
    const std::string osHdrFilename = CPLResetExtensionSafe(pszFilename, "hdr");

    // Open the file.
    fp = VSIFOpenL(osHdrFilename.c_str(), "wt");
    if (fp == nullptr)
    {
        CPLError(CE_Failure, CPLE_OpenFailed,
                 "Attempt to create file `%s' failed.", osHdrFilename.c_str());
        return nullptr;
    }

    // Decide how many bits the file should have.
    const char *pszNBITS = CSLFetchNameValue(papszParamList, "NBITS");
    const int nBits =
        pszNBITS ? atoi(pszNBITS) : GDALGetDataTypeSizeBits(eType);
    if (nBits <= 0 || nXSize > (static_cast<int64_t>(INT_MAX) * 8 - 7) / nBits)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Invalid NBITS or too large image width");
        CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
        return nullptr;
    }

    const int nRowBytes = (nBits * nXSize + 7) / 8;

    if (nRowBytes > INT_MAX / nBandsIn)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Too large band count and/or image width");
        CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
        return nullptr;
    }

    // Check for signed byte.
    const char *pszPixelType = CSLFetchNameValue(papszParamList, "PIXELTYPE");
    if (pszPixelType == nullptr)
        pszPixelType = "";

    // Write out the raw definition for the dataset as a whole.
    bOK &= VSIFPrintfL(fp, "BYTEORDER      I\n") >= 0;
    bOK &= VSIFPrintfL(fp, "LAYOUT         BIL\n") >= 0;
    bOK &= VSIFPrintfL(fp, "NROWS          %d\n", nYSize) >= 0;
    bOK &= VSIFPrintfL(fp, "NCOLS          %d\n", nXSize) >= 0;
    bOK &= VSIFPrintfL(fp, "NBANDS         %d\n", nBandsIn) >= 0;
    bOK &= VSIFPrintfL(fp, "NBITS          %d\n", nBits) >= 0;
    bOK &= VSIFPrintfL(fp, "BANDROWBYTES   %d\n", nRowBytes) >= 0;
    bOK &= VSIFPrintfL(fp, "TOTALROWBYTES  %d\n", nRowBytes * nBandsIn) >= 0;

    if (eType == GDT_Float32)
        bOK &= VSIFPrintfL(fp, "PIXELTYPE      FLOAT\n") >= 0;
    else if (eType == GDT_Int8 || eType == GDT_Int16 || eType == GDT_Int32)
        bOK &= VSIFPrintfL(fp, "PIXELTYPE      SIGNEDINT\n") >= 0;
    else if (eType == GDT_Byte && EQUAL(pszPixelType, "SIGNEDBYTE"))
        bOK &= VSIFPrintfL(fp, "PIXELTYPE      SIGNEDINT\n") >= 0;
    else
        bOK &= VSIFPrintfL(fp, "PIXELTYPE      UNSIGNEDINT\n") >= 0;

    if (VSIFCloseL(fp) != 0)
        bOK = false;

    if (!bOK)
        return nullptr;

    GDALOpenInfo oOpenInfo(pszFilename, GA_Update);
    return Open(&oOpenInfo, false);
}

/************************************************************************/
/*                             CreateCopy()                             */
/************************************************************************/

GDALDataset *EHdrDataset::CreateCopy(const char *pszFilename,
                                     GDALDataset *poSrcDS, int bStrict,
                                     char **papszOptions,
                                     GDALProgressFunc pfnProgress,
                                     void *pProgressData)

{
    const int nBands = poSrcDS->GetRasterCount();
    if (nBands == 0)
    {
        CPLError(CE_Failure, CPLE_NotSupported,
                 "EHdr driver does not support source dataset without any "
                 "bands.");
        return nullptr;
    }

    char **papszAdjustedOptions = CSLDuplicate(papszOptions);

    // Ensure we pass on NBITS and PIXELTYPE structure information.
    auto poSrcBand = poSrcDS->GetRasterBand(1);
    if (poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE") != nullptr &&
        CSLFetchNameValue(papszOptions, "NBITS") == nullptr)
    {
        papszAdjustedOptions = CSLSetNameValue(
            papszAdjustedOptions, "NBITS",
            poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE"));
    }

    if (poSrcBand->GetRasterDataType() == GDT_Byte &&
        CSLFetchNameValue(papszOptions, "PIXELTYPE") == nullptr)
    {
        poSrcBand->EnablePixelTypeSignedByteWarning(false);
        const char *pszPixelType =
            poSrcBand->GetMetadataItem("PIXELTYPE", "IMAGE_STRUCTURE");
        poSrcBand->EnablePixelTypeSignedByteWarning(true);
        if (pszPixelType != nullptr)
        {
            papszAdjustedOptions = CSLSetNameValue(papszAdjustedOptions,
                                                   "PIXELTYPE", pszPixelType);
        }
    }

    // Proceed with normal copying using the default createcopy  operators.
    GDALDriver *poDriver =
        reinterpret_cast<GDALDriver *>(GDALGetDriverByName("EHdr"));

    GDALDataset *poOutDS = poDriver->DefaultCreateCopy(
        pszFilename, poSrcDS, bStrict, papszAdjustedOptions, pfnProgress,
        pProgressData);
    CSLDestroy(papszAdjustedOptions);

    if (poOutDS != nullptr)
        poOutDS->FlushCache(false);

    return poOutDS;
}

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

double EHdrRasterBand::GetNoDataValue(int *pbSuccess)
{
    if (pbSuccess)
        *pbSuccess = bNoDataSet;

    if (bNoDataSet)
        return dfNoData;

    return RawRasterBand::GetNoDataValue(pbSuccess);
}

/************************************************************************/
/*                           GetMinimum()                               */
/************************************************************************/

double EHdrRasterBand::GetMinimum(int *pbSuccess)
{
    if (pbSuccess != nullptr)
        *pbSuccess = (minmaxmeanstddev & HAS_MIN_FLAG) != 0;

    if (minmaxmeanstddev & HAS_MIN_FLAG)
        return dfMin;

    return RawRasterBand::GetMinimum(pbSuccess);
}

/************************************************************************/
/*                           GetMaximum()                               */
/************************************************************************/

double EHdrRasterBand::GetMaximum(int *pbSuccess)
{
    if (pbSuccess != nullptr)
        *pbSuccess = (minmaxmeanstddev & HAS_MAX_FLAG) != 0;

    if (minmaxmeanstddev & HAS_MAX_FLAG)
        return dfMax;

    return RawRasterBand::GetMaximum(pbSuccess);
}

/************************************************************************/
/*                           GetStatistics()                            */
/************************************************************************/

CPLErr EHdrRasterBand::GetStatistics(int bApproxOK, int bForce, double *pdfMin,
                                     double *pdfMax, double *pdfMean,
                                     double *pdfStdDev)
{
    if (!(GetMetadataItem("STATISTICS_APPROXIMATE") && !bApproxOK))
    {
        if ((minmaxmeanstddev & HAS_ALL_FLAGS) == HAS_ALL_FLAGS)
        {
            if (pdfMin)
                *pdfMin = dfMin;
            if (pdfMax)
                *pdfMax = dfMax;
            if (pdfMean)
                *pdfMean = dfMean;
            if (pdfStdDev)
                *pdfStdDev = dfStdDev;
            return CE_None;
        }
    }

    const CPLErr eErr = RawRasterBand::GetStatistics(
        bApproxOK, bForce, &dfMin, &dfMax, &dfMean, &dfStdDev);
    if (eErr != CE_None)
        return eErr;

    EHdrDataset *poEDS = cpl::down_cast<EHdrDataset *>(poDS);

    minmaxmeanstddev = HAS_ALL_FLAGS;

    if (!bApproxOK && poEDS->RewriteSTX() != CE_None)
        RawRasterBand::SetStatistics(dfMin, dfMax, dfMean, dfStdDev);

    if (pdfMin)
        *pdfMin = dfMin;
    if (pdfMax)
        *pdfMax = dfMax;
    if (pdfMean)
        *pdfMean = dfMean;
    if (pdfStdDev)
        *pdfStdDev = dfStdDev;

    return CE_None;
}

/************************************************************************/
/*                           SetStatistics()                            */
/************************************************************************/

CPLErr EHdrRasterBand::SetStatistics(double dfMinIn, double dfMaxIn,
                                     double dfMeanIn, double dfStdDevIn)
{
    // Avoid churn if nothing is changing.
    if (dfMin == dfMinIn && dfMax == dfMaxIn && dfMean == dfMeanIn &&
        dfStdDev == dfStdDevIn)
        return CE_None;

    dfMin = dfMinIn;
    dfMax = dfMaxIn;
    dfMean = dfMeanIn;
    dfStdDev = dfStdDevIn;

    // marks stats valid
    minmaxmeanstddev = HAS_ALL_FLAGS;

    EHdrDataset *poEDS = cpl::down_cast<EHdrDataset *>(poDS);

    if (GetMetadataItem("STATISTICS_APPROXIMATE") == nullptr)
    {
        if (GetMetadataItem("STATISTICS_MINIMUM"))
        {
            SetMetadataItem("STATISTICS_MINIMUM", nullptr);
            SetMetadataItem("STATISTICS_MAXIMUM", nullptr);
            SetMetadataItem("STATISTICS_MEAN", nullptr);
            SetMetadataItem("STATISTICS_STDDEV", nullptr);
        }
        return poEDS->RewriteSTX();
    }

    return RawRasterBand::SetStatistics(dfMinIn, dfMaxIn, dfMeanIn, dfStdDevIn);
}

/************************************************************************/
/*                           GetColorTable()                            */
/************************************************************************/

GDALColorTable *EHdrRasterBand::GetColorTable()
{
    return m_poColorTable.get();
}

/************************************************************************/
/*                           SetColorTable()                            */
/************************************************************************/

CPLErr EHdrRasterBand::SetColorTable(GDALColorTable *poNewCT)
{
    if (poNewCT == nullptr)
        m_poColorTable.reset();
    else
        m_poColorTable.reset(poNewCT->Clone());

    cpl::down_cast<EHdrDataset *>(poDS)->bCLRDirty = true;

    return CE_None;
}

/************************************************************************/
/*                           GetDefaultRAT()                            */
/************************************************************************/

GDALRasterAttributeTable *EHdrRasterBand::GetDefaultRAT()
{
    return m_poRAT.get();
}

/************************************************************************/
/*                            SetDefaultRAT()                           */
/************************************************************************/

CPLErr EHdrRasterBand::SetDefaultRAT(const GDALRasterAttributeTable *poRAT)
{
    if (poRAT)
    {
        if (!(poRAT->GetColumnCount() == 4 &&
              poRAT->GetTypeOfCol(0) == GFT_Integer &&
              poRAT->GetTypeOfCol(1) == GFT_Integer &&
              poRAT->GetTypeOfCol(2) == GFT_Integer &&
              poRAT->GetTypeOfCol(3) == GFT_Integer &&
              poRAT->GetUsageOfCol(0) == GFU_Generic &&
              poRAT->GetUsageOfCol(1) == GFU_Red &&
              poRAT->GetUsageOfCol(2) == GFU_Green &&
              poRAT->GetUsageOfCol(3) == GFU_Blue))
        {
            CPLError(CE_Warning, CPLE_NotSupported,
                     "Unsupported type of RAT: "
                     "only value,R,G,B ones are supported");
            return CE_Failure;
        }
    }

    if (poRAT == nullptr)
        m_poRAT.reset();
    else
        m_poRAT.reset(poRAT->Clone());

    cpl::down_cast<EHdrDataset *>(poDS)->bCLRDirty = true;

    return CE_None;
}

/************************************************************************/
/*                         GDALRegister_EHdr()                          */
/************************************************************************/

void GDALRegister_EHdr()

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

    GDALDriver *poDriver = new GDALDriver();

    poDriver->SetDescription("EHdr");
    poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "ESRI .hdr Labelled");
    poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/ehdr.html");
    poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "bil");
    poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
                              "Byte Int8 Int16 UInt16 Int32 UInt32 Float32");

    poDriver->SetMetadataItem(
        GDAL_DMD_CREATIONOPTIONLIST,
        "<CreationOptionList>"
        "   <Option name='NBITS' type='int' description='Special pixel bits "
        "(1-7)'/>"
        "   <Option name='PIXELTYPE' type='string' description='By setting "
        "this to SIGNEDBYTE, a new Byte file can be forced to be written as "
        "signed byte'/>"
        "</CreationOptionList>");

    poDriver->SetMetadataItem(GDAL_DCAP_UPDATE, "YES");
    poDriver->SetMetadataItem(GDAL_DMD_UPDATE_ITEMS, "GeoTransform SRS NoData "
                                                     "RasterValues");

    poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    poDriver->pfnOpen = EHdrDataset::Open;
    poDriver->pfnCreate = EHdrDataset::Create;
    poDriver->pfnCreateCopy = EHdrDataset::CreateCopy;

    GetGDALDriverManager()->RegisterDriver(poDriver);
}