/******************************************************************************
 *
 * Project:  E00 grid driver
 * Purpose:  GDALDataset driver for E00 grid dataset.
 * Author:   Even Rouault, <even dot rouault at mines dash paris dot org>
 *
 ******************************************************************************
 * Copyright (c) 2011, Even Rouault <even dot rouault at mines-paris dot org>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ****************************************************************************/

#include "cpl_string.h"
#include "cpl_vsi_virtual.h"
#include "gdal_frmts.h"
#include "gdal_pam.h"
#include "ogr_spatialref.h"

/* Private import of e00read.c */
#define E00ReadOpen         GDALE00GRIDReadOpen
#define E00ReadCallbackOpen GDALE00GRIDReadCallbackOpen
#define E00ReadClose        GDALE00GRIDReadClose
#define E00ReadNextLine     GDALE00GRIDReadNextLine
#define E00ReadRewind       GDALE00GRIDReadRewind

CPL_CVSID("$Id: e00griddataset.cpp dc782e589a0abebc628dcb9a931d2424d1ec7d17 2018-09-03 19:46:58 +0200 Even Rouault $")

#undef NULL
#define NULL nullptr
#include "e00read.c"
#undef NULL

constexpr int E00_INT_SIZE = 10;
// constexpr int E00_INT14_SIZE = 14;
constexpr int E00_FLOAT_SIZE = 14;
constexpr int E00_DOUBLE_SIZE = 21;
constexpr int VALS_PER_LINE = 5;

/* g++ -fPIC -Wall -g frmts/e00grid/e00griddataset.cpp -shared -o gdal_E00GRID.so -Iport -Igcore -Iogr -L. -lgdal */

/* Test data ; (google for "EXP  0" "GRD  2")

ftp://msdis.missouri.edu/pub/dem/24k/county/
http://dusk.geo.orst.edu/djl/samoa/data/samoa_bathy.e00
http://dusk.geo.orst.edu/djl/samoa/FBNMS/RasterGrids-Metadata/ntae02_3m_utm.e00
http://www.navdat.org/coverages/elevation/iddem1.e00        (int32)
http://delta-vision.projects.atlas.ca.gov/lidar/bare_earth.grids/sac0165.e00
http://ag.arizona.edu/SRER/maps_e00/srer_dem.e00
http://ok.water.usgs.gov/projects/norlan/spatial/ntopo0408-10.e00 (compressed)
http://wrri.nmsu.edu/publish/techrpt/tr322/GIS/dem.e00 (compressed)
*/

/************************************************************************/
/* ==================================================================== */
/*                            E00GRIDDataset                            */
/* ==================================================================== */
/************************************************************************/

class E00GRIDRasterBand;

class E00GRIDDataset : public GDALPamDataset
{
    friend class E00GRIDRasterBand;

    E00ReadPtr  e00ReadPtr;
    VSILFILE   *fp;
    vsi_l_offset nDataStart;
    int         nBytesEOL;

    vsi_l_offset  nPosBeforeReadLine;
    vsi_l_offset* panOffsets;
    int         nLastYOff;
    int         nMaxYOffset;

    double      adfGeoTransform[6];
    CPLString   osProjection;

    double      dfNoData;

    char**      papszPrj;

    const char* ReadLine();

    int         bHasReadMetadata;
    void        ReadMetadata();

    int         bHasStats;
    double      dfMin, dfMax, dfMean, dfStddev;

    static const char* ReadNextLine(void * ptr);
    static void        Rewind(void* ptr);

  public:
    E00GRIDDataset();
    ~E00GRIDDataset() override;

    CPLErr GetGeoTransform( double * ) override;
    const char* GetProjectionRef() override;

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

/************************************************************************/
/* ==================================================================== */
/*                          E00GRIDRasterBand                           */
/* ==================================================================== */
/************************************************************************/

class E00GRIDRasterBand : public GDALPamRasterBand
{
    friend class E00GRIDDataset;

  public:
    E00GRIDRasterBand( E00GRIDDataset *, int, GDALDataType );

    CPLErr IReadBlock( int, int, void * ) override;

