EVOLUTION-MANAGER

Edit File: aaigriddataset.cpp

/******************************************************************************
 * $Id: aaigriddataset.cpp 27729 2014-09-24 00:40:16Z goatbar $
 *
 * Project:  GDAL
 * Purpose:  Implements Arc/Info ASCII Grid Format.
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 2001, Frank Warmerdam (warmerdam@pobox.com)
 * Copyright (c) 2007-2012, Even Rouault <even dot rouault at mines-paris dot org>
 * Copyright (c) 2014, Kyle Shannon <kyle at pobox dot com>
 *
 * 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 "gdal_pam.h"
#include <ctype.h>
#include <limits.h>
#include "cpl_string.h"
#include "ogr_spatialref.h"

CPL_CVSID("$Id: aaigriddataset.cpp 27729 2014-09-24 00:40:16Z goatbar $");

CPL_C_START
void    GDALRegister_AAIGrid(void);
void    GDALRegister_GRASSASCIIGrid(void);
CPL_C_END

static CPLString OSR_GDS( char **papszNV, const char * pszField, 
                           const char *pszDefaultValue );

typedef enum
{
    FORMAT_AAIG,
    FORMAT_GRASSASCII
} GridFormat;

/************************************************************************/
/* ==================================================================== */
/*                              AAIGDataset                             */
/* ==================================================================== */
/************************************************************************/

class AAIGRasterBand;

class CPL_DLL AAIGDataset : public GDALPamDataset
{
    friend class AAIGRasterBand;
    
    VSILFILE   *fp;

char        **papszPrj;
    CPLString   osPrjFilename;
    char        *pszProjection;

unsigned char achReadBuf[256];
    GUIntBig    nBufferOffset;
    int         nOffsetInBuffer;

char        Getc();
    GUIntBig    Tell();
    int         Seek( GUIntBig nOffset );

protected:
    GDALDataType eDataType;
    double      adfGeoTransform[6];
    int         bNoDataSet;
    double      dfNoDataValue;

virtual int ParseHeader(const char* pszHeader, const char* pszDataType);

public:
                AAIGDataset();
                ~AAIGDataset();

virtual char **GetFileList(void);

static GDALDataset *CommonOpen( GDALOpenInfo * poOpenInfo,
                                    GridFormat eFormat );

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

virtual CPLErr GetGeoTransform( double * );
    virtual const char *GetProjectionRef(void);
};

/************************************************************************/
/* ==================================================================== */
/*                        GRASSASCIIDataset                             */
/* ==================================================================== */
/************************************************************************/

