← Back to branch summary

~ubuntu-branches/debian/sid/gdal/sid

« back to all changes in this revision

Viewing changes to frmts/e00grid/e00griddataset.cpp

  • Committer: Package Import Robot
  • Author(s): Francesco Paolo Lovergine
  • Date: 2012-05-07 15:04:42 UTC
  • mfrom: (5.5.16 experimental)
  • Revision ID: package-import@ubuntu.com-20120507150442-2eks97loeh6rq005
Tags: 1.9.0-1
* Ready for sid, starting transition.
* All symfiles updated to latest builds.
* Added dh_numpy call in debian/rules to depend on numpy ABI.
* Policy bumped to 3.9.3, no changes required.

Show diffs side-by-side

added added

removed removed

Lines of Context:
frmts/e00grid/e00griddataset.cpp
 
1
/******************************************************************************
 
2
 * $Id: e00griddataset.cpp 23618 2011-12-20 22:27:21Z rouault $
 
3
 *
 
4
 * Project:  E00 grid driver
 
5
 * Purpose:  GDALDataset driver for E00 grid dataset.
 
6
 * Author:   Even Rouault, <even dot rouault at mines dash paris dot org>
 
7
 *
 
8
 ******************************************************************************
 
9
 * Copyright (c) 2011, Even Rouault
 
10
 *
 
11
 * Permission is hereby granted, free of charge, to any person obtaining a
 
12
 * copy of this software and associated documentation files (the "Software"),
 
13
 * to deal in the Software without restriction, including without limitation
 
14
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 
15
 * and/or sell copies of the Software, and to permit persons to whom the
 
16
 * Software is furnished to do so, subject to the following conditions:
 
17
 *
 
18
 * The above copyright notice and this permission notice shall be included
 
19
 * in all copies or substantial portions of the Software.
 
20
 *
 
21
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 
22
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 
23
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 
24
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 
25
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 
26
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 
27
 * DEALINGS IN THE SOFTWARE.
 
28
 ****************************************************************************/
 
29
 
 
30
#include "cpl_vsi_virtual.h"
 
31
#include "cpl_string.h"
 
32
#include "ogr_spatialref.h"
 
33
#include "gdal_pam.h"
 
34
 
 
35
/* Private import of e00read.c */
 
36
#define E00ReadOpen         GDALE00GRIDReadOpen
 
37
#define E00ReadCallbackOpen GDALE00GRIDReadCallbackOpen
 
38
#define E00ReadClose        GDALE00GRIDReadClose
 
39
#define E00ReadNextLine     GDALE00GRIDReadNextLine
 
40
#define E00ReadRewind       GDALE00GRIDReadRewind
 
41
#include "e00read.c"
 
42
 
 
43
#define E00_INT_SIZE    10
 
44
#define E00_INT14_SIZE  14
 
45
#define E00_FLOAT_SIZE  14
 
46
#define E00_DOUBLE_SIZE 21
 
47
#define VALS_PER_LINE   5
 
48
 
 
49
CPL_CVSID("$Id: e00griddataset.cpp 23618 2011-12-20 22:27:21Z rouault $");
 
50
 
 
51
CPL_C_START
 
52
void    GDALRegister_E00GRID(void);
 
53
CPL_C_END
 
54
 
 
55
/* g++ -fPIC -Wall -g frmts/e00grid/e00griddataset.cpp -shared -o gdal_E00GRID.so -Iport -Igcore -Iogr -L. -lgdal */
 
56
 
 
57
/* Test data ; (google for "EXP  0" "GRD  2")
 
58
 
 
59
ftp://msdis.missouri.edu/pub/dem/24k/county/
 
60
http://dusk.geo.orst.edu/djl/samoa/data/samoa_bathy.e00
 
61
http://dusk.geo.orst.edu/djl/samoa/FBNMS/RasterGrids-Metadata/ntae02_3m_utm.e00
 
62
http://www.navdat.org/coverages/elevation/iddem1.e00        (int32)
 
63
http://delta-vision.projects.atlas.ca.gov/lidar/bare_earth.grids/sac0165.e00
 
64
http://ag.arizona.edu/SRER/maps_e00/srer_dem.e00
 
65
http://ok.water.usgs.gov/projects/norlan/spatial/ntopo0408-10.e00 (compressed)
 
66
http://wrri.nmsu.edu/publish/techrpt/tr322/GIS/dem.e00 (compressed)
 
67
*/
 
68
 
 
69
/************************************************************************/
 
70
/* ==================================================================== */
 
71
/*                            E00GRIDDataset                            */
 
72
/* ==================================================================== */
 
73
/************************************************************************/
 
74
 
 
75
class E00GRIDRasterBand;
 
76
 
 
77
class E00GRIDDataset : public GDALPamDataset
 