    double GetNoDataValue( int *pbSuccess = nullptr ) override;
    const char *GetUnitType() override;
    double GetMinimum( int *pbSuccess = nullptr ) override;
    double GetMaximum( int *pbSuccess = nullptr ) override;
    CPLErr GetStatistics( int bApproxOK, int bForce,
                          double *pdfMin, double *pdfMax,
                          double *pdfMean, double *padfStdDev ) override;
};

/************************************************************************/
/*                         E00GRIDRasterBand()                          */
/************************************************************************/

E00GRIDRasterBand::E00GRIDRasterBand( E00GRIDDataset *poDSIn, int nBandIn,
                                      GDALDataType eDT )

{
    poDS = poDSIn;
    nBand = nBandIn;

    eDataType = eDT;

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

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

CPLErr E00GRIDRasterBand::IReadBlock( CPL_UNUSED int nBlockXOff,
                                      int nBlockYOff,
                                      void * pImage )
{
    E00GRIDDataset *poGDS = (E00GRIDDataset *) poDS;

    char szVal[E00_FLOAT_SIZE+1];
    szVal[E00_FLOAT_SIZE] = 0;

    float* pafImage = (float*)pImage;
    int* panImage = (int*)pImage;
    const float fNoData = static_cast<float>(poGDS->dfNoData);

    /* A new data line begins on a new text line. So if the xsize */
    /* is not a multiple of VALS_PER_LINE, there are padding values */
    /* that must be ignored */
    const int nRoundedBlockXSize = ((nBlockXSize + VALS_PER_LINE - 1) /
                                            VALS_PER_LINE) * VALS_PER_LINE;

    if (poGDS->e00ReadPtr)
    {
        if (poGDS->nLastYOff < 0)
        {
            E00ReadRewind(poGDS->e00ReadPtr);
            for(int i=0;i<6;i++)
                E00ReadNextLine(poGDS->e00ReadPtr);
        }

        if (nBlockYOff == poGDS->nLastYOff + 1)
        {
        }
        else if (nBlockYOff <= poGDS->nMaxYOffset)
        {
            //CPLDebug("E00GRID", "Skip to %d from %d", nBlockYOff, poGDS->nLastYOff);
            VSIFSeekL(poGDS->fp, poGDS->panOffsets[nBlockYOff], SEEK_SET);
            poGDS->nPosBeforeReadLine = poGDS->panOffsets[nBlockYOff];
            poGDS->e00ReadPtr->iInBufPtr = 0;
            poGDS->e00ReadPtr->szInBuf[0] = '\0';
        }
        else if (nBlockYOff > poGDS->nLastYOff + 1)
        {
            //CPLDebug("E00GRID", "Forward skip to %d from %d", nBlockYOff, poGDS->nLastYOff);
            for(int i=poGDS->nLastYOff + 1; i < nBlockYOff;i++)
            {
                if( IReadBlock(0, i, pImage) != CE_None )
                    return CE_Failure;
            }
        }

        if (nBlockYOff > poGDS->nMaxYOffset)
        {
            poGDS->panOffsets[nBlockYOff] = poGDS->nPosBeforeReadLine +
                                            poGDS->e00ReadPtr->iInBufPtr;
            poGDS->nMaxYOffset = nBlockYOff;
        }

        const char* pszLine = nullptr;
        for(int i=0;i<nBlockXSize;i++)
        {
            if ((i % VALS_PER_LINE) == 0)
            {
                pszLine = E00ReadNextLine(poGDS->e00ReadPtr);
                if (pszLine == nullptr || strlen(pszLine) < 5 * E00_FLOAT_SIZE)
                {
                    CPLError(CE_Failure, CPLE_FileIO,
                             "Could not find enough values for line %d",
                             nBlockYOff);
                    return CE_Failure;
                }
            }
            if (eDataType == GDT_Float32)
            {
                pafImage[i] = (float) CPLAtof(pszLine + (i%VALS_PER_LINE) * E00_FLOAT_SIZE);
                /* Workaround single vs double precision problems */
                if (fNoData != 0 && fabs((pafImage[i] - fNoData)/fNoData) < 1e-6)
                    pafImage[i] = fNoData;
            }
            else
            {
                panImage[i] = atoi(pszLine + (i%VALS_PER_LINE) * E00_FLOAT_SIZE);
            }
        }

        poGDS->nLastYOff = nBlockYOff;

        return CE_None;
    }

    vsi_l_offset nValsToSkip = (vsi_l_offset)nBlockYOff * nRoundedBlockXSize;
    vsi_l_offset nLinesToSkip = nValsToSkip / VALS_PER_LINE;
    int nBytesPerLine = VALS_PER_LINE * E00_FLOAT_SIZE + poGDS->nBytesEOL;
    vsi_l_offset nPos = poGDS->nDataStart + nLinesToSkip * nBytesPerLine;
    VSIFSeekL(poGDS->fp, nPos, SEEK_SET);

    for(int i=0;i<nBlockXSize;i++)
    {
        if (VSIFReadL(szVal, E00_FLOAT_SIZE, 1, poGDS->fp) != 1)
        {
            CPLError(CE_Failure, CPLE_FileIO,
                             "Could not find enough values for line %d",
                             nBlockYOff);
            return CE_Failure;
        }

        if (eDataType == GDT_Float32)
        {
            pafImage[i] = (float) CPLAtof(szVal);
            /* Workaround single vs double precision problems */
            if (fNoData != 0 && fabs((pafImage[i] - fNoData)/fNoData) < 1e-6)
                pafImage[i] = fNoData;
        }
        else
        {
            panImage[i] = atoi(szVal);
        }

        if (((i+1) % VALS_PER_LINE) == 0)
            VSIFReadL(szVal, poGDS->nBytesEOL, 1, poGDS->fp);
    }

    return CE_None;
}

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

double E00GRIDRasterBand::GetNoDataValue( int *pbSuccess )
{
    E00GRIDDataset *poGDS = (E00GRIDDataset *) poDS;

    if (pbSuccess)
        *pbSuccess = TRUE;

    if (eDataType == GDT_Float32)
        return (double)(float) poGDS->dfNoData;

    return (double)(int)poGDS->dfNoData;
}

/************************************************************************/
/*                             GetUnitType()                            */
/************************************************************************/

const char * E00GRIDRasterBand::GetUnitType()
{
    E00GRIDDataset *poGDS = (E00GRIDDataset *) poDS;

    poGDS->ReadMetadata();

    if (poGDS->papszPrj == nullptr)
        return GDALPamRasterBand::GetUnitType();

    char** papszIter = poGDS->papszPrj;
    const char* pszRet = "";
    while(*papszIter)
    {
        if (STARTS_WITH_CI(*papszIter, "Zunits"))
        {
            char** papszTokens = CSLTokenizeString(*papszIter);
            if (CSLCount(papszTokens) == 2)
            {
                if (EQUAL(papszTokens[1], "FEET"))
                    pszRet = "ft";
                else if (EQUAL(papszTokens[1], "METERS"))
                    pszRet = "m";
            }
            CSLDestroy(papszTokens);
            break;
        }
        papszIter ++;
    }

    return pszRet;
}

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

double E00GRIDRasterBand::GetMinimum( int *pbSuccess )
{
    E00GRIDDataset *poGDS = (E00GRIDDataset *) poDS;

    poGDS->ReadMetadata();

    if (poGDS->bHasStats)
    {
        if( pbSuccess != nullptr )
            *pbSuccess = TRUE;

        return poGDS->dfMin;
    }

    return GDALPamRasterBand::GetMinimum( pbSuccess );
}

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

double E00GRIDRasterBand::GetMaximum( int *pbSuccess )
{
    E00GRIDDataset *poGDS = (E00GRIDDataset *) poDS;

    poGDS->ReadMetadata();

    if (poGDS->bHasStats)
    {
        if( pbSuccess != nullptr )
            *pbSuccess = TRUE;

        return poGDS->dfMax;
    }

    return GDALPamRasterBand::GetMaximum( pbSuccess );
}

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

CPLErr E00GRIDRasterBand::GetStatistics( int bApproxOK, int bForce,
                                         double *pdfMin, double *pdfMax,
                                         double *pdfMean, double *pdfStdDev )
{
    E00GRIDDataset *poGDS = (E00GRIDDataset *) poDS;

    poGDS->ReadMetadata();

    if (poGDS->bHasStats)
    {
        if (pdfMin)
            *pdfMin = poGDS->dfMin;
        if (pdfMax)
            *pdfMax = poGDS->dfMax;
        if (pdfMean)
            *pdfMean = poGDS->dfMean;
        if (pdfStdDev)
            *pdfStdDev = poGDS->dfStddev;
        return CE_None;
    }

    return GDALPamRasterBand::GetStatistics(bApproxOK, bForce,
                                            pdfMin, pdfMax,
                                            pdfMean, pdfStdDev);
}

/************************************************************************/
/*                           E00GRIDDataset()                           */
/************************************************************************/

E00GRIDDataset::E00GRIDDataset() :
    e00ReadPtr(nullptr),
    fp(nullptr),
    nDataStart(0),
    nBytesEOL(1),
    nPosBeforeReadLine(0),
    panOffsets(nullptr),
    nLastYOff(-1),
    nMaxYOffset(-1),
    dfNoData(0),
    papszPrj(nullptr),
    bHasReadMetadata(FALSE),
    bHasStats(FALSE),
    dfMin(0),
    dfMax(0),
    dfMean(0),
    dfStddev(0)
{
    adfGeoTransform[0] = 0;
    adfGeoTransform[1] = 1;
    adfGeoTransform[2] = 0;
    adfGeoTransform[3] = 0;
    adfGeoTransform[4] = 0;
    adfGeoTransform[5] = 1;
}

/************************************************************************/
/*                           ~E00GRIDDataset()                          */
/************************************************************************/

E00GRIDDataset::~E00GRIDDataset()

{
    FlushCache();
    if (fp)
        VSIFCloseL(fp);
    CSLDestroy(papszPrj);
    E00ReadClose(e00ReadPtr);
    CPLFree(panOffsets);
}

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

int E00GRIDDataset::Identify( GDALOpenInfo * poOpenInfo )
{
    if (poOpenInfo->nHeaderBytes == 0)
        return FALSE;

    if (!(STARTS_WITH_CI((const char*)poOpenInfo->pabyHeader, "EXP  0") ||
          STARTS_WITH_CI((const char*)poOpenInfo->pabyHeader, "EXP  1")))
        return FALSE;

    /* FIXME: handle GRD  3 if that ever exists ? */
    if (strstr((const char*)poOpenInfo->pabyHeader, "GRD  2") == nullptr)
        return FALSE;

    return TRUE;
}

/************************************************************************/
/*                            ReadNextLine()                            */
/************************************************************************/

const char* E00GRIDDataset::ReadNextLine(void * ptr)
{
    E00GRIDDataset* poDS = (E00GRIDDataset*) ptr;
    poDS->nPosBeforeReadLine = VSIFTellL(poDS->fp);
    return CPLReadLine2L(poDS->fp, 256, nullptr);
}

/************************************************************************/
/*                                Rewind()                              */
/************************************************************************/

void E00GRIDDataset::Rewind(void * ptr)
{
    E00GRIDDataset* poDS = (E00GRIDDataset*) ptr;
    VSIRewindL(poDS->fp);
}

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

GDALDataset *E00GRIDDataset::Open( GDALOpenInfo * poOpenInfo )

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

/* -------------------------------------------------------------------- */
/*      Find dataset characteristics                                    */
/* -------------------------------------------------------------------- */
    VSILFILE* fp = poOpenInfo->fpL;
    poOpenInfo->fpL = nullptr;

    if (poOpenInfo->eAccess == GA_Update)
    {
        CPLError( CE_Failure, CPLE_NotSupported,
                  "The E00GRID driver does not support update access to existing"
                  " datasets.\n" );
        VSIFCloseL(fp);
        return nullptr;
    }

/* -------------------------------------------------------------------- */
/*      Create a corresponding GDALDataset.                             */
/* -------------------------------------------------------------------- */
    E00GRIDDataset *poDS = new E00GRIDDataset();
    if (strstr((const char*)poOpenInfo->pabyHeader, "\r\n") != nullptr)
        poDS->nBytesEOL = 2;
    poDS->fp = fp;

    /* read EXP  0 or EXP  1 line */
    const char* pszLine = CPLReadLine2L(fp, 81, nullptr);
    if (pszLine == nullptr)
    {
        CPLDebug("E00GRID", "Bad 1st line");
        delete poDS;
        return nullptr;
    }
    bool bCompressed = STARTS_WITH_CI(pszLine, "EXP  1");

    E00ReadPtr e00ReadPtr = nullptr;
    if (bCompressed)
    {
        VSIRewindL(fp);
        e00ReadPtr = E00ReadCallbackOpen(poDS,
                                         E00GRIDDataset::ReadNextLine,
                                         E00GRIDDataset::Rewind);
        if (e00ReadPtr == nullptr)
        {
            delete poDS;
            return nullptr;
        }
        E00ReadNextLine(e00ReadPtr);
        poDS->e00ReadPtr = e00ReadPtr;
    }

    /* skip GRD  2 line */
    if (e00ReadPtr)
        pszLine = E00ReadNextLine(e00ReadPtr);
    else
        pszLine = CPLReadLine2L(fp, 81, nullptr);
    if (pszLine == nullptr || !STARTS_WITH_CI(pszLine, "GRD  2"))
    {
        CPLDebug("E00GRID", "Bad 2nd line");
        delete poDS;
        return nullptr;
    }

    /* read ncols, nrows and nodata value */
    if (e00ReadPtr)
        pszLine = E00ReadNextLine(e00ReadPtr);
    else
        pszLine = CPLReadLine2L(fp, 81, nullptr);
    if (pszLine == nullptr || strlen(pszLine) <
                E00_INT_SIZE+E00_INT_SIZE+2+E00_DOUBLE_SIZE)
    {
        CPLDebug("E00GRID", "Bad 3rd line");
        delete poDS;
        return nullptr;
    }

    const int nRasterXSize = atoi(pszLine);
    const int nRasterYSize = atoi(pszLine + E00_INT_SIZE);

    if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize) ||
        /* ASCII format. Huge rasters do not make sense */
        nRasterXSize > 100000 ||
        nRasterYSize > 100000 )
    {
        delete poDS;
        return nullptr;
    }

    GDALDataType eDT = GDT_Float32;

    if (STARTS_WITH_CI(pszLine + E00_INT_SIZE + E00_INT_SIZE, " 1"))
        eDT = GDT_Int32;
    else if (STARTS_WITH_CI(pszLine + E00_INT_SIZE + E00_INT_SIZE, " 2"))
        eDT = GDT_Float32;
    else
    {
        CPLDebug("E00GRID", "Unknown data type : %s", pszLine);
    }

    const double dfNoData = CPLAtof(pszLine + E00_INT_SIZE + E00_INT_SIZE + 2);

    /* read pixel size */
    if (e00ReadPtr)
        pszLine = E00ReadNextLine(e00ReadPtr);
    else
        pszLine = CPLReadLine2L(fp, 81, nullptr);
    if (pszLine == nullptr || strlen(pszLine) < 2*E00_DOUBLE_SIZE)
    {
        CPLDebug("E00GRID", "Bad 4th line");
        delete poDS;
        return nullptr;
    }