class GRASSASCIIDataset : public AAIGDataset
{
    virtual int ParseHeader(const char* pszHeader, const char* pszDataType);

public:
                GRASSASCIIDataset() : AAIGDataset() {}

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

/************************************************************************/
/* ==================================================================== */
/*                            AAIGRasterBand                             */
/* ==================================================================== */
/************************************************************************/

class AAIGRasterBand : public GDALPamRasterBand
{
    friend class AAIGDataset;

GUIntBig      *panLineOffset;

public:

AAIGRasterBand( AAIGDataset *, int );
    virtual       ~AAIGRasterBand();

virtual double GetNoDataValue( int * );
    virtual CPLErr SetNoDataValue( double );
    virtual CPLErr IReadBlock( int, int, void * );
};

/************************************************************************/
/*                           AAIGRasterBand()                            */
/************************************************************************/

AAIGRasterBand::AAIGRasterBand( AAIGDataset *poDS, int nDataStart )

{
    this->poDS = poDS;

nBand = 1;
    eDataType = poDS->eDataType;

nBlockXSize = poDS->nRasterXSize;
    nBlockYSize = 1;

panLineOffset = 
        (GUIntBig *) VSICalloc( poDS->nRasterYSize, sizeof(GUIntBig) );
    if (panLineOffset == NULL)
    {
        CPLError(CE_Failure, CPLE_OutOfMemory,
                 "AAIGRasterBand::AAIGRasterBand : Out of memory (nRasterYSize = %d)",
                 poDS->nRasterYSize);
        return;
    }
    panLineOffset[0] = nDataStart;
}

/************************************************************************/
/*                          ~AAIGRasterBand()                           */
/************************************************************************/

AAIGRasterBand::~AAIGRasterBand()

{
    CPLFree( panLineOffset );
}

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

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

{
    AAIGDataset *poODS = (AAIGDataset *) poDS;
    int         iPixel;

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

if( panLineOffset[nBlockYOff] == 0 )
    {
        int iPrevLine;

for( iPrevLine = 1; iPrevLine <= nBlockYOff; iPrevLine++ )
            if( panLineOffset[iPrevLine] == 0 )
                IReadBlock( nBlockXOff, iPrevLine-1, NULL );
    }

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

if( poODS->Seek( panLineOffset[nBlockYOff] ) != 0 )
    {
        CPLError( CE_Failure, CPLE_FileIO,
                  "Can't seek to offset %lu in input file to read data.",
                  (long unsigned int)panLineOffset[nBlockYOff] );
        return CE_Failure;
    }

for( iPixel = 0; iPixel < poODS->nRasterXSize; )
    {
        char szToken[500];
        char chNext;
        int  iTokenChar = 0;

/* suck up any pre-white space. */
        do {
            chNext = poODS->Getc();
        } while( isspace( (unsigned char)chNext ) );

while( chNext != '\0' && !isspace((unsigned char)chNext)  )
        {
            if( iTokenChar == sizeof(szToken)-2 )
            {
                CPLError( CE_Failure, CPLE_FileIO, 
                          "Token too long at scanline %d.", 
                          nBlockYOff );
                return CE_Failure;
            }

szToken[iTokenChar++] = chNext;
            chNext = poODS->Getc();
        }

if( chNext == '\0' &&
            (iPixel != poODS->nRasterXSize - 1 ||
            nBlockYOff != poODS->nRasterYSize - 1) )
        {
            CPLError( CE_Failure, CPLE_FileIO, 
                      "File short, can't read line %d.",
                      nBlockYOff );
            return CE_Failure;
        }

szToken[iTokenChar] = '\0';

if( pImage != NULL )
        {
            if( eDataType == GDT_Float64 )
                ((double *) pImage)[iPixel] = CPLAtofM(szToken);
            else if( eDataType == GDT_Float32 )
                ((float *) pImage)[iPixel] = (float) CPLAtofM(szToken);
            else
                ((GInt32 *) pImage)[iPixel] = (GInt32) atoi(szToken);
        }
        
        iPixel++;
    }
    
    if( nBlockYOff < poODS->nRasterYSize - 1 )
        panLineOffset[nBlockYOff + 1] = poODS->Tell();

return CE_None;
}

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

double AAIGRasterBand::GetNoDataValue( int * pbSuccess )

{
    AAIGDataset *poODS = (AAIGDataset *) poDS;

if( pbSuccess )
        *pbSuccess = poODS->bNoDataSet;

return poODS->dfNoDataValue;
}

/************************************************************************/
/*                           SetNoDataValue()                           */
/************************************************************************/

CPLErr AAIGRasterBand::SetNoDataValue( double dfNoData )

{
    AAIGDataset *poODS = (AAIGDataset *) poDS;

poODS->bNoDataSet = TRUE;
    poODS->dfNoDataValue = dfNoData;

return CE_None;
}

/************************************************************************/
/* ==================================================================== */
/*                            AAIGDataset                               */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                            AAIGDataset()                            */
/************************************************************************/

AAIGDataset::AAIGDataset()

{
    papszPrj = NULL;
    pszProjection = CPLStrdup("");
    fp = NULL;
    eDataType = GDT_Int32;
    bNoDataSet = FALSE;
    dfNoDataValue = -9999.0;

adfGeoTransform[0] = 0.0;
    adfGeoTransform[1] = 1.0;
    adfGeoTransform[2] = 0.0;
    adfGeoTransform[3] = 0.0;
    adfGeoTransform[4] = 0.0;
    adfGeoTransform[5] = 1.0;

nOffsetInBuffer = 256;
    nBufferOffset = 0;
}

/************************************************************************/
/*                           ~AAIGDataset()                            */
/************************************************************************/

AAIGDataset::~AAIGDataset()

{
    FlushCache();

if( fp != NULL )
        VSIFCloseL( fp );

CPLFree( pszProjection );
    CSLDestroy( papszPrj );
}

/************************************************************************/
/*                                Tell()                                */
/************************************************************************/

GUIntBig AAIGDataset::Tell()

{
    return nBufferOffset + nOffsetInBuffer;
}

/************************************************************************/
/*                                Seek()                                */
/************************************************************************/

int AAIGDataset::Seek( GUIntBig nNewOffset )

{
    nOffsetInBuffer = sizeof(achReadBuf);
    return VSIFSeekL( fp, nNewOffset, SEEK_SET );
}

/************************************************************************/
/*                                Getc()                                */
/*                                                                      */
/*      Read a single character from the input file (efficiently we     */
/*      hope).                                                          */
/************************************************************************/

char AAIGDataset::Getc()

{
    if( nOffsetInBuffer < (int) sizeof(achReadBuf) )
        return achReadBuf[nOffsetInBuffer++];

nBufferOffset = VSIFTellL( fp );
    int nRead = VSIFReadL( achReadBuf, 1, sizeof(achReadBuf), fp );
    unsigned int i;
    for(i=nRead;i<sizeof(achReadBuf);i++)
        achReadBuf[i] = '\0';

nOffsetInBuffer = 0;

return achReadBuf[nOffsetInBuffer++];
}

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

char **AAIGDataset::GetFileList()

{
    char **papszFileList = GDALPamDataset::GetFileList();

if( papszPrj != NULL )
        papszFileList = CSLAddString( papszFileList, osPrjFilename );

return papszFileList;
}

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

int AAIGDataset::Identify( GDALOpenInfo * poOpenInfo )

{
/* -------------------------------------------------------------------- */
/*      Does this look like an AI grid file?                            */
/* -------------------------------------------------------------------- */
    if( poOpenInfo->nHeaderBytes < 40
        || !( EQUALN((const char *) poOpenInfo->pabyHeader,"ncols",5) ||
              EQUALN((const char *) poOpenInfo->pabyHeader,"nrows",5) ||
              EQUALN((const char *) poOpenInfo->pabyHeader,"xllcorner",9)||
              EQUALN((const char *) poOpenInfo->pabyHeader,"yllcorner",9)||
              EQUALN((const char *) poOpenInfo->pabyHeader,"xllcenter",9)||
              EQUALN((const char *) poOpenInfo->pabyHeader,"yllcenter",9)||
              EQUALN((const char *) poOpenInfo->pabyHeader,"dx",2)||
              EQUALN((const char *) poOpenInfo->pabyHeader,"dy",2)||
              EQUALN((const char *) poOpenInfo->pabyHeader,"cellsize",8)) )
        return FALSE;

return TRUE;
}

int GRASSASCIIDataset::Identify( GDALOpenInfo * poOpenInfo )

{
/* -------------------------------------------------------------------- */
/*      Does this look like a GRASS ASCII grid file?                    */
/* -------------------------------------------------------------------- */
    if( poOpenInfo->nHeaderBytes < 40
        || !( EQUALN((const char *) poOpenInfo->pabyHeader,"north:",6) ||
              EQUALN((const char *) poOpenInfo->pabyHeader,"south:",6) ||
              EQUALN((const char *) poOpenInfo->pabyHeader,"east:",5)||
              EQUALN((const char *) poOpenInfo->pabyHeader,"west:",5)||
              EQUALN((const char *) poOpenInfo->pabyHeader,"rows:",5)||
              EQUALN((const char *) poOpenInfo->pabyHeader,"cols:",5) ) )
        return FALSE;

return TRUE;
}

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

GDALDataset *AAIGDataset::Open( GDALOpenInfo * poOpenInfo )
{
    if (!Identify(poOpenInfo))
        return NULL;

return CommonOpen(poOpenInfo, FORMAT_AAIG);
}

/************************************************************************/
/*                          ParseHeader()                               */
/************************************************************************/

int AAIGDataset::ParseHeader(const char* pszHeader, const char* pszDataType)
{
    int i, j;
    char** papszTokens =
        CSLTokenizeString2( pszHeader,  " \n\r\t" , 0 );
    int nTokens = CSLCount(papszTokens);
    double dfCellDX = 0;
    double dfCellDY = 0;

if ( (i = CSLFindString( papszTokens, "ncols" )) < 0 ||
         i + 1 >= nTokens)
    {
        CSLDestroy( papszTokens );
        return FALSE;
    }
    nRasterXSize = atoi(papszTokens[i + 1]);
    if ( (i = CSLFindString( papszTokens, "nrows" )) < 0 ||
         i + 1 >= nTokens)
    {
        CSLDestroy( papszTokens );
        return FALSE;
    }
    nRasterYSize = atoi(papszTokens[i + 1]);

if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
    {
        CSLDestroy( papszTokens );
        return FALSE;
    }

if ( (i = CSLFindString( papszTokens, "cellsize" )) < 0 )
    {
        int iDX, iDY;
        if( (iDX = CSLFindString(papszTokens,"dx")) < 0
            || (iDY = CSLFindString(papszTokens,"dy")) < 0
            || iDX+1 >= nTokens
            || iDY+1 >= nTokens)
        {
            CSLDestroy( papszTokens );
            return FALSE;
        }

dfCellDX = CPLAtofM( papszTokens[iDX+1] );
        dfCellDY = CPLAtofM( papszTokens[iDY+1] );
    }
    else
    {
        if (i + 1 >= nTokens)
        {
            CSLDestroy( papszTokens );
            return FALSE;
        }
        dfCellDX = dfCellDY = CPLAtofM( papszTokens[i + 1] );
    }

if ((i = CSLFindString( papszTokens, "xllcorner" )) >= 0 &&
        (j = CSLFindString( papszTokens, "yllcorner" )) >= 0 &&
        i + 1 < nTokens && j + 1 < nTokens)
    {
        adfGeoTransform[0] = CPLAtofM( papszTokens[i + 1] );
       
        /* Small hack to compensate from insufficient precision in cellsize */
        /* parameter in datasets of http://ccafs-climate.org/data/A2a_2020s/hccpr_hadcm3 */
        if ((nRasterXSize % 360) == 0 &&
            fabs(adfGeoTransform[0] - (-180.0)) < 1e-12 &&
            dfCellDX == dfCellDY &&
            fabs(dfCellDX - (360.0 / nRasterXSize)) < 1e-9)
        {
            dfCellDX = dfCellDY = 360.0 / nRasterXSize;
        }
            
        adfGeoTransform[1] = dfCellDX;
        adfGeoTransform[2] = 0.0;
        adfGeoTransform[3] = CPLAtofM( papszTokens[j + 1] )
            + nRasterYSize * dfCellDY;
        adfGeoTransform[4] = 0.0;
        adfGeoTransform[5] = - dfCellDY;
    }
    else if ((i = CSLFindString( papszTokens, "xllcenter" )) >= 0 &&
             (j = CSLFindString( papszTokens, "yllcenter" )) >= 0  &&
             i + 1 < nTokens && j + 1 < nTokens)
    {
        SetMetadataItem( GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT );

adfGeoTransform[0] = CPLAtofM(papszTokens[i + 1]) - 0.5 * dfCellDX;
        adfGeoTransform[1] = dfCellDX;
        adfGeoTransform[2] = 0.0;
        adfGeoTransform[3] = CPLAtofM( papszTokens[j + 1] )
            - 0.5 * dfCellDY
            + nRasterYSize * dfCellDY;
        adfGeoTransform[4] = 0.0;
        adfGeoTransform[5] = - dfCellDY;
    }
    else
    {
        adfGeoTransform[0] = 0.0;
        adfGeoTransform[1] = dfCellDX;
        adfGeoTransform[2] = 0.0;
        adfGeoTransform[3] = 0.0;
        adfGeoTransform[4] = 0.0;
        adfGeoTransform[5] = - dfCellDY;
    }

if( (i = CSLFindString( papszTokens, "NODATA_value" )) >= 0 &&
        i + 1 < nTokens)
    {
        const char* pszNoData = papszTokens[i + 1];

bNoDataSet = TRUE;
        dfNoDataValue = CPLAtofM(pszNoData);
        if( pszDataType == NULL &&
            (strchr( pszNoData, '.' ) != NULL ||
             strchr( pszNoData, ',' ) != NULL ||
             INT_MIN > dfNoDataValue || dfNoDataValue > INT_MAX) )
        {
            eDataType = GDT_Float32;
        }
        if( eDataType == GDT_Float32 )
        {
            dfNoDataValue = (double) (float) dfNoDataValue;
        }
    }

CSLDestroy( papszTokens );

return TRUE;
}

GDALDataset *GRASSASCIIDataset::Open( GDALOpenInfo * poOpenInfo )
{
    if (!Identify(poOpenInfo))
        return NULL;

return CommonOpen(poOpenInfo, FORMAT_GRASSASCII);
}

int GRASSASCIIDataset::ParseHeader(const char* pszHeader, const char* pszDataType)
{
    int i;
    char** papszTokens =
        CSLTokenizeString2( pszHeader,  " \n\r\t:" , 0 );
    int nTokens = CSLCount(papszTokens);
    if ( (i = CSLFindString( papszTokens, "cols" )) < 0 ||
         i + 1 >= nTokens)
    {
        CSLDestroy( papszTokens );
        return FALSE;
    }
    nRasterXSize = atoi(papszTokens[i + 1]);
    if ( (i = CSLFindString( papszTokens, "rows" )) < 0 ||
         i + 1 >= nTokens)
    {
        CSLDestroy( papszTokens );
        return FALSE;
    }
    nRasterYSize = atoi(papszTokens[i + 1]);

if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
    {
        CSLDestroy( papszTokens );
        return FALSE;
    }

int iNorth = CSLFindString( papszTokens, "north" );
    int iSouth = CSLFindString( papszTokens, "south" );
    int iEast = CSLFindString( papszTokens, "east" );
    int iWest = CSLFindString( papszTokens, "west" );

if (iNorth == -1 || iSouth == -1 || iEast == -1 || iWest == -1 ||
        MAX(MAX(iNorth, iSouth), MAX(iEast, iWest)) + 1 >= nTokens)
    {
        CSLDestroy( papszTokens );
        return FALSE;
    }

double dfNorth = CPLAtofM( papszTokens[iNorth + 1] );
    double dfSouth = CPLAtofM( papszTokens[iSouth + 1] );
    double dfEast = CPLAtofM( papszTokens[iEast + 1] );
    double dfWest = CPLAtofM( papszTokens[iWest + 1] );
    double dfPixelXSize = (dfEast - dfWest) / nRasterXSize;
    double dfPixelYSize = (dfNorth - dfSouth) / nRasterYSize;

adfGeoTransform[0] = dfWest;
    adfGeoTransform[1] = dfPixelXSize;
    adfGeoTransform[2] = 0.0;
    adfGeoTransform[3] = dfNorth;
    adfGeoTransform[4] = 0.0;
    adfGeoTransform[5] = - dfPixelYSize;

if( (i = CSLFindString( papszTokens, "null" )) >= 0 &&
        i + 1 < nTokens)
    {
        const char* pszNoData = papszTokens[i + 1];

if( (i = CSLFindString( papszTokens, "type" )) >= 0 &&
        i + 1 < nTokens)
    {
        const char* pszType = papszTokens[i + 1];
        if (EQUAL(pszType, "int"))
            eDataType = GDT_Int32;
        else if (EQUAL(pszType, "float"))
            eDataType = GDT_Float32;
        else if (EQUAL(pszType, "double"))
            eDataType = GDT_Float64;
        else
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                     "Invalid value for type parameter : %s", pszType);
        }
    }

CSLDestroy(papszTokens);

return TRUE;
}