78
{
 
79
    friend class E00GRIDRasterBand;
 
80
 
 
81
    E00ReadPtr  e00ReadPtr;
 
82
    VSILFILE   *fp;
 
83
    vsi_l_offset nDataStart;
 
84
    int         nBytesEOL;
 
85
 
 
86
    vsi_l_offset  nPosBeforeReadLine;
 
87
    vsi_l_offset* panOffsets;
 
88
    int         nLastYOff;
 
89
    int         nMaxYOffset;
 
90
 
 
91
    double      adfGeoTransform[6];
 
92
    CPLString   osProjection;
 
93
 
 
94
    double      dfNoData;
 
95
 
 
96
    char**      papszPrj;
 
97
 
 
98
    const char* ReadLine();
 
99
 
 
100
    int         bHasReadMetadata;
 
101
    void        ReadMetadata();
 
102
 
 
103
    int         bHasStats;
 
104
    double      dfMin, dfMax, dfMean, dfStddev;
 
105
 
 
106
    static const char* ReadNextLine(void * ptr);
 
107
    static void        Rewind(void* ptr);
 
108
 
 
109
  public:
 
110
                 E00GRIDDataset();
 
111
    virtual     ~E00GRIDDataset();
 
112
    
 
113
    virtual CPLErr GetGeoTransform( double * );
 
114
    virtual const char* GetProjectionRef();
 
115
    
 
116
    static GDALDataset *Open( GDALOpenInfo * );
 
117
    static int          Identify( GDALOpenInfo * );
 
118
};
 
119
 
 
120
/************************************************************************/
 
121
/* ==================================================================== */
 
122
/*                          E00GRIDRasterBand                           */
 
123
/* ==================================================================== */
 
124
/************************************************************************/
 
125
 
 
126
class E00GRIDRasterBand : public GDALPamRasterBand
 
127
{
 
128
    friend class E00GRIDDataset;
 
129
 
 
130
  public:
 
131
 
 
132
                E00GRIDRasterBand( E00GRIDDataset *, int, GDALDataType );
 
133
 
 
134
    virtual CPLErr      IReadBlock( int, int, void * );
 
135
 
 
136
    virtual double      GetNoDataValue( int *pbSuccess = NULL );
 
137
    virtual const char *GetUnitType();
 
138
    virtual double      GetMinimum( int *pbSuccess = NULL );
 
139
    virtual double      GetMaximum( int *pbSuccess = NULL );
 
140
    virtual CPLErr      GetStatistics( int bApproxOK, int bForce,
 
141
                                       double *pdfMin, double *pdfMax,
 
142
                                       double *pdfMean, double *padfStdDev );
 
143
};
 
144
 
 
145
 
 
146
/************************************************************************/
 
147
/*                         E00GRIDRasterBand()                          */
 
148
/************************************************************************/
 
149
 
 
150
E00GRIDRasterBand::E00GRIDRasterBand( E00GRIDDataset *poDS, int nBand,
 
151
                                      GDALDataType eDT )
 
152
 
 
153
{
 
154
    this->poDS = poDS;
 
155
    this->nBand = nBand;
 
156
 
 
157
    eDataType = eDT;
 
158
 
 
159
    nBlockXSize = poDS->GetRasterXSize();
 
160
    nBlockYSize = 1;
 
161
}
 
162
 
 
163
/************************************************************************/
 
164
/*                             IReadBlock()                             */
 
165
/************************************************************************/
 
166
 
 
167
CPLErr E00GRIDRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
 
168
                                      void * pImage )
 
