Ticket #2191: gsagdataset.cpp

/******************************************************************************
 * $Id: gsagdataset.cpp 13500 2008-01-08 22:17:42Z rouault $
 *
 * Project:  GDAL
 * Purpose:  Implements the Golden Software ASCII Grid Format.
 * Author:   Kevin Locke, kwl7@cornell.edu
 *           (Based largely on aaigriddataset.cpp by Frank Warmerdam)
 *
 ******************************************************************************
 * Copyright (c) 2006, Kevin Locke <kwl7@cornell.edu>
 *
 * 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_conv.h"

#include <sstream>
#include <float.h>
#include <limits.h>
#include <assert.h>

#include "gdal_pam.h"

#ifndef DBL_MAX
# ifdef __DBL_MAX__
#  define DBL_MAX __DBL_MAX__
# else
#  define DBL_MAX 1.7976931348623157E+308
# endif /* __DBL_MAX__ */
#endif /* DBL_MAX */

#ifndef INT_MAX
# define INT_MAX 2147483647
#endif /* INT_MAX */

CPL_CVSID("$Id: gsagdataset.cpp 13500 2008-01-08 22:17:42Z rouault $");

CPL_C_START
void    GDALRegister_GSAG(void);
CPL_C_END

/************************************************************************/
/* ==================================================================== */
/*                              GSAGDataset                             */
/* ==================================================================== */
/************************************************************************/

class GSAGRasterBand;

class GSAGDataset : public GDALPamDataset
{
    friend class GSAGRasterBand;

    static const double dfNODATA_VALUE;
    static const int nFIELD_PRECISION;
    static const size_t nMAX_HEADER_SIZE;

    static CPLErr ShiftFileContents( FILE *, vsi_l_offset, int, const char * );

    FILE        *fp;
    size_t       nMinMaxZOffset;
    char         szEOL[3];

    CPLErr UpdateHeader();

  public:
                GSAGDataset( const char *pszEOL = "\x0D\x0A" );
                ~GSAGDataset();

    static GDALDataset *Open( GDALOpenInfo * );
    static GDALDataset *CreateCopy( const char *pszFilename,
                                    GDALDataset *poSrcDS,
                                    int bStrict, char **papszOptions,
                                    GDALProgressFunc pfnProgress,
                                    void *pProgressData );
    static CPLErr Delete( const char *pszFilename );

    CPLErr GetGeoTransform( double *padfGeoTransform );
    CPLErr SetGeoTransform( double *padfGeoTransform );
};

/* NOTE:  This is not mentioned in the spec, but Surfer 8 uses this value */
const double GSAGDataset::dfNODATA_VALUE = 1.70141E+38;

const int GSAGDataset::nFIELD_PRECISION = 14;
const size_t GSAGDataset::nMAX_HEADER_SIZE = 200;

/************************************************************************/
/* ==================================================================== */
/*                            GSAGRasterBand                            */
/* ==================================================================== */
/************************************************************************/

class GSAGRasterBand : public GDALPamRasterBand
{
    friend class GSAGDataset;

    double dfMinX;
    double dfMaxX;
    double dfMinY;
    double dfMaxY;
    double dfMinZ;
    double dfMaxZ;

    vsi_l_offset *panLineOffset;
        int nLastReadLine;

    double *padfRowMinZ;
    double *padfRowMaxZ;
    int nMinZRow;
    int nMaxZRow;

    CPLErr ScanForMinMaxZ();

  public:

                GSAGRasterBand( GSAGDataset *, int, vsi_l_offset );
                ~GSAGRasterBand();
    
    CPLErr IReadBlock( int, int, void * );
    CPLErr IWriteBlock( int, int, void * );

    double GetNoDataValue( int *pbSuccess = NULL );
    double GetMinimum( int *pbSuccess = NULL );
    double GetMaximum( int *pbSuccess = NULL );
};

/************************************************************************/
/*                            AlmostEqual()                             */
/* This function is needed because in release mode "1.70141E+38" is not */
/* parsed as 1.70141E+38 in the last bit of the mantissa.               */
/* See http://gcc.gnu.org/ml/gcc/2003-08/msg01195.html for some         */
/* explanation.                                                         */
/************************************************************************/

bool AlmostEqual( double dfVal1, double dfVal2 )