/************************************************************************/
/*                           CommonOpen()                               */
/************************************************************************/

GDALDataset *AAIGDataset::CommonOpen( GDALOpenInfo * poOpenInfo,
                                            GridFormat eFormat )
{
    int i = 0;

/* -------------------------------------------------------------------- */
/*      Create a corresponding GDALDataset.                             */
/* -------------------------------------------------------------------- */
    AAIGDataset         *poDS;

if (eFormat == FORMAT_AAIG)
        poDS = new AAIGDataset();
    else
        poDS = new GRASSASCIIDataset();

const char* pszDataTypeOption = (eFormat == FORMAT_AAIG) ? "AAIGRID_DATATYPE":
                                                               "GRASSASCIIGRID_DATATYPE";
    const char* pszDataType = CPLGetConfigOption(pszDataTypeOption, NULL);
    if (pszDataType != NULL)
    {
        poDS->eDataType = GDALGetDataTypeByName(pszDataType);
        if (!(poDS->eDataType == GDT_Int32 || poDS->eDataType == GDT_Float32 ||
              poDS->eDataType == GDT_Float64))
        {
            CPLError(CE_Warning, CPLE_NotSupported,
                     "Unsupported value for %s : %s", pszDataTypeOption, pszDataType);
            poDS->eDataType = GDT_Int32;
            pszDataType = NULL;
        }
    }
    
/* -------------------------------------------------------------------- */
/*      Parse the header.                                               */
/* -------------------------------------------------------------------- */
    if (!poDS->ParseHeader((const char *) poOpenInfo->pabyHeader, pszDataType))
    {
        delete poDS;
        return NULL;
    }

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

poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "r" );
    if( poDS->fp == NULL )
    {
        CPLError( CE_Failure, CPLE_OpenFailed, 
                  "VSIFOpenL(%s) failed unexpectedly.", 
                  poOpenInfo->pszFilename );
        delete poDS;
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Find the start of real data.                                    */
/* -------------------------------------------------------------------- */
    int         nStartOfData;

for( i = 2; TRUE ; i++ )
    {
        if( poOpenInfo->pabyHeader[i] == '\0' )
        {
            CPLError( CE_Failure, CPLE_AppDefined, 
                      "Couldn't find data values in ASCII Grid file.\n" );
            delete poDS;
            return NULL;
        }
        
        if( poOpenInfo->pabyHeader[i-1] == '\n' 
            || poOpenInfo->pabyHeader[i-2] == '\n' 
            || poOpenInfo->pabyHeader[i-1] == '\r' 
            || poOpenInfo->pabyHeader[i-2] == '\r' )
        {
            if( !isalpha(poOpenInfo->pabyHeader[i]) 
                && poOpenInfo->pabyHeader[i] != '\n' 
                && poOpenInfo->pabyHeader[i] != '\r')
            {
                nStartOfData = i;

/* Beginning of real data found. */
                break;
            }
        }
    }

/* -------------------------------------------------------------------- */
/*      Recognize the type of data.										*/
/* -------------------------------------------------------------------- */
    CPLAssert( NULL != poDS->fp );

if( pszDataType == NULL && poDS->eDataType != GDT_Float32)
    {
        /* Allocate 100K chunk + 1 extra byte for NULL character. */
        const size_t nChunkSize = 1024 * 100;
        GByte* pabyChunk = (GByte *) VSICalloc( nChunkSize + 1, sizeof(GByte) );
        if (pabyChunk == NULL)
        {
            CPLError( CE_Failure, CPLE_OutOfMemory, "Out of memory");
            delete poDS;
            return NULL;
        }
        pabyChunk[nChunkSize] = '\0';

VSIFSeekL( poDS->fp, nStartOfData, SEEK_SET );

/* Scan for dot in subsequent chunks of data. */
        while( !VSIFEofL( poDS->fp) )
        {
            VSIFReadL( pabyChunk, sizeof(GByte), nChunkSize, poDS->fp );
            CPLAssert( pabyChunk[nChunkSize] == '\0' );

for(i = 0; i < (int)nChunkSize; i++)
            {
                GByte ch = pabyChunk[i];
                if (ch == '.' || ch == ',' || ch == 'e' || ch == 'E')
                {
                    poDS->eDataType = GDT_Float32;
                    break;
                }
            }
        }

/* Deallocate chunk. */
        VSIFree( pabyChunk );
    }

/* -------------------------------------------------------------------- */
/*      Create band information objects.                                */
/* -------------------------------------------------------------------- */
    AAIGRasterBand* band = new AAIGRasterBand( poDS, nStartOfData );
    poDS->SetBand( 1, band );
    if (band->panLineOffset == NULL)
    {
        delete poDS;
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Try to read projection file.                                    */
/* -------------------------------------------------------------------- */
    char        *pszDirname, *pszBasename;
    VSIStatBufL   sStatBuf;

pszDirname = CPLStrdup(CPLGetPath(poOpenInfo->pszFilename));
    pszBasename = CPLStrdup(CPLGetBasename(poOpenInfo->pszFilename));

poDS->osPrjFilename = CPLFormFilename( pszDirname, pszBasename, "prj" );
    int nRet = VSIStatL( poDS->osPrjFilename, &sStatBuf );

if( nRet != 0 && VSIIsCaseSensitiveFS(poDS->osPrjFilename) )
    {
        poDS->osPrjFilename = CPLFormFilename( pszDirname, pszBasename, "PRJ" );
        nRet = VSIStatL( poDS->osPrjFilename, &sStatBuf );
    }

if( nRet == 0 )
    {
        OGRSpatialReference     oSRS;

poDS->papszPrj = CSLLoad( poDS->osPrjFilename );

CPLDebug( "AAIGrid", "Loaded SRS from %s", 
                  poDS->osPrjFilename.c_str() );

if( oSRS.importFromESRI( poDS->papszPrj ) == OGRERR_NONE )
        {
            // If geographic values are in seconds, we must transform. 
            // Is there a code for minutes too? 
            if( oSRS.IsGeographic() 
                && EQUAL(OSR_GDS( poDS->papszPrj, "Units", ""), "DS") )
            {
                poDS->adfGeoTransform[0] /= 3600.0;
                poDS->adfGeoTransform[1] /= 3600.0;
                poDS->adfGeoTransform[2] /= 3600.0;
                poDS->adfGeoTransform[3] /= 3600.0;
                poDS->adfGeoTransform[4] /= 3600.0;
                poDS->adfGeoTransform[5] /= 3600.0;
            }

CPLFree( poDS->pszProjection );
            oSRS.exportToWkt( &(poDS->pszProjection) );
        }
    }

CPLFree( pszDirname );
    CPLFree( pszBasename );

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

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

return( poDS );
}

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

CPLErr AAIGDataset::GetGeoTransform( double * padfTransform )

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

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

const char *AAIGDataset::GetProjectionRef()

{
    return pszProjection;
}

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

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

{
    int  nBands = poSrcDS->GetRasterCount();
    int  nXSize = poSrcDS->GetRasterXSize();
    int  nYSize = poSrcDS->GetRasterYSize();

/* -------------------------------------------------------------------- */
/*      Some rudimentary checks                                         */
/* -------------------------------------------------------------------- */
    if( nBands != 1 )
    {
        CPLError( CE_Failure, CPLE_NotSupported, 
                  "AAIG driver doesn't support %d bands.  Must be 1 band.\n",
                  nBands );

return NULL;
    }

if( !pfnProgress( 0.0, NULL, pProgressData ) )
        return NULL;

/* -------------------------------------------------------------------- */
/*      Create the dataset.                                             */
/* -------------------------------------------------------------------- */
    VSILFILE        *fpImage;

fpImage = VSIFOpenL( pszFilename, "wt" );
    if( fpImage == NULL )
    {
        CPLError( CE_Failure, CPLE_OpenFailed, 
                  "Unable to create file %s.\n", 
                  pszFilename );
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Write ASCII Grid file header                                    */
/* -------------------------------------------------------------------- */
    double      adfGeoTransform[6];
    char        szHeader[2000];
    const char *pszForceCellsize = 
        CSLFetchNameValue( papszOptions, "FORCE_CELLSIZE" );

poSrcDS->GetGeoTransform( adfGeoTransform );

if( ABS(adfGeoTransform[1]+adfGeoTransform[5]) < 0.0000001 
        || ABS(adfGeoTransform[1]-adfGeoTransform[5]) < 0.0000001 
        || (pszForceCellsize && CSLTestBoolean(pszForceCellsize)) )
        sprintf( szHeader, 
                 "ncols        %d\n" 
                 "nrows        %d\n"
                 "xllcorner    %.12f\n"
                 "yllcorner    %.12f\n"
                 "cellsize     %.12f\n", 
                 nXSize, nYSize, 
                 adfGeoTransform[0], 
                 adfGeoTransform[3] - nYSize * adfGeoTransform[1],
                 adfGeoTransform[1] );
    else
    {
        if( pszForceCellsize == NULL )
            CPLError( CE_Warning, CPLE_AppDefined, 
                      "Producing a Golden Surfer style file with DX and DY instead\n"
                      "of CELLSIZE since the input pixels are non-square.  Use the\n"
                      "FORCE_CELLSIZE=TRUE creation option to force use of DX for\n"
                      "even though this will be distorted.  Most ASCII Grid readers\n"
                      "(ArcGIS included) do not support the DX and DY parameters.\n" );
        sprintf( szHeader, 
                 "ncols        %d\n" 
                 "nrows        %d\n"
                 "xllcorner    %.12f\n"
                 "yllcorner    %.12f\n"
                 "dx           %.12f\n"
                 "dy           %.12f\n", 
                 nXSize, nYSize, 
                 adfGeoTransform[0], 
                 adfGeoTransform[3] + nYSize * adfGeoTransform[5],
                 adfGeoTransform[1],
                 fabs(adfGeoTransform[5]) );
    }

/* -------------------------------------------------------------------- */
/*     Builds the format string used for printing float values.         */
/* -------------------------------------------------------------------- */
    char szFormatFloat[32];
    strcpy(szFormatFloat, " %.20g");
    const char *pszDecimalPrecision = 
        CSLFetchNameValue( papszOptions, "DECIMAL_PRECISION" );
    const char *pszSignificantDigits =
        CSLFetchNameValue( papszOptions, "SIGNIFICANT_DIGITS" );
    int bIgnoreSigDigits = FALSE;
    if( pszDecimalPrecision && pszSignificantDigits )
    {
        CPLError( CE_Warning, CPLE_AppDefined,
                  "Conflicting precision arguments, using DECIMAL_PRECISION" );
        bIgnoreSigDigits = TRUE;
    }
    int nPrecision;
    if ( pszSignificantDigits && !bIgnoreSigDigits )
    {
        nPrecision = atoi( pszSignificantDigits );
        if (nPrecision >= 0)
            sprintf( szFormatFloat, " %%.%dg", nPrecision );
        CPLDebug( "AAIGrid", "Setting precision format: %s", szFormatFloat );
    }
    else if( pszDecimalPrecision )
    {
        nPrecision = atoi( pszDecimalPrecision );
        if ( nPrecision >= 0 )
            sprintf( szFormatFloat, " %%.%dlf", nPrecision );
        CPLDebug( "AAIGrid", "Setting precision format: %s", szFormatFloat );
    }

/* -------------------------------------------------------------------- */
/*      Handle nodata (optionally).                                     */
/* -------------------------------------------------------------------- */
    GDALRasterBand * poBand = poSrcDS->GetRasterBand( 1 );
    double dfNoData;
    int bSuccess;
    int bReadAsInt;
    bReadAsInt = ( poBand->GetRasterDataType() == GDT_Byte
                || poBand->GetRasterDataType() == GDT_Int16
                || poBand->GetRasterDataType() == GDT_UInt16
                || poBand->GetRasterDataType() == GDT_Int32 );

// Write `nodata' value to header if it is exists in source dataset
    dfNoData = poBand->GetNoDataValue( &bSuccess );
    if ( bSuccess )
    {
        sprintf( szHeader+strlen( szHeader ), "NODATA_value " );
        if( bReadAsInt )
            sprintf( szHeader+strlen( szHeader ), "%d", (int)dfNoData );
        else
            sprintf( szHeader+strlen( szHeader ), szFormatFloat, dfNoData );
        sprintf( szHeader+strlen( szHeader ), "\n" );
    }

VSIFWriteL( szHeader, 1, strlen(szHeader), fpImage );

/* -------------------------------------------------------------------- */
/*      Loop over image, copying image data.                            */
/* -------------------------------------------------------------------- */
    int         *panScanline = NULL;
    double      *padfScanline = NULL;
    int         iLine, iPixel;
    CPLErr      eErr = CE_None;

// Write scanlines to output file
    if (bReadAsInt)
        panScanline = (int *) CPLMalloc( nXSize *
                                GDALGetDataTypeSize(GDT_Int32) / 8 );
    else
        padfScanline = (double *) CPLMalloc( nXSize *
                                    GDALGetDataTypeSize(GDT_Float64) / 8 );

for( iLine = 0; eErr == CE_None && iLine < nYSize; iLine++ )
    {
        CPLString osBuf;
        eErr = poBand->RasterIO( GF_Read, 0, iLine, nXSize, 1, 
                                 (bReadAsInt) ? (void*)panScanline : (void*)padfScanline,
                                 nXSize, 1, (bReadAsInt) ? GDT_Int32 : GDT_Float64,
                                 0, 0 );

if( bReadAsInt )
        {
            for ( iPixel = 0; iPixel < nXSize; iPixel++ )
            {
                sprintf( szHeader, " %d", panScanline[iPixel] );
                osBuf += szHeader;
                if( (iPixel & 1023) == 0 || iPixel == nXSize - 1 )
                {
                    if ( VSIFWriteL( osBuf, (int)osBuf.size(), 1, fpImage ) != 1 )
                    {
                        eErr = CE_Failure;
                        CPLError( CE_Failure, CPLE_AppDefined, 
                                  "Write failed, disk full?\n" );
                        break;
                    }
                    osBuf = "";
                }
            }
        }
        else
        {
            for ( iPixel = 0; iPixel < nXSize; iPixel++ )
            {
                sprintf( szHeader, szFormatFloat, padfScanline[iPixel] );
                osBuf += szHeader;
                if( (iPixel & 1023) == 0 || iPixel == nXSize - 1 )
                {
                    if ( VSIFWriteL( osBuf, (int)osBuf.size(), 1, fpImage ) != 1 )
                    {
                        eErr = CE_Failure;
                        CPLError( CE_Failure, CPLE_AppDefined, 
                                  "Write failed, disk full?\n" );
                        break;
                    }
                    osBuf = "";
                }
            }
        }
        VSIFWriteL( (void *) "\n", 1, 1, fpImage );

if( eErr == CE_None &&
            !pfnProgress((iLine + 1) / ((double) nYSize), NULL, pProgressData) )
        {
            eErr = CE_Failure;
            CPLError( CE_Failure, CPLE_UserInterrupt, 
                      "User terminated CreateCopy()" );
        }
    }

CPLFree( panScanline );
    CPLFree( padfScanline );
    VSIFCloseL( fpImage );

if( eErr != CE_None )
        return NULL;

/* -------------------------------------------------------------------- */
/*      Try to write projection file.                                   */
/* -------------------------------------------------------------------- */
    const char  *pszOriginalProjection;

pszOriginalProjection = (char *)poSrcDS->GetProjectionRef();
    if( !EQUAL( pszOriginalProjection, "" ) )
    {
        char                    *pszDirname, *pszBasename;
        char                    *pszPrjFilename;
        char                    *pszESRIProjection = NULL;
        VSILFILE                *fp;
        OGRSpatialReference     oSRS;

pszDirname = CPLStrdup( CPLGetPath(pszFilename) );
        pszBasename = CPLStrdup( CPLGetBasename(pszFilename) );

pszPrjFilename = CPLStrdup( CPLFormFilename( pszDirname, pszBasename, "prj" ) );
        fp = VSIFOpenL( pszPrjFilename, "wt" );
        if (fp != NULL)
        {
            oSRS.importFromWkt( (char **) &pszOriginalProjection );
            oSRS.morphToESRI();
            oSRS.exportToWkt( &pszESRIProjection );
            VSIFWriteL( pszESRIProjection, 1, strlen(pszESRIProjection), fp );

VSIFCloseL( fp );
            CPLFree( pszESRIProjection );
        }
        else
        {
            CPLError( CE_Failure, CPLE_FileIO, 
                      "Unable to create file %s.\n", pszPrjFilename );
        }
        CPLFree( pszDirname );
        CPLFree( pszBasename );
        CPLFree( pszPrjFilename );
    }
    
/* -------------------------------------------------------------------- */
/*      Re-open dataset, and copy any auxilary pam information.         */
/* -------------------------------------------------------------------- */

/* If outputing to stdout, we can't reopen it, so we'll return */
    /* a fake dataset to make the caller happy */
    CPLPushErrorHandler(CPLQuietErrorHandler);
    GDALPamDataset* poDS = (GDALPamDataset*) GDALOpen(pszFilename, GA_ReadOnly);
    CPLPopErrorHandler();
    if (poDS)
    {
        poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
        return poDS;
    }

CPLErrorReset();

AAIGDataset* poAAIG_DS = new AAIGDataset();
    poAAIG_DS->nRasterXSize = nXSize;
    poAAIG_DS->nRasterYSize = nYSize;
    poAAIG_DS->nBands = 1;
    poAAIG_DS->SetBand( 1, new AAIGRasterBand( poAAIG_DS, 1 ) );
    return poAAIG_DS;
}

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

static CPLString OSR_GDS( char **papszNV, const char * pszField, 
                           const char *pszDefaultValue )

{
    int         iLine;

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

for( iLine = 0; 
         papszNV[iLine] != NULL && 
             !EQUALN(papszNV[iLine],pszField,strlen(pszField));
         iLine++ ) {}

if( papszNV[iLine] == NULL )
        return pszDefaultValue;
    else
    {
        CPLString osResult;
        char    **papszTokens;
        
        papszTokens = CSLTokenizeString(papszNV[iLine]);

if( CSLCount(papszTokens) > 1 )
            osResult = papszTokens[1];
        else
            osResult = pszDefaultValue;
        
        CSLDestroy( papszTokens );
        return osResult;
    }
}

/************************************************************************/
/*                        GDALRegister_AAIGrid()                        */
/************************************************************************/

void GDALRegister_AAIGrid()

{
    GDALDriver  *poDriver;

if( GDALGetDriverByName( "AAIGrid" ) == NULL )
    {
        poDriver = new GDALDriver();
        
        poDriver->SetDescription( "AAIGrid" );
        poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, 
                                   "Arc/Info ASCII Grid" );
        poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, 
                                   "frmt_various.html#AAIGrid" );
        poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "asc" );
        poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, 
                                   "Byte UInt16 Int16 Int32 Float32" );

poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
        poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST, 
"<CreationOptionList>\n"
"   <Option name='FORCE_CELLSIZE' type='boolean' description='Force use of CELLSIZE, default is FALSE.'/>\n"
"   <Option name='DECIMAL_PRECISION' type='int' description='Number of decimal when writing floating-point numbers(%f).'/>\n"
"   <Option name='SIGNIFICANT_DIGITS' type='int' description='Number of significant digits when writing floating-point numbers(%g).'/>\n"
"</CreationOptionList>\n" );

poDriver->pfnOpen = AAIGDataset::Open;
        poDriver->pfnIdentify = AAIGDataset::Identify;
        poDriver->pfnCreateCopy = AAIGDataset::CreateCopy;
        
        GetGDALDriverManager()->RegisterDriver( poDriver );
    }
}

/************************************************************************/
/*                   GDALRegister_GRASSASCIIGrid()                      */
/************************************************************************/

void GDALRegister_GRASSASCIIGrid()

{
    GDALDriver  *poDriver;

if( GDALGetDriverByName( "GRASSASCIIGrid" ) == NULL )
    {
        poDriver = new GDALDriver();

poDriver->SetDescription( "GRASSASCIIGrid" );
        poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
                                   "GRASS ASCII Grid" );
        poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
                                   "frmt_various.html#GRASSASCIIGrid" );

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

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

GetGDALDriverManager()->RegisterDriver( poDriver );
    }
}