169
 
 
170
{
 
171
    E00GRIDDataset *poGDS = (E00GRIDDataset *) poDS;
 
172
 
 
173
    char szVal[E00_FLOAT_SIZE+1];
 
174
    szVal[E00_FLOAT_SIZE] = 0;
 
175
 
 
176
    int i;
 
177
    float* pafImage = (float*)pImage;
 
178
    int* panImage = (int*)pImage;
 
179
    const float fNoData = (const float)poGDS->dfNoData;
 
180
 
 
181
    /* A new data line begins on a new text line. So if the xsize */
 
182
    /* is not a multiple of VALS_PER_LINE, there are padding values */
 
183
    /* that must be ignored */
 
184
    const int nRoundedBlockXSize = ((nBlockXSize + VALS_PER_LINE - 1) /
 
185
                                            VALS_PER_LINE) * VALS_PER_LINE;
 
186
 
 
187
    if (poGDS->e00ReadPtr)
 
188
    {
 
189
        if (poGDS->nLastYOff < 0)
 
190
        {
 
191
            E00ReadRewind(poGDS->e00ReadPtr);
 
192
            for(i=0;i<6;i++)
 
193
                E00ReadNextLine(poGDS->e00ReadPtr);
 
194
        }
 
195
 
 
196
        if (nBlockYOff == poGDS->nLastYOff + 1)
 
197
        {
 
198
        }
 
199
        else if (nBlockYOff <= poGDS->nMaxYOffset)
 
200
        {
 
201
            //CPLDebug("E00GRID", "Skip to %d from %d", nBlockYOff, poGDS->nLastYOff);
 
202
            VSIFSeekL(poGDS->fp, poGDS->panOffsets[nBlockYOff], SEEK_SET);
 
203
            poGDS->nPosBeforeReadLine = poGDS->panOffsets[nBlockYOff];
 
204
            poGDS->e00ReadPtr->iInBufPtr = 0;
 
205
            poGDS->e00ReadPtr->szInBuf[0] = '\0';
 
206
        }
 
207
        else if (nBlockYOff > poGDS->nLastYOff + 1)
 
208
        {
 
209
            //CPLDebug("E00GRID", "Forward skip to %d from %d", nBlockYOff, poGDS->nLastYOff);
 
210
            for(i=poGDS->nLastYOff + 1; i < nBlockYOff;i++)
 
211
                IReadBlock(0, i, pImage);
 
212
        }
 
213
 
 
214
        if (nBlockYOff > poGDS->nMaxYOffset)
 
215
        {
 
216
            poGDS->panOffsets[nBlockYOff] = poGDS->nPosBeforeReadLine +
 
217
                                            poGDS->e00ReadPtr->iInBufPtr;
 
218
            poGDS->nMaxYOffset = nBlockYOff;
 
219
        }
 
220
 
 
221
        const char* pszLine = NULL;
 
222
        for(i=0;i<nBlockXSize;i++)
 
223
        {
 
224
            if ((i % VALS_PER_LINE) == 0)
 
225
            {
 
226
                pszLine = E00ReadNextLine(poGDS->e00ReadPtr);
 
227
                if (pszLine == NULL || strlen(pszLine) < 5 * E00_FLOAT_SIZE)
 
228
                    return CE_Failure;
 
229
            }
 
230
            if (eDataType == GDT_Float32)
 
231
            {
 
232
                pafImage[i] = (float) atof(pszLine + (i%VALS_PER_LINE) * E00_FLOAT_SIZE);
 
233
                /* Workaround single vs double precision problems */
 
234
                if (fNoData != 0 && fabs((pafImage[i] - fNoData)/fNoData) < 1e-6)
 
235
                    pafImage[i] = fNoData;
 
236
            }
 
237
            else
 
238
            {
 
239
                panImage[i] = atoi(pszLine + (i%VALS_PER_LINE) * E00_FLOAT_SIZE);
 
240
            }
 
241
        }
 
242
 
 
243
        poGDS->nLastYOff = nBlockYOff;
 
244
 
 
245
        return CE_None;
 
246
    }
 
247
 
 
248
    vsi_l_offset nValsToSkip = (vsi_l_offset)nBlockYOff * nRoundedBlockXSize;
 
249
    vsi_l_offset nLinesToSkip = nValsToSkip / VALS_PER_LINE;
 
250
    int nBytesPerLine = VALS_PER_LINE * E00_FLOAT_SIZE + poGDS->nBytesEOL;
 
251
    vsi_l_offset nPos = poGDS->nDataStart + nLinesToSkip * nBytesPerLine;
 
252
    VSIFSeekL(poGDS->fp, nPos, SEEK_SET);
 
253
 
 
254
    for(i=0;i<nBlockXSize;i++)
 
255
    {
 
256
        if (VSIFReadL(szVal, E00_FLOAT_SIZE, 1, poGDS->fp) != 1)
 
257
            return CE_Failure;
 
258
 
 
259
        if (eDataType == GDT_Float32)
 
260
        {
 
261
            pafImage[i] = (float) atof(szVal);
 
262
            /* Workaround single vs double precision problems */
 
263
            if (fNoData != 0 && fabs((pafImage[i] - fNoData)/fNoData) < 1e-6)
 
264
                pafImage[i] = fNoData;
 
265
        }
 
266
        else
 
267
        {
 
268
            panImage[i] = atoi(szVal);
 
269
        }
 
270
 
 
271
        if (((i+1) % VALS_PER_LINE) == 0)
 
272
            VSIFReadL(szVal, poGDS->nBytesEOL, 1, poGDS->fp);
 
273
    }
 
274
 
 
275
    return CE_None;
 
276
}
 
277
 
 
278
/************************************************************************/
 
279
/*                           GetNoDataValue()                           */
 
280
/************************************************************************/
 
281
 
 
282
double E00GRIDRasterBand::GetNoDataValue( int *pbSuccess )
 
283
{
 
284
    E00GRIDDataset *poGDS = (E00GRIDDataset *) poDS;
 
285
 
 
286
    if (pbSuccess)
 
287
        *pbSuccess = TRUE;
 
288
 
 
289
    if (eDataType == GDT_Float32)
 
290
        return (double)(float) poGDS->dfNoData;
 
291
    else
 
292
        return (double)(int)poGDS->dfNoData;
 
293
}
 
294
 
 
295
/************************************************************************/
 
296
/*                             GetUnitType()                            */
 
297
/************************************************************************/
 
298
 
 
299
const char * E00GRIDRasterBand::GetUnitType()
 