{
    const double dfTOLERANCE = 0.0000000001;
    if( dfVal1 == 0.0 || dfVal2 == 0.0 )
        return fabs(dfVal1 - dfVal2) < dfTOLERANCE;
    return fabs((dfVal1 - dfVal2)/dfVal1) < dfTOLERANCE;
}

/************************************************************************/
/*                           GSAGRasterBand()                           */
/************************************************************************/

GSAGRasterBand::GSAGRasterBand( GSAGDataset *poDS, int nBand,
                                vsi_l_offset nDataStart ) :
    padfRowMinZ(NULL),
    padfRowMaxZ(NULL),
    nMinZRow(-1),
    nMaxZRow(-1)

{
    this->poDS = poDS;
    nBand = nBand;
    
    eDataType = GDT_Float64;

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

    panLineOffset =
        (vsi_l_offset *)CPLCalloc( poDS->nRasterYSize+1, sizeof(vsi_l_offset) );
    if( panLineOffset == NULL )
        return;

        panLineOffset[poDS->nRasterYSize-1] = nDataStart;
        nLastReadLine = poDS->nRasterYSize;
}

/************************************************************************/
/*                          ~GSAGRasterBand()                           */
/************************************************************************/

GSAGRasterBand::~GSAGRasterBand()
{
    CPLFree( panLineOffset );
    if( padfRowMinZ != NULL )
        CPLFree( padfRowMinZ );
    if( padfRowMaxZ != NULL )
        CPLFree( padfRowMaxZ );
}

/************************************************************************/
/*                           ScanForMinMaxZ()                           */
/************************************************************************/

CPLErr GSAGRasterBand::ScanForMinMaxZ()

{
    double *padfRowValues = (double *)VSIMalloc2( nBlockXSize, sizeof(double) );
    if( padfRowValues == NULL )
    {
        CPLError( CE_Failure, CPLE_OutOfMemory,
                  "Unable to allocate memory for grid row values.\n" );
        return CE_Failure;
    }

    double dfNewMinZ = DBL_MAX;
    double dfNewMaxZ = -DBL_MAX;
    int nNewMinZRow = 0;
    int nNewMaxZRow = 0;

    /* Since we have to scan, lets calc. statistics too */
    double dfSum = 0.0;
    double dfSum2 = 0.0;
    unsigned long nValuesRead = 0;
    for( int iRow=0; iRow<nRasterYSize; iRow++ )
    {
        CPLErr eErr = IReadBlock( 0, iRow, padfRowValues );
        if( eErr != CE_None )
        {
            VSIFree( padfRowValues );
            return eErr;
        }

        padfRowMinZ[iRow] = DBL_MAX;
        padfRowMaxZ[iRow] = -DBL_MAX;
        for( int iCell=0; iCell<nRasterXSize; iCell++ )
        {
            if( AlmostEqual(padfRowValues[iCell], GSAGDataset::dfNODATA_VALUE) )
                continue;

            if( padfRowValues[iCell] < padfRowMinZ[iRow] )
                padfRowMinZ[iRow] = padfRowValues[iCell];

            if( padfRowValues[iCell] > padfRowMaxZ[iRow] )
                padfRowMaxZ[iRow] = padfRowValues[iCell];

            dfSum += padfRowValues[iCell];
            dfSum2 += padfRowValues[iCell] * padfRowValues[iCell];
            nValuesRead++;
        }

        if( padfRowMinZ[iRow] < dfNewMinZ )
        {
            dfNewMinZ = padfRowMinZ[iRow];
            nNewMinZRow = iRow;
        }

        if( padfRowMaxZ[iRow] > dfNewMaxZ )
        {
            dfNewMaxZ = padfRowMaxZ[iRow];
            nNewMaxZRow = iRow;
        }
    }

    VSIFree( padfRowValues );

    if( nValuesRead == 0 )
    {
        dfMinZ = 0.0;
        dfMaxZ = 0.0;
        nMinZRow = 0;
        nMaxZRow = 0;
        return CE_None;
    }

    dfMinZ = dfNewMinZ;
    dfMaxZ = dfNewMaxZ;
    nMinZRow = nNewMinZRow;
    nMaxZRow = nNewMaxZRow;

    double dfMean = dfSum / nValuesRead;
    double dfStdDev = sqrt((dfSum2 / nValuesRead) - (dfMean * dfMean));
    SetStatistics( dfMinZ, dfMaxZ, dfMean, dfStdDev );

    return CE_None;
}

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