/*
    double dfPixelX = CPLAtof(pszLine);
    double dfPixelY = CPLAtof(pszLine + E00_DOUBLE_SIZE);
*/

    /* read xmin, ymin */
    if (e00ReadPtr)
        pszLine = E00ReadNextLine(e00ReadPtr);
    else
        pszLine = CPLReadLine2L(fp, 81, nullptr);
    if (pszLine == nullptr || strlen(pszLine) < 2*E00_DOUBLE_SIZE)
    {
        CPLDebug("E00GRID", "Bad 5th line");
        delete poDS;
        return nullptr;
    }
    const double dfMinX = CPLAtof(pszLine);
    const double dfMinY = CPLAtof(pszLine + E00_DOUBLE_SIZE);

    /* read xmax, ymax */
    if (e00ReadPtr)
        pszLine = E00ReadNextLine(e00ReadPtr);
    else
        pszLine = CPLReadLine2L(fp, 81, nullptr);
    if (pszLine == nullptr || strlen(pszLine) < 2*E00_DOUBLE_SIZE)
    {
        CPLDebug("E00GRID", "Bad 6th line");
        delete poDS;
        return nullptr;
    }
    const double dfMaxX = CPLAtof(pszLine);
    const double dfMaxY = CPLAtof(pszLine + E00_DOUBLE_SIZE);

    poDS->nRasterXSize = nRasterXSize;
    poDS->nRasterYSize = nRasterYSize;
    poDS->dfNoData = dfNoData;
    poDS->adfGeoTransform[0] = dfMinX;
    poDS->adfGeoTransform[1] = (dfMaxX - dfMinX) / nRasterXSize;
    poDS->adfGeoTransform[2] = 0;
    poDS->adfGeoTransform[3] = dfMaxY;
    poDS->adfGeoTransform[4] = 0;
    poDS->adfGeoTransform[5] = - (dfMaxY - dfMinY) / nRasterYSize;
    poDS->nDataStart = VSIFTellL(fp);
    if (bCompressed)
    {
        poDS->panOffsets = (vsi_l_offset*)
                        VSIMalloc2(sizeof(vsi_l_offset), nRasterYSize);
        if (poDS->panOffsets == nullptr)
        {
            delete poDS;
            return nullptr;
        }
    }