300
{
 
301
    E00GRIDDataset *poGDS = (E00GRIDDataset *) poDS;
 
302
 
 
303
    poGDS->ReadMetadata();
 
304
 
 
305
    if (poGDS->papszPrj == NULL)
 
306
        return GDALPamRasterBand::GetUnitType();
 
307
 
 
308
    char** papszIter = poGDS->papszPrj;
 
309
    const char* pszRet = "";
 
310
    while(*papszIter)
 
311
    {
 
312
        if (EQUALN(*papszIter, "Zunits", 6))
 
313
        {
 
314
            char** papszTokens = CSLTokenizeString(*papszIter);
 
315
            if (CSLCount(papszTokens) == 2)
 
316
            {
 
317
                if (EQUAL(papszTokens[1], "FEET"))
 
318
                    pszRet = "ft";
 
319
                else if (EQUAL(papszTokens[1], "METERS"))
 
320
                    pszRet = "m";
 
321
            }
 
322
            CSLDestroy(papszTokens);
 
323
            break;
 
324
        }
 
325
        papszIter ++;
 
326
    }
 
327
 
 
328
    return pszRet;
 
329
}
 
330
 
 
331
/************************************************************************/
 
332
/*                           GetMinimum()                               */
 
333
/************************************************************************/
 
334
 
 
335
double E00GRIDRasterBand::GetMinimum( int *pbSuccess )
 
336
{
 
337
    E00GRIDDataset *poGDS = (E00GRIDDataset *) poDS;
 
338
 
 
339
    poGDS->ReadMetadata();
 
340
 
 
341
    if (poGDS->bHasStats)
 
342
    {
 
343
        if( pbSuccess != NULL )
 
344
            *pbSuccess = TRUE;
 
345
 
 
346
        return poGDS->dfMin;
 
347
    }
 
348
 
 
349
    return GDALPamRasterBand::GetMinimum( pbSuccess );
 
350
}
 
351
 
 
352
/************************************************************************/
 
353
/*                           GetMaximum()                               */
 
354
/************************************************************************/
 
355
 
 
356
double E00GRIDRasterBand::GetMaximum( int *pbSuccess )
 
357
{
 
358
    E00GRIDDataset *poGDS = (E00GRIDDataset *) poDS;
 
359
 
 
360
    poGDS->ReadMetadata();
 
361
 
 
362
    if (poGDS->bHasStats)
 
363
    {
 
364
        if( pbSuccess != NULL )
 
365
            *pbSuccess = TRUE;
 
366
 
 
367
        return poGDS->dfMax;
 
368
    }
 
369
 
 
370
    return GDALPamRasterBand::GetMaximum( pbSuccess );
 
371
}
 
372
 
 
373
/************************************************************************/
 
374
/*                            GetStatistics()                           */
 
375
/************************************************************************/
 
376
 
 
377
CPLErr E00GRIDRasterBand::GetStatistics( int bApproxOK, int bForce,
 
378
                                         double *pdfMin, double *pdfMax,
 
379
                                         double *pdfMean, double *pdfStdDev )
 
380
{
 
381
    E00GRIDDataset *poGDS = (E00GRIDDataset *) poDS;
 
382
 
 
383
    poGDS->ReadMetadata();
 
384
 
 
385
    if (poGDS->bHasStats)
 
386
    {
 
387
        if (pdfMin)
 
388
            *pdfMin = poGDS->dfMin;
 
389
        if (pdfMax)
 
390
            *pdfMax = poGDS->dfMax;
 
391
        if (pdfMean)
 
392
            *pdfMean = poGDS->dfMean;
 
393
        if (pdfStdDev)
 
394
            *pdfStdDev = poGDS->dfStddev;
 
395
        return CE_None;
 
396
    }
 
397
 
 
398
    return GDALPamRasterBand::GetStatistics(bApproxOK, bForce,
 
399
                                            pdfMin, pdfMax,
 
400
                                            pdfMean, pdfStdDev);
 
401
}
 
402
 
 
403
/************************************************************************/
 
404
/*                           E00GRIDDataset()                           */
 
405
/************************************************************************/
 
406
 
 
407
E00GRIDDataset::E00GRIDDataset()
 
408
{
 
409
    e00ReadPtr = NULL;
 
410
    fp = NULL;
 
411
    nDataStart = 0;
 
412
    nBytesEOL = 1;
 
413
 
 
414
    nPosBeforeReadLine = 0;
 
415
    panOffsets = NULL;
 
416
    nLastYOff = -1;
 
417
    nMaxYOffset = -1;
 
418
 
 
419
    adfGeoTransform[0] = 0;
 
420
    adfGeoTransform[1] = 1;
 
421
    adfGeoTransform[2] = 0;
 
422
    adfGeoTransform[3] = 0;
 
423
    adfGeoTransform[4] = 0;
 
424
    adfGeoTransform[5] = 1;
 
425
 
 
426
    dfNoData = 0;
 
427
 
 
428
    papszPrj = NULL;
 
429
 
 
430
    bHasReadMetadata = FALSE;
 
431
 
 
432
    bHasStats = FALSE;
 
433
    dfMin = 0;
 
434
    dfMax = 0;
 
435
    dfMean = 0;
 
436
    dfStddev = 0;
 
437
}
 
438
 
 
439
/************************************************************************/
 
440
/*                           ~E00GRIDDataset()                          */
 
441
/************************************************************************/
 
442
 
 
443
E00GRIDDataset::~E00GRIDDataset()
 
444
 
 
445
{
 
446
    FlushCache();
 
447
    if (fp)
 
448
        VSIFCloseL(fp);
 
449
    CSLDestroy(papszPrj);
 
450
    E00ReadClose(e00ReadPtr);
 
451
    CPLFree(panOffsets);
 
452
}
 