CPLErr GSAGRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
                                   void * pImage )
{
    static size_t nMaxLineSize = 128;
    double *pdfImage = (double *)pImage;
    GSAGDataset *poGDS = (GSAGDataset *)poDS;

    assert( poGDS != NULL );

    if( nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0 )
        return CE_Failure;

    if( panLineOffset[nBlockYOff] == 0 )
    {
                // Discover the last read block
                for ( int iFoundLine = nLastReadLine - 1; iFoundLine > nBlockYOff; iFoundLine--)
                {
                        IReadBlock( nBlockXOff, iFoundLine, NULL);
                }
    }

    if( panLineOffset[nBlockYOff] == 0 )
        return CE_Failure;
    if( VSIFSeekL( poGDS->fp, panLineOffset[nBlockYOff], SEEK_SET ) != 0 )
    {
        CPLError( CE_Failure, CPLE_FileIO,
              "Can't seek to offset %ld to read grid row %d.",
              panLineOffset[nBlockYOff], nBlockYOff );
        return CE_Failure;
    }

    size_t nLineBufSize;
    char *szLineBuf = NULL;
    size_t nCharsRead;
    size_t nCharsExamined = 0;
    /* If we know the offsets, we can just read line directly */
    if( (nBlockYOff > 0) && ( panLineOffset[nBlockYOff-1] != 0 ) )
    {
        assert(panLineOffset[nBlockYOff-1] > panLineOffset[nBlockYOff]);
        nLineBufSize = panLineOffset[nBlockYOff-1]
                              - panLineOffset[nBlockYOff] + 1;
    }
    else
    {
        nLineBufSize = nMaxLineSize;
    }

    szLineBuf = (char *)VSIMalloc( nLineBufSize );
    if( szLineBuf == NULL )
    {
        CPLError( CE_Failure, CPLE_OutOfMemory,
                  "Unable to read block, unable to allocate line buffer.\n" );
        return CE_Failure;
    }

    nCharsRead = VSIFReadL( szLineBuf, 1, nLineBufSize - 1, poGDS->fp );
    if( nCharsRead == 0 )
    {
        VSIFree( szLineBuf );
        CPLError( CE_Failure, CPLE_FileIO,
                  "Can't read grid row %d at offset %ld.\n",
                  nBlockYOff, panLineOffset[nBlockYOff] );
        return CE_Failure;
    }
    szLineBuf[nCharsRead] = '\0';

    char *szStart = szLineBuf;
    char *szEnd = szStart;
    for( int iCell=0; iCell<nBlockXSize; szStart = szEnd )
    {
        double dfValue = CPLStrtod( szStart, &szEnd );
        if( szStart == szEnd )
        {
            /* No number found */

            /* Check if this was an expected failure */
            while( isspace( *szStart ) )
                szStart++;

            /* Found sign at end of input, seek back to re-read it */
            if ( (*szStart == '-' || *szStart == '+') && *(szStart+1) == '\0' )
            {
                if( VSIFSeekL( poGDS->fp, 
                               VSIFTellL( poGDS->fp)-1, 
                               SEEK_SET ) != 0 )
                {
                    VSIFree( szLineBuf );
                    CPLError( CE_Failure, CPLE_FileIO,
                              "Unable to seek in grid row %d "
                              "(offset %ld, seek %d).\n",
                              nBlockYOff, VSIFTellL(poGDS->fp),
                              -1 );

                    return CE_Failure;
                }
            }
            else if( *szStart != '\0' )
            {
                szEnd = szStart;
                while( !isspace( *szEnd ) && *szEnd != '\0' )
                    szEnd++;
                char cOldEnd = *szEnd;
                *szEnd = '\0';

                CPLError( CE_Warning, CPLE_FileIO,
                          "Unexpected value in grid row %d (expected floating "
                          "point value, found \"%s\").\n",
                          nBlockYOff, szStart );

                *szEnd = cOldEnd;

                szEnd = szStart;
                while( !isdigit( *szEnd ) && *szEnd != '.' && *szEnd != '\0' )
                    szEnd++;

                continue;
            }
            else if( static_cast<size_t>(szStart - szLineBuf) != nCharsRead )
            {
                CPLError( CE_Warning, CPLE_FileIO,
                          "Unexpected ASCII null-character in grid row %d at "
                          "offset %ld.\n", nBlockYOff, szStart - szLineBuf );

                while( *szStart == '\0' &&
                        static_cast<size_t>(szStart - szLineBuf) < nCharsRead )
                    szStart++;

                szEnd = szStart;
                continue;
            }

            nCharsExamined += szStart - szLineBuf;
            nCharsRead = VSIFReadL( szLineBuf, 1, nLineBufSize - 1, poGDS->fp );
            if( nCharsRead == 0 )
            {
                VSIFree( szLineBuf );
                CPLError( CE_Failure, CPLE_FileIO,
                          "Can't read portion of grid row %d at offset %ld.",
                          nBlockYOff, panLineOffset[nBlockYOff] );
                return CE_Failure;
            }
            szLineBuf[nCharsRead] = '\0';
            szStart = szEnd = szLineBuf;
            continue;
        }
        else if( *szEnd == '\0'
                || (*szEnd == '.' && *(szEnd+1) == '\0')
                || (*szEnd == '-' && *(szEnd+1) == '\0')
                || (*szEnd == '+' && *(szEnd+1) == '\0')
                || (*szEnd == 'E' && *(szEnd+1) == '\0')
                || (*szEnd == 'E' && *(szEnd+1) == '-' && *(szEnd+2) == '\0')
                || (*szEnd == 'E' && *(szEnd+1) == '+' && *(szEnd+2) == '\0')
                || (*szEnd == 'e' && *(szEnd+1) == '\0')
                || (*szEnd == 'e' && *(szEnd+1) == '-' && *(szEnd+2) == '\0')
                || (*szEnd == 'e' && *(szEnd+1) == '+' && *(szEnd+2) == '\0'))
        {
            /* Number was interrupted by a nul character */
            while( *szEnd != '\0' )
                szEnd++;

            if( static_cast<size_t>(szEnd - szLineBuf) != nCharsRead )
            {
                CPLError( CE_Warning, CPLE_FileIO,
                          "Unexpected ASCII null-character in grid row %d at "
                          "offset %ld.\n", nBlockYOff, szStart - szLineBuf );

                while( *szEnd == '\0' &&
                       static_cast<size_t>(szStart - szLineBuf) < nCharsRead )
                    szEnd++;

                continue;
            }

            /* End of buffer, could be interrupting a number */
            if( VSIFSeekL( poGDS->fp, szStart - szEnd, SEEK_CUR ) != 0 )
            {
                VSIFree( szLineBuf );
                CPLError( CE_Failure, CPLE_FileIO,
                          "Unable to seek in grid row %d (offset %ld, seek %d)"
                          ".\n", nBlockYOff, VSIFTellL(poGDS->fp),
                          szStart - szEnd );

                return CE_Failure;
            }
            nCharsExamined += szStart - szLineBuf;
            nCharsRead = VSIFReadL( szLineBuf, 1, nLineBufSize - 1, poGDS->fp );
            szLineBuf[nCharsRead] = '\0';

            if( nCharsRead == 0 )
            {
                VSIFree( szLineBuf );
                CPLError( CE_Failure, CPLE_FileIO,
                          "Can't read portion of grid row %d at offset %ld.",
                          nBlockYOff, panLineOffset[nBlockYOff] );
                return CE_Failure;
            }
            else if( nCharsRead > static_cast<size_t>(szEnd - szStart) )
            {
                /* Read new data, this was not really the end */
                szEnd = szStart = szLineBuf;
                continue;
            }

            /* This is really the last value and has no tailing newline */
            szStart = szLineBuf;
            szEnd = szLineBuf + nCharsRead;
        }

        if( pdfImage != NULL )
        {
            *(pdfImage+iCell) = dfValue;
        }

        iCell++;
    }

    while( *szEnd == ' ' )
        szEnd++;

    if( *szEnd != '\0' && *szEnd != poGDS->szEOL[0] )
        CPLDebug( "GSAG", "Grid row %d does not end with a newline.  "
                         "Possible skew.\n", nBlockYOff );

    while( isspace( *szEnd ) )
        szEnd++;

    nCharsExamined += szEnd - szLineBuf;

    if( nCharsExamined >= nMaxLineSize )
        nMaxLineSize = nCharsExamined + 1;

    panLineOffset[nBlockYOff - 1] = panLineOffset[nBlockYOff] + nCharsExamined;
        nLastReadLine = nBlockYOff;

    VSIFree( szLineBuf );

    return CE_None;
}

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

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