/* -------------------------------------------------------------------- */
/*      Create band information objects.                                */
/* -------------------------------------------------------------------- */
    poDS->nBands = 1;
    for( int i = 0; i < poDS->nBands; i++ )
        poDS->SetBand( i+1, new E00GRIDRasterBand( poDS, i+1, eDT ) );

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

/* -------------------------------------------------------------------- */
/*      Support overviews.                                              */
/* -------------------------------------------------------------------- */
    poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
    return poDS;
}

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

CPLErr E00GRIDDataset::GetGeoTransform( double * padfTransform )

{
    memcpy(padfTransform, adfGeoTransform, 6 * sizeof(double));

    return CE_None;
}

/************************************************************************/
/*                             ReadLine()                               */
/************************************************************************/

const char* E00GRIDDataset::ReadLine()
{
    if (e00ReadPtr)
        return E00ReadNextLine(e00ReadPtr);

    return CPLReadLine2L(fp, 81, nullptr);
}

/************************************************************************/
/*                         GetProjectionRef()                           */
/************************************************************************/

const char* E00GRIDDataset::GetProjectionRef()

{
    ReadMetadata();
    return osProjection.c_str();
}

/************************************************************************/
/*                          ReadMetadata()                              */
/************************************************************************/

void E00GRIDDataset::ReadMetadata()