453
 
 
454
/************************************************************************/
 
455
/*                             Identify()                               */
 
456
/************************************************************************/
 
457
 
 
458
int E00GRIDDataset::Identify( GDALOpenInfo * poOpenInfo )
 
459
{
 
460
    if (poOpenInfo->nHeaderBytes == 0)
 
461
        return FALSE;
 
462
 
 
463
    if (!(EQUALN((const char*)poOpenInfo->pabyHeader, "EXP  0", 6) ||
 
464
          EQUALN((const char*)poOpenInfo->pabyHeader, "EXP  1", 6)))
 
465
        return FALSE;
 
466
 
 
467
    /* FIXME: handle GRD  3 if that ever exists ? */
 
468
    if (strstr((const char*)poOpenInfo->pabyHeader, "GRD  2") == NULL)
 
469
        return FALSE;
 
470
 
 
471
    return TRUE;
 
472
}
 
473
 
 
474
/************************************************************************/
 
475
/*                            ReadNextLine()                            */
 
476
/************************************************************************/
 
477
 
 
478
const char* E00GRIDDataset::ReadNextLine(void * ptr)
 
479
{
 
480
    E00GRIDDataset* poDS = (E00GRIDDataset*) ptr;
 
481
    poDS->nPosBeforeReadLine = VSIFTellL(poDS->fp);
 
482
    return CPLReadLine2L(poDS->fp, 256, NULL);
 
483
}
 
484
 
 
485
/************************************************************************/
 
486
/*                                Rewind()                              */
 
487
/************************************************************************/
 
488
 
 
489
void E00GRIDDataset::Rewind(void * ptr)
 
490
{
 
491
    E00GRIDDataset* poDS = (E00GRIDDataset*) ptr;
 
492
    VSIRewindL(poDS->fp);
 
493
}
 
494
 
 
495
/************************************************************************/
 
496
/*                                Open()                                */
 
497
/************************************************************************/
 
498
 
 
499
GDALDataset *E00GRIDDataset::Open( GDALOpenInfo * poOpenInfo )
 