{
    if( eAccess == GA_ReadOnly )
    {
        CPLError( CE_Failure, CPLE_NoWriteAccess,
                  "Unable to write block, dataset opened read only.\n" );
        return CE_Failure;
    }

    if( nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0 )
        return CE_Failure;

    GSAGDataset *poGDS = (GSAGDataset *)poDS;
    assert( poGDS != NULL );

    if( padfRowMinZ == NULL || padfRowMaxZ == NULL
        || nMinZRow < 0 || nMaxZRow < 0 )
    {
        padfRowMinZ = (double *)VSIMalloc2( nRasterYSize,sizeof(double) );
        if( padfRowMinZ == NULL )
        {
            CPLError( CE_Failure, CPLE_OutOfMemory,
                      "Unable to allocate space for row minimums array.\n" );
            return CE_Failure;
        }

        padfRowMaxZ = (double *)VSIMalloc2( nRasterYSize,sizeof(double) );
        if( padfRowMaxZ == NULL )
        {
            VSIFree( padfRowMinZ );
            padfRowMinZ = NULL;
            CPLError( CE_Failure, CPLE_OutOfMemory,
                      "Unable to allocate space for row maximums array.\n" );
            return CE_Failure;
        }

        CPLErr eErr = ScanForMinMaxZ();
        if( eErr != CE_None )
            return eErr;
    }

    if( panLineOffset[nBlockYOff+1] == 0 )
        IReadBlock( nBlockXOff, nBlockYOff, NULL );

    if( panLineOffset[nBlockYOff+1] == 0 || panLineOffset[nBlockYOff] == 0 )
        return CE_Failure;

    std::ostringstream ssOutBuf;
    ssOutBuf.precision( GSAGDataset::nFIELD_PRECISION );
    ssOutBuf.setf( std::ios::uppercase );

    double *pdfImage = (double *)pImage;
    padfRowMinZ[nBlockYOff] = DBL_MAX;
    padfRowMaxZ[nBlockYOff] = -DBL_MAX;
    for( int iCell=0; iCell<nBlockXSize; )
    {
        for( int iCol=0; iCol<10 && iCell<nBlockXSize; iCol++, iCell++ )
        {
            if( AlmostEqual( pdfImage[iCell], GSAGDataset::dfNODATA_VALUE ) )
            {
                if( pdfImage[iCell] < padfRowMinZ[nBlockYOff] )
                    padfRowMinZ[nBlockYOff] = pdfImage[iCell];

                if( pdfImage[iCell] > padfRowMaxZ[nBlockYOff] )
                    padfRowMaxZ[nBlockYOff] = pdfImage[iCell];
            }

            ssOutBuf << pdfImage[iCell] << " ";
        }
        ssOutBuf << poGDS->szEOL;
    }
    ssOutBuf << poGDS->szEOL;

    CPLString sOut = ssOutBuf.str();
    if( sOut.length() != panLineOffset[nBlockYOff+1]-panLineOffset[nBlockYOff] )
    {
        int nShiftSize = sOut.length() - (panLineOffset[nBlockYOff+1]
                                          - panLineOffset[nBlockYOff]);
        if( nBlockYOff != poGDS->nRasterYSize
            && GSAGDataset::ShiftFileContents( poGDS->fp,
                                               panLineOffset[nBlockYOff+1],
                                               nShiftSize,
                                               poGDS->szEOL ) != CE_None )
        {
            CPLError( CE_Failure, CPLE_FileIO,
                      "Failure writing block, "
                      "unable to shift file contents.\n" );
            return CE_Failure;
        }

        for( size_t iLine=nBlockYOff+1;
             iLine < static_cast<unsigned>(poGDS->nRasterYSize+1)
                && panLineOffset[iLine] != 0; iLine++ )
            panLineOffset[iLine] += nShiftSize;
    }

    if( VSIFSeekL( poGDS->fp, panLineOffset[nBlockYOff], SEEK_SET ) != 0 )
    {
        CPLError( CE_Failure, CPLE_FileIO, "Unable to seek to grid line.\n" );
        return CE_Failure;
    }

    if( VSIFWriteL( sOut.c_str(), 1, sOut.length(),
                    poGDS->fp ) != sOut.length() )
    {
        CPLError( CE_Failure, CPLE_FileIO, "Unable to write grid block.\n" );
        return CE_Failure;
    }

    /* Update header as needed */
    bool bHeaderNeedsUpdate = false;
    if( nMinZRow == nBlockYOff && padfRowMinZ[nBlockYOff] > dfMinZ )
    {
        double dfNewMinZ = -DBL_MAX;
        for( int iRow=0; iRow<nRasterYSize; iRow++ )
        {
            if( padfRowMinZ[iRow] < dfNewMinZ )
            {
                dfNewMinZ = padfRowMinZ[iRow];
                nMinZRow = iRow;
            }
        }

        if( dfNewMinZ != dfMinZ )
        {
            dfMinZ = dfNewMinZ;
            bHeaderNeedsUpdate = true;
        }
    }

    if( nMaxZRow == nBlockYOff && padfRowMaxZ[nBlockYOff] < dfMaxZ )
    {
        double dfNewMaxZ = -DBL_MAX;
        for( int iRow=0; iRow<nRasterYSize; iRow++ )
        {
            if( padfRowMaxZ[iRow] > dfNewMaxZ )
            {
                dfNewMaxZ = padfRowMaxZ[iRow];
                nMaxZRow = iRow;
            }
        }

        if( dfNewMaxZ != dfMaxZ )
        {
            dfMaxZ = dfNewMaxZ;
            bHeaderNeedsUpdate = true;
        }
    }

    if( padfRowMinZ[nBlockYOff] < dfMinZ || padfRowMaxZ[nBlockYOff] > dfMaxZ )
    {
        if( padfRowMinZ[nBlockYOff] < dfMinZ )
        {
            dfMinZ = padfRowMinZ[nBlockYOff];
            nMinZRow = nBlockYOff;
        }

        if( padfRowMaxZ[nBlockYOff] > dfMaxZ )
        {
            dfMaxZ = padfRowMaxZ[nBlockYOff];
            nMaxZRow = nBlockYOff;
        }

        bHeaderNeedsUpdate = true;
    }

    if( bHeaderNeedsUpdate && dfMaxZ > dfMinZ )
    {
        CPLErr eErr = poGDS->UpdateHeader();
        return eErr;
    }

    return CE_None;
}

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