{
    if (bHasReadMetadata)
        return;

    bHasReadMetadata = TRUE;

    if (e00ReadPtr == nullptr)
    {
        const int nRoundedBlockXSize =
            DIV_ROUND_UP(nRasterXSize, VALS_PER_LINE) * VALS_PER_LINE;
        if( static_cast<vsi_l_offset>(nRasterYSize) >
                                    GUINTBIG_MAX / nRoundedBlockXSize )
        {
            return;
        }
        const vsi_l_offset nValsToSkip =
                               (vsi_l_offset)nRasterYSize * nRoundedBlockXSize;
        const vsi_l_offset nLinesToSkip = nValsToSkip / VALS_PER_LINE;
        const int nBytesPerLine = VALS_PER_LINE * E00_FLOAT_SIZE + nBytesEOL;
        const vsi_l_offset nPos = nDataStart + nLinesToSkip * nBytesPerLine;
        VSIFSeekL(fp, nPos, SEEK_SET);
    }
    else
    {
        nLastYOff = -1;

        const unsigned int BUFFER_SIZE = 65536;
        const unsigned int NEEDLE_SIZE = 3*5;
        const unsigned int nToRead = BUFFER_SIZE - NEEDLE_SIZE;
        char* pabyBuffer = (char*)CPLCalloc(1, BUFFER_SIZE+NEEDLE_SIZE);
        int nRead;
        int bEOGFound = FALSE;

        VSIFSeekL(fp, 0, SEEK_END);
        vsi_l_offset nEndPos = VSIFTellL(fp);
        if (nEndPos > BUFFER_SIZE)
            nEndPos -= BUFFER_SIZE;
        else
            nEndPos = 0;
        VSIFSeekL(fp, nEndPos, SEEK_SET);

#define GOTO_NEXT_CHAR() \
    i ++; \
    if (pabyBuffer[i] == 13 || pabyBuffer[i] == 10) \
    { \
        i++; \
        if (pabyBuffer[i] == 10) \
            i++; \
    } \

        while ((nRead = static_cast<int>(VSIFReadL(pabyBuffer, 1, nToRead, fp))) != 0)
        {
            int i;
            for(i = 0; i < nRead; i++)
            {
                if (pabyBuffer[i] == 'E')
                {
                    GOTO_NEXT_CHAR();
                    if (pabyBuffer[i] == 'O')
                    {
                        GOTO_NEXT_CHAR();
                        if (pabyBuffer[i] == 'G')
                        {
                            GOTO_NEXT_CHAR();
                            if (pabyBuffer[i] == '~')
                            {
                                GOTO_NEXT_CHAR();
                                if (pabyBuffer[i] == '}')
                                {
                                    bEOGFound = TRUE;
                                    break;
                                }
                            }
                        }
                    }
                }
            }

            if (bEOGFound)
            {
                VSIFSeekL(fp, VSIFTellL(fp) - nRead + i + 1, SEEK_SET);
                e00ReadPtr->iInBufPtr = 0;
                e00ReadPtr->szInBuf[0] = '\0';
                break;
            }

            if (nEndPos == 0)
                break;

            if ((unsigned int)nRead == nToRead)
            {
                memmove(pabyBuffer + nToRead, pabyBuffer, NEEDLE_SIZE);
                if (nEndPos >= (vsi_l_offset)nToRead)
                    nEndPos -= nToRead;
                else
                    nEndPos = 0;
                VSIFSeekL(fp, nEndPos, SEEK_SET);
            }
            else
                break;
        }
        CPLFree(pabyBuffer);
        if (!bEOGFound)
            return;
    }

    const char* pszLine = nullptr;
    bool bPRJFound = false;
    bool bStatsFound = false;
    while((pszLine = ReadLine()) != nullptr)
    {
        if (STARTS_WITH_CI(pszLine, "PRJ  2"))
        {
            bPRJFound = true;
            while((pszLine = ReadLine()) != nullptr)
            {
                if (EQUAL(pszLine, "EOP"))
                {
                    break;
                }
                if (!EQUAL(pszLine, "~") )
                {
                    papszPrj = CSLAddString(papszPrj, pszLine);
                }
            }

            OGRSpatialReference oSRS;
            if( oSRS.importFromESRI( papszPrj ) != OGRERR_NONE )
            {
                CPLError( CE_Warning, CPLE_AppDefined,
                            "Failed to parse PRJ section, ignoring." );
            }
            else
            {
                char* pszWKT = nullptr;
                if (oSRS.exportToWkt(&pszWKT) == OGRERR_NONE && pszWKT != nullptr)
                    osProjection = pszWKT;
                CPLFree(pszWKT);
            }
            if (bStatsFound)
                break;
        }
        else if (strcmp(pszLine, "STDV              8-1  254-1  15 3 60-1  -1  -1-1                   4-") == 0)
        {
            bStatsFound = true;
            pszLine = ReadLine();
            if (pszLine)
            {
                CPLString osStats = pszLine;
                pszLine = ReadLine();
                if (pszLine)
                {
                    osStats += pszLine;
                    char** papszTokens = CSLTokenizeString(osStats);
                    if (CSLCount(papszTokens) == 4)
                    {
                        dfMin = CPLAtof(papszTokens[0]);
                        dfMax = CPLAtof(papszTokens[1]);
                        dfMean = CPLAtof(papszTokens[2]);
                        dfStddev = CPLAtof(papszTokens[3]);
                        bHasStats = TRUE;
                    }
                    CSLDestroy(papszTokens);
                }
            }
            if (bPRJFound)
                break;
        }
    }
}

/************************************************************************/
/*                       GDALRegister_E00GRID()                         */
/************************************************************************/

void GDALRegister_E00GRID()

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

    GDALDriver *poDriver = new GDALDriver();

    poDriver->SetDescription( "E00GRID" );
    poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
    poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
                               "Arc/Info Export E00 GRID" );
    poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
                               "frmt_various.html#E00GRID" );
    poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "e00" );

    poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );

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

    GetGDALDriverManager()->RegisterDriver( poDriver );
}