500
 
 
501
{
 
502
    int         i;
 
503
    GDALDataType eDT = GDT_Float32;
 
504
 
 
505
    if (!Identify(poOpenInfo))
 
506
        return NULL;
 
507
 
 
508
/* -------------------------------------------------------------------- */
 
509
/*      Find dataset characteristics                                    */
 
510
/* -------------------------------------------------------------------- */
 
511
    VSILFILE* fp = VSIFOpenL(poOpenInfo->pszFilename, "rb");
 
512
    if (fp == NULL)
 
513
        return NULL;
 
514
 
 
515
    if (poOpenInfo->eAccess == GA_Update)
 
516
    {
 
517
        CPLError( CE_Failure, CPLE_NotSupported, 
 
518
                  "The E00GRID driver does not support update access to existing"
 
519
                  " datasets.\n" );
 
520
        VSIFCloseL(fp);
 
521
        return NULL;
 
522
    }
 
523
 
 
524
/* -------------------------------------------------------------------- */
 
525
/*      Create a corresponding GDALDataset.                             */
 
526
/* -------------------------------------------------------------------- */
 
527
    E00GRIDDataset         *poDS;
 
528
 
 
529
    poDS = new E00GRIDDataset();
 
530
    if (strstr((const char*)poOpenInfo->pabyHeader, "\r\n") != NULL)
 
531
        poDS->nBytesEOL = 2;
 
532
    poDS->fp = fp;
 
533
 
 
534
    const char* pszLine;
 
535
    /* read EXP  0 or EXP  1 line */
 
536
    pszLine = CPLReadLine2L(fp, 81, NULL);
 
537
    if (pszLine == NULL)
 
538
    {
 
539
        CPLDebug("E00GRID", "Bad 1st line");
 
540
        delete poDS;
 
541
        return NULL;
 
542
    }
 
543
    int bCompressed = EQUALN(pszLine, "EXP  1", 6);
 
544
 
 
545
    E00ReadPtr e00ReadPtr = NULL;
 
546
    if (bCompressed)
 
547
    {
 
548
        VSIRewindL(fp);
 
549
        e00ReadPtr = E00ReadCallbackOpen(poDS,
 
550
                                         E00GRIDDataset::ReadNextLine,
 
551
                                         E00GRIDDataset::Rewind);
 
552
        if (e00ReadPtr == NULL)
 
553
        {
 
554
            delete poDS;
 
555
            return NULL;
 
556
        }
 
557
        E00ReadNextLine(e00ReadPtr);
 
558
        poDS->e00ReadPtr = e00ReadPtr;
 
559
    }
 
560
 
 
561
    /* skip GRD  2 line */
 
562
    if (e00ReadPtr)
 
563
        pszLine = E00ReadNextLine(e00ReadPtr);
 
564
    else
 
565
        pszLine = CPLReadLine2L(fp, 81, NULL);
 
566
    if (pszLine == NULL || !EQUALN(pszLine, "GRD  2", 6))
 
567
    {
 
568
        CPLDebug("E00GRID", "Bad 2nd line");
 
569
        delete poDS;
 
570
        return NULL;
 
571
    }
 
572
 
 
573
    /* read ncols, nrows and nodata value */
 
574
    if (e00ReadPtr)
 
575
        pszLine = E00ReadNextLine(e00ReadPtr);
 
576
    else
 
577
        pszLine = CPLReadLine2L(fp, 81, NULL);
 
578
    if (pszLine == NULL || strlen(pszLine) <
 
579
                E00_INT_SIZE+E00_INT_SIZE+2+E00_DOUBLE_SIZE)
 
580
    {
 
581
        CPLDebug("E00GRID", "Bad 3rd line");
 
582
        delete poDS;
 
583
        return NULL;
 
584
    }
 
585
 
 
586
    int nRasterXSize = atoi(pszLine);
 
587
    int nRasterYSize = atoi(pszLine + E00_INT_SIZE);
 
588
 
 
589
    if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
 
590
    {
 
591
        delete poDS;
 
592
        return NULL;
 
593
    }
 
594
 
 
595
    if (EQUALN(pszLine + E00_INT_SIZE + E00_INT_SIZE, " 1", 2))
 
596
        eDT = GDT_Int32;
 
597
    else if (EQUALN(pszLine + E00_INT_SIZE + E00_INT_SIZE, " 2", 2))
 
598
        eDT = GDT_Float32;
 
599
    else
 
600
    {
 
601
        CPLDebug("E00GRID", "Unknown data type : %s", pszLine);
 
602
    }
 
603
 
 
604
    double dfNoData = atof(pszLine + E00_INT_SIZE + E00_INT_SIZE + 2);
 
605
 
 
606
    /* read pixel size */
 
607
    if (e00ReadPtr)
 
608
        pszLine = E00ReadNextLine(e00ReadPtr);
 
609
    else
 
610
        pszLine = CPLReadLine2L(fp, 81, NULL);
 
611
    if (pszLine == NULL || strlen(pszLine) < 2*E00_DOUBLE_SIZE)
 
612
    {
 
613
        CPLDebug("E00GRID", "Bad 4th line");
 
614
        delete poDS;
 
615
        return NULL;
 
616
    }
 
617
/*
 
618
    double dfPixelX = atof(pszLine);
 
619
    double dfPixelY = atof(pszLine + E00_DOUBLE_SIZE);
 
620
*/
 
621
 
 
622
    /* read xmin, ymin */
 
623
    if (e00ReadPtr)
 
624
        pszLine = E00ReadNextLine(e00ReadPtr);
 
625
    else
 
626
        pszLine = CPLReadLine2L(fp, 81, NULL);
 
627
    if (pszLine == NULL || strlen(pszLine) < 2*E00_DOUBLE_SIZE)
 
628
    {
 
629
        CPLDebug("E00GRID", "Bad 5th line");
 
630
        delete poDS;
 
631
        return NULL;
 
632
    }
 
633
    double dfMinX = atof(pszLine);
 
634
    double dfMinY = atof(pszLine + E00_DOUBLE_SIZE);
 
635
 
 
636
    /* read xmax, ymax */
 
637
    if (e00ReadPtr)
 
638
        pszLine = E00ReadNextLine(e00ReadPtr);
 
639
    else
 
640
        pszLine = CPLReadLine2L(fp, 81, NULL);
 
641
    if (pszLine == NULL || strlen(pszLine) < 2*E00_DOUBLE_SIZE)
 
642
    {
 
643
        CPLDebug("E00GRID", "Bad 6th line");
 
644
        delete poDS;
 
645
        return NULL;
 
646
    }
 
647
    double dfMaxX = atof(pszLine);
 
648
    double dfMaxY = atof(pszLine + E00_DOUBLE_SIZE);
 
649
 
 
650
    poDS->nRasterXSize = nRasterXSize;
 
651
    poDS->nRasterYSize = nRasterYSize;
 
652
    poDS->dfNoData = dfNoData;
 
653
    poDS->adfGeoTransform[0] = dfMinX;
 
654
    poDS->adfGeoTransform[1] = (dfMaxX - dfMinX) / nRasterXSize;
 
655
    poDS->adfGeoTransform[2] = 0;
 
656
    poDS->adfGeoTransform[3] = dfMaxY;
 
657
    poDS->adfGeoTransform[4] = 0;
 
658
    poDS->adfGeoTransform[5] = - (dfMaxY - dfMinY) / nRasterYSize;
 
659
    poDS->nDataStart = VSIFTellL(fp);
 
660
    if (bCompressed)
 
661
    {
 
662
        poDS->panOffsets = (vsi_l_offset*)
 
663
                        VSIMalloc2(sizeof(vsi_l_offset), nRasterYSize);
 
664
        if (poDS->panOffsets == NULL)
 
665
        {
 
666
            delete poDS;
 
667
            return NULL;
 
668
        }
 
669
    }
 
670
/* -------------------------------------------------------------------- */
 
671
/*      Create band information objects.                                */
 
672
/* -------------------------------------------------------------------- */
 
673
    poDS->nBands = 1;
 
674
    for( i = 0; i < poDS->nBands; i++ )
 
675
        poDS->SetBand( i+1, new E00GRIDRasterBand( poDS, i+1, eDT ) );
 
676
 
 
677
/* -------------------------------------------------------------------- */
 
678
/*      Initialize any PAM information.                                 */
 
679
/* -------------------------------------------------------------------- */
 
680
    poDS->SetDescription( poOpenInfo->pszFilename );
 
681
    poDS->TryLoadXML();
 
682
 
 
683
/* -------------------------------------------------------------------- */
 
684
/*      Support overviews.                                              */
 
685
/* -------------------------------------------------------------------- */
 
686
    poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
 
687
    return( poDS );
 
688
}
 