double GSAGRasterBand::GetNoDataValue( int * pbSuccess )
{
    if( pbSuccess )
        *pbSuccess = TRUE;

    return GSAGDataset::dfNODATA_VALUE;
}

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

double GSAGRasterBand::GetMinimum( int *pbSuccess )
{
    if( pbSuccess )
        *pbSuccess = TRUE;

    return dfMinZ;
}

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

double GSAGRasterBand::GetMaximum( int *pbSuccess )
{
    if( pbSuccess )
        *pbSuccess = TRUE;

    return dfMaxZ;
}

/************************************************************************/
/* ==================================================================== */
/*                              GSAGDataset                             */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                             GSAGDataset()                            */
/************************************************************************/

GSAGDataset::GSAGDataset( const char *pszEOL )

{
    if( pszEOL == NULL || EQUAL(pszEOL, "") )
    {
        CPLDebug( "GSAG", "GSAGDataset() created with invalid EOL string.\n" );
        this->szEOL[0] = '\x0D';
        this->szEOL[1] = '\x0A';
        this->szEOL[2] = '\0';
    }
    else
    {
        strncpy(this->szEOL, pszEOL, sizeof(this->szEOL));
        this->szEOL[sizeof(this->szEOL) - 1] = '\0';
    }
}

/************************************************************************/
/*                            ~GSAGDataset()                            */
/************************************************************************/

GSAGDataset::~GSAGDataset()

{
    FlushCache();
    if( fp != NULL )
        VSIFCloseL( fp );
}

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

GDALDataset *GSAGDataset::Open( GDALOpenInfo * poOpenInfo )

{
    /* Check for signature */
    if( poOpenInfo->nHeaderBytes < 5
        || !EQUALN((const char *) poOpenInfo->pabyHeader,"DSAA",4)
        || ( poOpenInfo->pabyHeader[4] != '\x0D'
             && poOpenInfo->pabyHeader[4] != '\x0A' ))
    {
        return NULL;
    }

    /* Identify the end of line marker (should be \x0D\x0A, but try some others)
     * (note that '\x0D' == '\r' and '\x0A' == '\n' on most systems) */
    char szEOL[3];
    szEOL[0] = poOpenInfo->pabyHeader[4];
    szEOL[1] = poOpenInfo->pabyHeader[5];
    szEOL[2] = '\0';
    if( szEOL[1] != '\x0D' && szEOL[1] != '\x0A' )
        szEOL[1] = '\0';

/* -------------------------------------------------------------------- */
/*      Create a corresponding GDALDataset.                             */
/* -------------------------------------------------------------------- */
    GSAGDataset *poDS = new GSAGDataset( szEOL );

/* -------------------------------------------------------------------- */
/*      Open file with large file API.                                  */
/* -------------------------------------------------------------------- */

    poDS->eAccess = poOpenInfo->eAccess;
    if( poOpenInfo->eAccess == GA_ReadOnly )
        poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
    else
        poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "r+b" );

    if( poDS->fp == NULL )
    {
        CPLError( CE_Failure, CPLE_OpenFailed, 
                  "VSIFOpenL(%s) failed unexpectedly.", 
                  poOpenInfo->pszFilename );
        delete poDS;
        return NULL;
    }
 