689
 
 
690
/************************************************************************/
 
691
/*                          GetGeoTransform()                           */
 
692
/************************************************************************/
 
693
 
 
694
CPLErr E00GRIDDataset::GetGeoTransform( double * padfTransform )
 
695
 
 
696
{
 
697
    memcpy(padfTransform, adfGeoTransform, 6 * sizeof(double));
 
698
 
 
699
    return( CE_None );
 
700
}
 
701
 
 
702
 
 
703
/************************************************************************/
 
704
/*                             ReadLine()                               */
 
705
/************************************************************************/
 
706
 
 
707
const char* E00GRIDDataset::ReadLine()
 
708
{
 
709
    if (e00ReadPtr)
 
710
        return E00ReadNextLine(e00ReadPtr);
 
711
    else
 
712
        return CPLReadLine2L(fp, 81, NULL);
 
713
}
 
714
 
 
715
/************************************************************************/
 
716
/*                         GetProjectionRef()                           */
 
717
/************************************************************************/
 
718
 
 
719
const char* E00GRIDDataset::GetProjectionRef()
 
720
 
 
721
{
 
722
    ReadMetadata();
 
723
    return osProjection.c_str();
 
724
}
 
725
 
 
726
/************************************************************************/
 
727
/*                          ReadMetadata()                              */
 
728
/************************************************************************/
 
729
 
 
730
void E00GRIDDataset::ReadMetadata()
 