/* -------------------------------------------------------------------- */
/*      Read the header.                                                */
/* -------------------------------------------------------------------- */
    char *pabyHeader;
    bool bMustFreeHeader = false;
    if( poOpenInfo->nHeaderBytes >= static_cast<int>(nMAX_HEADER_SIZE) )
    {
        pabyHeader = (char *)poOpenInfo->pabyHeader;
    }
    else
    {
        bMustFreeHeader = true;
        pabyHeader = (char *)VSIMalloc( nMAX_HEADER_SIZE );
        if( pabyHeader == NULL )
        {
            CPLError( CE_Failure, CPLE_OutOfMemory,
                      "Unable to open dataset, unable to header buffer.\n" );
            return NULL;
        }

        size_t nRead = VSIFReadL( pabyHeader, 1, nMAX_HEADER_SIZE-1, poDS->fp );
        pabyHeader[nRead+1] = '\0';
    }

    const char *szErrorMsg;
    const char *szStart = pabyHeader + 5;
    char *szEnd;
    double dfTemp;
    double dfMinX;
    double dfMaxX;
    double dfMinY;
    double dfMaxY;
    double dfMinZ;
    double dfMaxZ;

    /* Parse number of X axis grid rows */
    long nTemp = strtol( szStart, &szEnd, 10 );
    if( szStart == szEnd || nTemp < 0l )
    {
        szErrorMsg = "Unable to parse the number of X axis grid columns.\n";
        goto error;
    }
    else if( nTemp > INT_MAX )
    {
        CPLError( CE_Warning, CPLE_AppDefined,
                  "Number of X axis grid columns not representable.\n" );
        poDS->nRasterXSize = INT_MAX;
    }
    else
    {
        poDS->nRasterXSize = static_cast<int>(nTemp);
    }
    szStart = szEnd;

    /* Parse number of Y axis grid rows */
    nTemp = strtol( szStart, &szEnd, 10 );
    if( szStart == szEnd || nTemp < 0l )
    {
        szErrorMsg = "Unable to parse the number of Y axis grid rows.\n";
        goto error;
    }
    else if( nTemp > INT_MAX )
    {
        CPLError( CE_Warning, CPLE_AppDefined,
                  "Number of Y axis grid rows not representable.\n" );
        poDS->nRasterYSize = INT_MAX;
    }
    else
    {
        poDS->nRasterYSize = static_cast<int>(nTemp);
    }
    szStart = szEnd;

    /* Parse the minimum X value of the grid */
    dfTemp = CPLStrtod( szStart, &szEnd );
    if( szStart == szEnd )
    {
        szErrorMsg = "Unable to parse the minimum X value.\n";
        goto error;
    }
    else
    {
        dfMinX = dfTemp;
    }
    szStart = szEnd;

    /* Parse the maximum X value of the grid */
    dfTemp = CPLStrtod( szStart, &szEnd );
    if( szStart == szEnd )
    {
        szErrorMsg = "Unable to parse the maximum X value.\n";
        goto error;
    }
    else