731
 
 
732
{
 
733
    if (bHasReadMetadata)
 
734
        return;
 
735
 
 
736
    bHasReadMetadata = TRUE;
 
737
 
 
738
    if (e00ReadPtr == NULL)
 
739
    {
 
740
        int nRoundedBlockXSize = ((nRasterXSize + VALS_PER_LINE - 1) /
 
741
                                                VALS_PER_LINE) * VALS_PER_LINE;
 
742
        vsi_l_offset nValsToSkip =
 
743
                               (vsi_l_offset)nRasterYSize * nRoundedBlockXSize;
 
744
        vsi_l_offset nLinesToSkip = nValsToSkip / VALS_PER_LINE;
 
745
        int nBytesPerLine = VALS_PER_LINE * E00_FLOAT_SIZE + nBytesEOL;
 
746
        vsi_l_offset nPos = nDataStart + nLinesToSkip * nBytesPerLine;
 
747
        VSIFSeekL(fp, nPos, SEEK_SET);
 
748
    }
 
749
    else
 
750
    {
 
751
        nLastYOff = -1;
 
752
 
 
753
        const unsigned int BUFFER_SIZE = 65536;
 
754
        const unsigned int NEEDLE_SIZE = 3*5;
 
755
        const unsigned int nToRead = BUFFER_SIZE - NEEDLE_SIZE;
 
756
        char* pabyBuffer = (char*)CPLCalloc(1, BUFFER_SIZE+NEEDLE_SIZE);
 
757
        int nRead;
 
758
        int bEOGFound = FALSE;
 
759
 
 
760
        VSIFSeekL(fp, 0, SEEK_END);
 
761
        vsi_l_offset nEndPos = VSIFTellL(fp);
 
762
        if (nEndPos > BUFFER_SIZE)
 
763
            nEndPos -= BUFFER_SIZE;
 
764
        else
 
765
            nEndPos = 0;
 
766
        VSIFSeekL(fp, nEndPos, SEEK_SET);
 
767
 
 
768
#define GOTO_NEXT_CHAR() \
 
769
    i ++; \
 
770
    if (pabyBuffer[i] == 13 || pabyBuffer[i] == 10) \
 
771
    { \
 
772
        i++; \
 
773
        if (pabyBuffer[i] == 10) \
 
774
            i++; \
 
775
    } \
 
776
 
 
777
        while ((nRead = VSIFReadL(pabyBuffer, 1, nToRead, fp)) != 0)
 
778
        {
 
779
            int i;
 
780
            for(i = 0; i < nRead; i++)
 
781
            {
 
782
                if (pabyBuffer[i] == 'E')
 
783
                {
 
784
                    GOTO_NEXT_CHAR();
 
785
                    if (pabyBuffer[i] == 'O')
 
786
                    {
 
787
                        GOTO_NEXT_CHAR();
 
788
                        if (pabyBuffer[i] == 'G')
 
789
                        {
 
790
                            GOTO_NEXT_CHAR();
 
791
                            if (pabyBuffer[i] == '~')
 
792
                            {
 
793
                                GOTO_NEXT_CHAR();
 
794
                                if (pabyBuffer[i] == '}')
 
795
                                {
 
796
                                    bEOGFound = TRUE;
 
797
                                    break;
 
798
                                }
 
799
                            }
 
800
                        }
 
801
                    }
 
802
                }
 
803
            }
 
804
 
 
805
            if (bEOGFound)
 
806
            {
 
807
                VSIFSeekL(fp, VSIFTellL(fp) - nRead + i + 1, SEEK_SET);
 
808
                e00ReadPtr->iInBufPtr = 0;
 
809
                e00ReadPtr->szInBuf[0] = '\0';
 
810
                break;
 
811
            }
 
812
 
 
813
            if (nEndPos == 0)
 
814
                break;
 
815
 
 
816
            if ((unsigned int)nRead == nToRead)
 
817
            {
 
818
                memmove(pabyBuffer + nToRead, pabyBuffer, NEEDLE_SIZE);
 
819
                if (nEndPos >= (vsi_l_offset)nToRead)
 
820
                    nEndPos -= nToRead;
 
821
                else
 
822
                    nEndPos = 0;
 
823
                VSIFSeekL(fp, nEndPos, SEEK_SET);
 
824
            }
 
825
            else
 
826
                break;
 
827
        }
 
828
        CPLFree(pabyBuffer);
 
829
        if (!bEOGFound)
 
830
            return;
 
831
    }
 
832
 
 
833
    const char* pszLine;
 
834
    int bPRJFound = FALSE;
 
835
    int bStatsFound = FALSE;
 
836
    while((pszLine = ReadLine()) != NULL)
 
837
    {
 
838
        if (EQUALN(pszLine, "PRJ  2", 6))
 
839
        {
 
840
            bPRJFound = TRUE;
 
841
            while((pszLine = ReadLine()) != NULL)
 
842
            {
 
843
                if (EQUAL(pszLine, "EOP"))
 
844
                {
 
845
                    break;
 
846
                }
 
847
                papszPrj = CSLAddString(papszPrj, pszLine);
 
848
            }
 
849
 
 
850
            OGRSpatialReference oSRS;
 
851
            if( oSRS.importFromESRI( papszPrj ) != OGRERR_NONE )
 
852
            {
 
853
                CPLError( CE_Warning, CPLE_AppDefined,
 
854
                            "Failed to parse PRJ section, ignoring." );
 
855
            }
 
856
            else
 
857
            {
 
858
                char* pszWKT = NULL;
 
859
                if (oSRS.exportToWkt(&pszWKT) == OGRERR_NONE && pszWKT != NULL)
 
860
                    osProjection = pszWKT;
 
861
                CPLFree(pszWKT);
 
862
            }
 
863
            if (bStatsFound)
 
864
                break;
 
865
        }
 
866
        else if (strcmp(pszLine, "STDV              8-1  254-1  15 3 60-1  -1  -1-1                   4-") == 0)
 
867
        {
 
868
            bStatsFound = TRUE;
 
869
            pszLine = ReadLine();
 
870
            if (pszLine)
 
871
            {
 
872
                CPLString osStats = pszLine;
 
873
                pszLine = ReadLine();
 
874
                if (pszLine)
 
875
                {
 
876
                    osStats += pszLine;
 
877
                    char** papszTokens = CSLTokenizeString(osStats);
 
878
                    if (CSLCount(papszTokens) == 4)
 
879
                    {
 
880
                        dfMin = atof(papszTokens[0]);
 
881
                        dfMax = atof(papszTokens[1]);
 
882
                        dfMean = atof(papszTokens[2]);
 
883
                        dfStddev = atof(papszTokens[3]);
 
884
                        bHasStats = TRUE;
 
885
                    }
 
886
                    CSLDestroy(papszTokens);
 
887
                }
 
888
            }
 
889
            if (bPRJFound)
 
890
                break;
 
891
        }
 
892
    }
 
893
}
 
894
 
 
895
/************************************************************************/
 
896
/*                       GDALRegister_E00GRID()                         */
 
897
/************************************************************************/
 
898
 
 
899
void GDALRegister_E00GRID()
 
900
 
 
901
{
 
902
    GDALDriver  *poDriver;
 
903
 
 
904
    if( GDALGetDriverByName( "E00GRID" ) == NULL )
 
905
    {
 
906
        poDriver = new GDALDriver();
 
907
        
 
908
        poDriver->SetDescription( "E00GRID" );
 
909
        poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, 
 
910
                                   "Arc/Info Export E00 GRID" );
 
911
        poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, 
 
912
                                   "frmt_various.html#E00GRID" );
 
913
        poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "e00" );
 
914
 
 
915
 
 
916
        poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
 
917
 
 
918
        poDriver->pfnOpen = E00GRIDDataset::Open;
 
919
        poDriver->pfnIdentify = E00GRIDDataset::Identify;
 
920
 
 
921
        GetGDALDriverManager()->RegisterDriver( poDriver );
 
922
    }
 
923
}
 
924