← Back to branch summary

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

« back to all changes in this revision

Viewing changes to frmts/hf2/hf2dataset.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/hf2/hf2dataset.cpp
 
1
/******************************************************************************
 
2
 * $Id: hf2dataset.cpp 21680 2011-02-11 21:12:07Z warmerdam $
 
3
 *
 
4
 * Project:  HF2 driver
 
5
 * Purpose:  GDALDataset driver for HF2/HFZ dataset.
 
6
 * Author:   Even Rouault, <even dot rouault at mines dash paris dot org>
 
7
 *
 
8
 ******************************************************************************
 
9
 * Copyright (c) 2010, 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_string.h"
 
31
#include "gdal_pam.h"
 
32
#include "ogr_spatialref.h"
 
33
 
 
34
CPL_CVSID("$Id: hf2dataset.cpp 21680 2011-02-11 21:12:07Z warmerdam $");
 
35
 
 
36
CPL_C_START
 
37
void    GDALRegister_HF2(void);
 
38
CPL_C_END
 
39
 
 
40
/************************************************************************/
 
41
/* ==================================================================== */
 
42
/*                              HF2Dataset                              */
 
43
/* ==================================================================== */
 
44
/************************************************************************/
 
45
 
 
46
class HF2RasterBand;
 
47
 
 
48
class HF2Dataset : public GDALPamDataset
 
49
{
 
50
    friend class HF2RasterBand;
 
51
    
 
52
    VSILFILE   *fp;
 
53
    double      adfGeoTransform[6];
 
54
    char       *pszWKT;
 
55
    vsi_l_offset    *panBlockOffset;
 
56
 
 
57
    int         nTileSize;
 
58
    int         bHasLoaderBlockMap;
 
59
    int         LoadBlockMap();
 
60
 
 
61
  public:
 
62
                 HF2Dataset();
 
63
    virtual     ~HF2Dataset();
 
64
    
 
65
    virtual CPLErr GetGeoTransform( double * );
 
66
    virtual const char* GetProjectionRef();
 
67
    
 
68
    static GDALDataset *Open( GDALOpenInfo * );
 
69
    static int          Identify( GDALOpenInfo * );
 
70
    static GDALDataset *CreateCopy( const char * pszFilename, GDALDataset *poSrcDS, 
 
71
                                    int bStrict, char ** papszOptions, 
 
72
                                    GDALProgressFunc pfnProgress, void * pProgressData );
 
73
};
 
74
 
 
75
/************************************************************************/
 
76
/* ==================================================================== */
 
77
/*                            HF2RasterBand                             */
 
78
/* ==================================================================== */
 
79
/************************************************************************/
 
80
 
 
81
class HF2RasterBand : public GDALPamRasterBand
 
82
{
 
83
    friend class HF2Dataset;
 
84
 
 
85
    float*  pafBlockData;
 
86
    int     nLastBlockYOff;
 
87
 
 
88
  public:
 
89
 
 
90
                HF2RasterBand( HF2Dataset *, int, GDALDataType );
 
91
               ~HF2RasterBand();
 
92
 
 
93
    virtual CPLErr IReadBlock( int, int, void * );
 
94
};
 
95
 
 
96
 
 
97
/************************************************************************/
 
98
/*                           HF2RasterBand()                            */
 
99
/************************************************************************/
 
100
 
 
101
HF2RasterBand::HF2RasterBand( HF2Dataset *poDS, int nBand, GDALDataType eDT )
 
102
 
 
103
{
 
104
    this->poDS = poDS;
 
105
    this->nBand = nBand;
 
106
 
 
107
    eDataType = eDT;
 
108
 
 
109
    nBlockXSize = poDS->nTileSize;
 
110
    nBlockYSize = 1;
 
111
 
 
112
    pafBlockData = NULL;
 
113
    nLastBlockYOff = -1;
 
114
}
 
115
 
 
116
/************************************************************************/
 
117
/*                          ~HF2RasterBand()                            */
 
118
/************************************************************************/
 
119
 
 
120
HF2RasterBand::~HF2RasterBand()
 
121
{
 
122
    CPLFree(pafBlockData);
 
123
}
 
124
 
 
125
/************************************************************************/
 
126
/*                             IReadBlock()                             */
 
127
/************************************************************************/
 
128
 
 
129
CPLErr HF2RasterBand::IReadBlock( int nBlockXOff, int nLineYOff,
 
130
                                  void * pImage )
 
131
 
 
132
{
 
133
    HF2Dataset *poGDS = (HF2Dataset *) poDS;
 
134
 
 
135
    int nXBlocks = (nRasterXSize + nBlockXSize - 1) / nBlockXSize;
 
136
    int nYBlocks = (nRasterYSize + nBlockXSize - 1) / nBlockXSize;
 
137
    
 
138
    if (!poGDS->LoadBlockMap())
 
139
        return CE_Failure;
 
140
    
 
141
    if (pafBlockData == NULL)
 
142
    {
 
143
        pafBlockData = (float*)VSIMalloc3(nXBlocks * sizeof(float), poGDS->nTileSize, poGDS->nTileSize);
 
144
        if (pafBlockData == NULL)
 
145
            return CE_Failure;
 
146
    }
 
147
    
 
148
    nLineYOff = nRasterYSize - 1 - nLineYOff;
 
149
 
 
150
    int nBlockYOff = nLineYOff / nBlockXSize;
 
151
    int nYOffInTile = nLineYOff % nBlockXSize;
 
152
 
 
153
    if (nBlockYOff != nLastBlockYOff)
 
154
    {
 
155
        nLastBlockYOff = nBlockYOff;
 
156
 
 
157
        memset(pafBlockData, 0, nXBlocks * sizeof(float) * nBlockXSize * nBlockXSize);
 
158
 
 
159
        /* 4 * nBlockXSize is the upper bound */
 
160
        void* pabyData = CPLMalloc( 4 * nBlockXSize );
 
161
 
 
162
        int nxoff;
 
163
        for(nxoff = 0; nxoff < nXBlocks; nxoff++)
 
164
        {
 
165
            VSIFSeekL(poGDS->fp, poGDS->panBlockOffset[(nYBlocks - 1 - nBlockYOff) * nXBlocks + nxoff], SEEK_SET);
 
166
            float fScale, fOff;
 
167
            VSIFReadL(&fScale, 4, 1, poGDS->fp);
 
168
            VSIFReadL(&fOff, 4, 1, poGDS->fp);
 
169
            CPL_LSBPTR32(&fScale);
 
170
            CPL_LSBPTR32(&fOff);
 
171
    
 
172
            int nTileWidth = MIN(nBlockXSize, nRasterXSize - nxoff * nBlockXSize);
 
173
            int nTileHeight = MIN(nBlockXSize, nRasterYSize - nBlockYOff * nBlockXSize);
 
174
            
 
175
            int j;
 
176
            for(j=0;j<nTileHeight;j++)
 
177
            {
 
178
                GByte nWordSize;
 
179
                VSIFReadL(&nWordSize, 1, 1, poGDS->fp);
 
180
                if (nWordSize != 1 && nWordSize != 2 && nWordSize != 4)
 
181
                {
 
182
                    CPLError(CE_Failure, CPLE_AppDefined,
 
183
                             "Unexpected word size : %d", (int)nWordSize);
 
184
                    break;
 
185
                }
 
186
 
 
187
                GInt32 nVal;
 
188
                VSIFReadL(&nVal, 4, 1, poGDS->fp);
 
189
                CPL_LSBPTR32(&nVal);
 
190
                VSIFReadL(pabyData, nWordSize * (nTileWidth - 1), 1, poGDS->fp);
 
191
#if defined(CPL_MSB)
 
192
                if (nWordSize > 1)
 
193
                    GDALSwapWords(pabyData, nWordSize, nTileWidth - 1, nWordSize);
 
194
#endif
 
195
 
 
196
                pafBlockData[nxoff * nBlockXSize * nBlockXSize + j * nBlockXSize + 0] = nVal * fScale + fOff;
 
197
                int i;
 
198
                for(i=1;i<nTileWidth;i++)
 
199
                {
 
200
                    if (nWordSize == 1)
 
201
                        nVal += ((char*)pabyData)[i-1];
 
202
                    else if (nWordSize == 2)
 
203
                        nVal += ((GInt16*)pabyData)[i-1];
 
204
                    else
 
205
                        nVal += ((GInt32*)pabyData)[i-1];
 
206
                    pafBlockData[nxoff * nBlockXSize * nBlockXSize + j * nBlockXSize + i] = nVal * fScale + fOff;
 
207
                }
 
208
            }
 
209
        }
 
210
 
 
211
        CPLFree(pabyData);
 
212
    }
 
213
 
 
214
    int nTileWidth = MIN(nBlockXSize, nRasterXSize - nBlockXOff * nBlockXSize);
 
215
    memcpy(pImage, pafBlockData + nBlockXOff * nBlockXSize * nBlockXSize +
 
216
                                  nYOffInTile * nBlockXSize,
 
217
           nTileWidth * sizeof(float));
 
218
 
 
219
    return CE_None;
 
220
}
 
221
 
 
222
/************************************************************************/
 
223
/*                            ~HF2Dataset()                            */
 
224
/************************************************************************/
 
225
 
 
226
HF2Dataset::HF2Dataset()
 
227
{
 
228
    fp = NULL;
 
229
    pszWKT = NULL;
 
230
    panBlockOffset = NULL;
 
231
    adfGeoTransform[0] = 0;
 
232
    adfGeoTransform[1] = 1;
 
233
    adfGeoTransform[2] = 0;
 
234
    adfGeoTransform[3] = 0;
 
235
    adfGeoTransform[4] = 0;
 
236
    adfGeoTransform[5] = 1;
 
237
    bHasLoaderBlockMap = FALSE;
 
238
    nTileSize = 0;
 
239
}
 
240
 
 
241
/************************************************************************/
 
242
/*                            ~HF2Dataset()                            */
 
243
/************************************************************************/
 
244
 
 
245
HF2Dataset::~HF2Dataset()
 
246
 
 
247
{
 
248
    FlushCache();
 
249
    CPLFree(pszWKT);
 
250
    CPLFree(panBlockOffset);
 
251
    if (fp)
 
252
        VSIFCloseL(fp);
 
253
}
 
254
 
 
255
/************************************************************************/
 
256
/*                            LoadBlockMap()                            */
 
257
/************************************************************************/
 
258
 
 
259
int HF2Dataset::LoadBlockMap()
 
260
{
 
261
    if (bHasLoaderBlockMap)
 
262
        return panBlockOffset != NULL;
 
263
 
 
264
    bHasLoaderBlockMap = TRUE;
 
265
 
 
266
    int nXBlocks = (nRasterXSize + nTileSize - 1) / nTileSize;
 
267
    int nYBlocks = (nRasterYSize + nTileSize - 1) / nTileSize;
 
268
    panBlockOffset = (vsi_l_offset*) VSIMalloc3(sizeof(vsi_l_offset), nXBlocks, nYBlocks);
 
269
    if (panBlockOffset == NULL)
 
270
    {
 
271
        return FALSE;
 
272
    }
 
273
    int i, j;
 
274
    for(j = 0; j < nYBlocks; j++)
 
275
    {
 
276
        for(i = 0; i < nXBlocks; i++)
 
277
        {
 
278
            vsi_l_offset nOff = VSIFTellL(fp);
 
279
            panBlockOffset[(nYBlocks - 1 - j) * nXBlocks + i] = nOff;
 
280
            //VSIFSeekL(fp, 4 + 4, SEEK_CUR);
 
281
            float fScale, fOff;
 
282
            VSIFReadL(&fScale, 4, 1, fp);
 
283
            VSIFReadL(&fOff, 4, 1, fp);
 
284
            CPL_LSBPTR32(&fScale);
 
285
            CPL_LSBPTR32(&fOff);
 
286
            //printf("fScale = %f, fOff = %f\n", fScale, fOff);
 
287
            int k;
 
288
            int nCols = MIN(nTileSize, nRasterXSize - nTileSize *i);
 
289
            int nLines = MIN(nTileSize, nRasterYSize - nTileSize *j);
 
290
            for(k = 0; k < nLines; k++)
 
291
            {
 
292
                GByte nWordSize;
 
293
                VSIFReadL(&nWordSize, 1, 1, fp);
 
294
                //printf("nWordSize=%d\n", nWordSize);
 
295
                if (nWordSize == 1 || nWordSize == 2 || nWordSize == 4)
 
296
                    VSIFSeekL(fp, 4 + nWordSize * (nCols - 1), SEEK_CUR);
 
297
                else
 
298
                {
 
299
                    CPLError(CE_Failure, CPLE_AppDefined,
 
300
                            "Got unexpected byte depth (%d) for block (%d, %d) line %d",
 
301
                            (int)nWordSize, i, j, k);
 
302
                    VSIFree(panBlockOffset);
 
303
                    panBlockOffset = NULL;
 
304
                    return FALSE;
 
305
                }
 
306
            }
 
307
        }
 
308
    }
 
309
 
 
310
    return TRUE;
 
311
}
 
312
 
 
313
/************************************************************************/
 
314
/*                     GetProjectionRef()                               */
 
315
/************************************************************************/
 
316
 
 
317
const char* HF2Dataset::GetProjectionRef()
 
318
{
 
319
    if (pszWKT)
 
320
        return pszWKT;
 
321
    return GDALPamDataset::GetProjectionRef();
 
322
}
 
323
 
 
324
/************************************************************************/
 
325
/*                             Identify()                               */
 
326
/************************************************************************/
 
327
 
 
328
int HF2Dataset::Identify( GDALOpenInfo * poOpenInfo)
 
329
{
 
330
 
 
331
    GDALOpenInfo* poOpenInfoToDelete = NULL;
 
332
    /*  GZipped .hf2 files are common, so automagically open them */
 
333
    /*  if the /vsigzip/ has not been explicitely passed */
 
334
    CPLString osFilename(poOpenInfo->pszFilename);
 
335
    if ((EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "hfz") ||
 
336
        (strlen(poOpenInfo->pszFilename) > 6 &&
 
337
         EQUAL(poOpenInfo->pszFilename + strlen(poOpenInfo->pszFilename) - 6, "hf2.gz"))) &&
 
338
         !EQUALN(poOpenInfo->pszFilename, "/vsigzip/", 9))
 
339
    {
 
340
        osFilename = "/vsigzip/";
 
341
        osFilename += poOpenInfo->pszFilename;
 
342
        poOpenInfo = poOpenInfoToDelete =
 
343
                new GDALOpenInfo(osFilename.c_str(), GA_ReadOnly,
 
344
                                 poOpenInfo->papszSiblingFiles);
 
345
    }
 
346
 
 
347
    if (poOpenInfo->nHeaderBytes < 28)
 
348
    {
 
349
        delete poOpenInfoToDelete;
 
350
        return FALSE;
 
351
    }
 
352
 
 
353
    if (memcmp(poOpenInfo->pabyHeader, "HF2\0\0\0\0", 6) != 0)
 
354
    {
 
355
        delete poOpenInfoToDelete;
 
356
        return FALSE;
 
357
    }
 
358
 
 
359
    delete poOpenInfoToDelete;
 
360
    return TRUE;
 
361
}
 
362
 
 
363
/************************************************************************/
 
364
/*                                Open()                                */
 
365
/************************************************************************/
 
366
 
 
367
GDALDataset *HF2Dataset::Open( GDALOpenInfo * poOpenInfo )
 
368
 
 
369
{
 
370
    CPLString osOriginalFilename(poOpenInfo->pszFilename);
 
371
 
 
372
    if (!Identify(poOpenInfo))
 
373
        return NULL;
 
374
 
 
375
    GDALOpenInfo* poOpenInfoToDelete = NULL;
 
376
    /*  GZipped .hf2 files are common, so automagically open them */
 
377
    /*  if the /vsigzip/ has not been explicitely passed */
 
378
    CPLString osFilename(poOpenInfo->pszFilename);
 
379
    if ((EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "hfz") ||
 
380
        (strlen(poOpenInfo->pszFilename) > 6 &&
 
381
         EQUAL(poOpenInfo->pszFilename + strlen(poOpenInfo->pszFilename) - 6, "hf2.gz"))) &&
 
382
         !EQUALN(poOpenInfo->pszFilename, "/vsigzip/", 9))
 
383
    {
 
384
        osFilename = "/vsigzip/";
 
385
        osFilename += poOpenInfo->pszFilename;
 
386
        poOpenInfo = poOpenInfoToDelete =
 
387
                new GDALOpenInfo(osFilename.c_str(), GA_ReadOnly,
 
388
                                 poOpenInfo->papszSiblingFiles);
 
389
    }
 
390
 
 
391
/* -------------------------------------------------------------------- */
 
392
/*      Parse header                                                    */
 
393
/* -------------------------------------------------------------------- */
 
394
 
 
395
    int nXSize, nYSize;
 
396
    memcpy(&nXSize, poOpenInfo->pabyHeader + 6, 4);
 
397
    CPL_LSBPTR32(&nXSize);
 
398
    memcpy(&nYSize, poOpenInfo->pabyHeader + 10, 4);
 
399
    CPL_LSBPTR32(&nYSize);
 
400
 
 
401
    GUInt16 nTileSize;
 
402
    memcpy(&nTileSize, poOpenInfo->pabyHeader + 14, 2);
 
403
    CPL_LSBPTR16(&nTileSize);
 
404
 
 
405
    float fVertPres, fHorizScale;
 
406
    memcpy(&fVertPres, poOpenInfo->pabyHeader + 16, 4);
 
407
    CPL_LSBPTR32(&fVertPres);
 
408
    memcpy(&fHorizScale, poOpenInfo->pabyHeader + 20, 4);
 
409
    CPL_LSBPTR32(&fHorizScale);
 
410
 
 
411
    GUInt32 nExtendedHeaderLen;
 
412
    memcpy(&nExtendedHeaderLen, poOpenInfo->pabyHeader + 24, 4);
 
413
    CPL_LSBPTR32(&nExtendedHeaderLen);
 
414
 
 
415
    delete poOpenInfoToDelete;
 
416
    poOpenInfoToDelete = NULL;
 
417
 
 
418
    if (nTileSize < 8)
 
419
        return NULL;
 
420
    if (nXSize <= 0 || nXSize > INT_MAX - nTileSize ||
 
421
        nYSize <= 0 || nYSize > INT_MAX - nTileSize)
 
422
        return NULL;
 
423
    /* To avoid later potential int overflows */
 
424
    if (nExtendedHeaderLen > 1024 * 65536)
 
425
        return NULL;
 
426
 
 
427
    if (!GDALCheckDatasetDimensions(nXSize, nYSize))
 
428
    {
 
429
        return NULL;
 
430
    }
 
431
 
 
432
/* -------------------------------------------------------------------- */
 
433
/*      Parse extended blocks                                           */
 
434
/* -------------------------------------------------------------------- */
 
435
 
 
436
    VSILFILE* fp = VSIFOpenL(osFilename.c_str(), "rb");
 
437
    if (fp == NULL)
 
438
        return NULL;
 
439
 
 
440
    VSIFSeekL(fp, 28, SEEK_SET);
 
441
 
 
442
    int bHasExtent = FALSE;
 
443
    double dfMinX = 0, dfMaxX = 0, dfMinY = 0, dfMaxY = 0;
 
444
    int bHasUTMZone = FALSE;
 
445
    GInt16 nUTMZone = 0;
 
446
    int bHasEPSGDatumCode = FALSE;
 
447
    GInt16 nEPSGDatumCode = 0;
 
448
    int bHasEPSGCode = FALSE;
 
449
    GInt16 nEPSGCode = 0;
 
450
    int bHasRelativePrecision = FALSE;
 
451
    float fRelativePrecision = 0;
 
452
    char szApplicationName[256];
 
453
    szApplicationName[0] = 0;
 
454
 
 
455
    GUInt32 nExtendedHeaderOff = 0;
 
456
    while(nExtendedHeaderOff < nExtendedHeaderLen)
 
457
    {
 
458
        char pabyBlockHeader[24];
 
459
        VSIFReadL(pabyBlockHeader, 24, 1, fp);
 
460
 
 
461
        char szBlockName[16 + 1];
 
462
        memcpy(szBlockName, pabyBlockHeader + 4, 16);
 
463
        szBlockName[16] = 0;
 
464
        GUInt32 nBlockSize;
 
465
        memcpy(&nBlockSize, pabyBlockHeader + 20, 4);
 
466
        CPL_LSBPTR32(&nBlockSize);
 
467
        if (nBlockSize > 65536)
 
468
            break;
 
469
 
 
470
        nExtendedHeaderOff += 24 + nBlockSize;
 
471
 
 
472
        if (strcmp(szBlockName, "georef-extents") == 0 &&
 
473
            nBlockSize == 34)
 
474
        {
 
475
            char pabyBlockData[34];
 
476
            VSIFReadL(pabyBlockData, 34, 1, fp);
 
477
 
 
478
            memcpy(&dfMinX, pabyBlockData + 2, 8);
 
479
            CPL_LSBPTR64(&dfMinX);
 
480
            memcpy(&dfMaxX, pabyBlockData + 2 + 8, 8);
 
481
            CPL_LSBPTR64(&dfMaxX);
 
482
            memcpy(&dfMinY, pabyBlockData + 2 + 8 + 8, 8);
 
483
            CPL_LSBPTR64(&dfMinY);
 
484
            memcpy(&dfMaxY, pabyBlockData + 2 + 8 + 8 + 8, 8);
 
485
            CPL_LSBPTR64(&dfMaxY);
 
486
 
 
487
            bHasExtent = TRUE;
 
488
        }
 
489
        else if (strcmp(szBlockName, "georef-utm") == 0 &&
 
490
                nBlockSize == 2)
 
491
        {
 
492
            VSIFReadL(&nUTMZone, 2, 1, fp);
 
493
            CPL_LSBPTR16(&nUTMZone);
 
494
            CPLDebug("HF2", "UTM Zone = %d", nUTMZone);
 
495
 
 
496
            bHasUTMZone = TRUE;
 
497
        }
 
498
        else if (strcmp(szBlockName, "georef-datum") == 0 &&
 
499
                 nBlockSize == 2)
 
500
        {
 
501
            VSIFReadL(&nEPSGDatumCode, 2, 1, fp);
 
502
            CPL_LSBPTR16(&nEPSGDatumCode);
 
503
            CPLDebug("HF2", "EPSG Datum Code = %d", nEPSGDatumCode);
 
504
 
 
505
            bHasEPSGDatumCode = TRUE;
 
506
        }
 
507
        else if (strcmp(szBlockName, "georef-epsg-prj") == 0 &&
 
508
                 nBlockSize == 2)
 
509
        {
 
510
            VSIFReadL(&nEPSGCode, 2, 1, fp);
 
511
            CPL_LSBPTR16(&nEPSGCode);
 
512
            CPLDebug("HF2", "EPSG Code = %d", nEPSGCode);
 
513
 
 
514
            bHasEPSGCode = TRUE;
 
515
        }
 
516
        else if (strcmp(szBlockName, "precis-rel") == 0 &&
 
517
                 nBlockSize == 4)
 
518
        {
 
519
            VSIFReadL(&fRelativePrecision, 4, 1, fp);
 
520
            CPL_LSBPTR32(&fRelativePrecision);
 
521
 
 
522
            bHasRelativePrecision = TRUE;
 
523
        }
 
524
        else if (strcmp(szBlockName, "app-name") == 0 &&
 
525
                 nBlockSize < 256)
 
526
        {
 
527
            VSIFReadL(szApplicationName, nBlockSize, 1, fp);
 
528
            szApplicationName[nBlockSize] = 0;
 
529
        }
 
530
        else
 
531
        {
 
532
            CPLDebug("HF2", "Skipping block %s", szBlockName);
 
533
            VSIFSeekL(fp, nBlockSize, SEEK_CUR);
 
534
        }
 
535
    }
 
536
 
 
537
/* -------------------------------------------------------------------- */
 
538
/*      Create a corresponding GDALDataset.                             */
 
539
/* -------------------------------------------------------------------- */
 
540
    HF2Dataset         *poDS;
 
541
 
 
542
    poDS = new HF2Dataset();
 
543
    poDS->fp = fp;
 
544
    poDS->nRasterXSize = nXSize;
 
545
    poDS->nRasterYSize = nYSize;
 
546
    poDS->nTileSize = nTileSize;
 
547
    CPLDebug("HF2", "nXSize = %d, nYSize = %d, nTileSize = %d", nXSize, nYSize, nTileSize);
 
548
    if (bHasExtent)
 
549
    {
 
550
        poDS->adfGeoTransform[0] = dfMinX;
 
551
        poDS->adfGeoTransform[3] = dfMaxY;
 
552
        poDS->adfGeoTransform[1] = (dfMaxX - dfMinX) / nXSize;
 
553
        poDS->adfGeoTransform[5] = -(dfMaxY - dfMinY) / nYSize;
 
554
    }
 
555
    else
 
556
    {
 
557
        poDS->adfGeoTransform[1] = fHorizScale;
 
558
        poDS->adfGeoTransform[5] = fHorizScale;
 
559
    }
 
560
 
 
561
    if (bHasEPSGCode)
 
562
    {
 
563
        OGRSpatialReference oSRS;
 
564
        if (oSRS.importFromEPSG(nEPSGCode) == OGRERR_NONE)
 
565
            oSRS.exportToWkt(&poDS->pszWKT);
 
566
    }
 
567
    else
 
568
    {
 
569
        int bHasSRS = FALSE;
 
570
        OGRSpatialReference oSRS;
 
571
        oSRS.SetGeogCS("unknown", "unknown", "unknown", SRS_WGS84_SEMIMAJOR, SRS_WGS84_INVFLATTENING); 
 
572
        if (bHasEPSGDatumCode)
 
573
        {
 
574
            if (nEPSGDatumCode == 23 || nEPSGDatumCode == 6326)
 
575
            {
 
576
                bHasSRS = TRUE;
 
577
                oSRS.SetWellKnownGeogCS("WGS84");
 
578
            }
 
579
            else if (nEPSGDatumCode >= 6000)
 
580
            {
 
581
                char szName[32];
 
582
                sprintf( szName, "EPSG:%d", nEPSGDatumCode-2000 );
 
583
                oSRS.SetWellKnownGeogCS( szName );
 
584
                bHasSRS = TRUE;
 
585
            }
 
586
        }
 
587
 
 
588
        if (bHasUTMZone && ABS(nUTMZone) >= 1 && ABS(nUTMZone) <= 60)
 
589
        {
 
590
            bHasSRS = TRUE;
 
591
            oSRS.SetUTM(ABS(nUTMZone), nUTMZone > 0);
 
592
        }
 
593
        if (bHasSRS)
 
594
            oSRS.exportToWkt(&poDS->pszWKT);
 
595
    }
 
596
 
 
597
/* -------------------------------------------------------------------- */
 
598
/*      Create band information objects.                                */
 
599
/* -------------------------------------------------------------------- */
 
600
    poDS->nBands = 1;
 
601
    int i;
 
602
    for( i = 0; i < poDS->nBands; i++ )
 
603
    {
 
604
        poDS->SetBand( i+1, new HF2RasterBand( poDS, i+1, GDT_Float32 ) );
 
605
        poDS->GetRasterBand(i+1)->SetUnitType("m");
 
606
    }
 
607
 
 
608
    if (szApplicationName[0] != '\0')
 
609
        poDS->SetMetadataItem("APPLICATION_NAME", szApplicationName);
 
610
    poDS->SetMetadataItem("VERTICAL_PRECISION", CPLString().Printf("%f", fVertPres));
 
611
    if (bHasRelativePrecision)
 
612
        poDS->SetMetadataItem("RELATIVE_VERTICAL_PRECISION", CPLString().Printf("%f", fRelativePrecision));
 
613
 
 
614
/* -------------------------------------------------------------------- */
 
615
/*      Initialize any PAM information.                                 */
 
616
/* -------------------------------------------------------------------- */
 
617
    poDS->SetDescription( osOriginalFilename.c_str() );
 
618
    poDS->TryLoadXML();
 
619
 
 
620
/* -------------------------------------------------------------------- */
 
621
/*      Support overviews.                                              */
 
622
/* -------------------------------------------------------------------- */
 
623
    poDS->oOvManager.Initialize( poDS, osOriginalFilename.c_str() );
 
624
    return( poDS );
 
625
}
 
626
 
 
627
/************************************************************************/
 
628
/*                          GetGeoTransform()                           */
 
629
/************************************************************************/
 
630
 
 
631
CPLErr HF2Dataset::GetGeoTransform( double * padfTransform )
 
632
 
 
633
{
 
634
    memcpy(padfTransform, adfGeoTransform, 6 * sizeof(double));
 
635
 
 
636
    return( CE_None );
 
637
}
 
638
 
 
639
 
 
640
static void WriteShort(VSILFILE* fp, GInt16 val)
 
641
{
 
642
    CPL_LSBPTR16(&val);
 
643
    VSIFWriteL(&val, 2, 1, fp);
 
644
}
 
645
 
 
646
static void WriteInt(VSILFILE* fp, GInt32 val)
 
647
{
 
648
    CPL_LSBPTR32(&val);
 
649
    VSIFWriteL(&val, 4, 1, fp);
 
650
}
 
651
 
 
652
static void WriteFloat(VSILFILE* fp, float val)
 
653
{
 
654
    CPL_LSBPTR32(&val);
 
655
    VSIFWriteL(&val, 4, 1, fp);
 
656
}
 
657
 
 
658
static void WriteDouble(VSILFILE* fp, double val)
 
659
{
 
660
    CPL_LSBPTR64(&val);
 
661
    VSIFWriteL(&val, 8, 1, fp);
 
662
}
 
663
 
 
664
/************************************************************************/
 
665
/*                             CreateCopy()                             */
 
666
/************************************************************************/
 
667
 
 
668
GDALDataset* HF2Dataset::CreateCopy( const char * pszFilename,
 
669
                                     GDALDataset *poSrcDS, 
 
670
                                     int bStrict, char ** papszOptions, 
 
671
                                     GDALProgressFunc pfnProgress,
 
672
                                     void * pProgressData )
 
673
{
 
674
/* -------------------------------------------------------------------- */
 
675
/*      Some some rudimentary checks                                    */
 
676
/* -------------------------------------------------------------------- */
 
677
    int nBands = poSrcDS->GetRasterCount();
 
678
    if (nBands == 0)
 
679
    {
 
680
        CPLError( CE_Failure, CPLE_NotSupported, 
 
681
                  "HF2 driver does not support source dataset with zero band.\n");
 
682
        return NULL;
 
683
    }
 
684
 
 
685
    if (nBands != 1)
 
686
    {
 
687
        CPLError( (bStrict) ? CE_Failure : CE_Warning, CPLE_NotSupported, 
 
688
                  "HF2 driver only uses the first band of the dataset.\n");
 
689
        if (bStrict)
 
690
            return NULL;
 
691
    }
 
692
 
 
693
    if( pfnProgress && !pfnProgress( 0.0, NULL, pProgressData ) )
 
694
        return NULL;
 
695
 
 
696
/* -------------------------------------------------------------------- */
 
697
/*      Get source dataset info                                         */
 
698
/* -------------------------------------------------------------------- */
 
699
 
 
700
    int nXSize = poSrcDS->GetRasterXSize();
 
701
    int nYSize = poSrcDS->GetRasterYSize();
 
702
    double adfGeoTransform[6];
 
703
    poSrcDS->GetGeoTransform(adfGeoTransform);
 
704
    int bHasGeoTransform = !(adfGeoTransform[0] == 0 &&
 
705
                             adfGeoTransform[1] == 1 &&
 
706
                             adfGeoTransform[2] == 0 &&
 
707
                             adfGeoTransform[3] == 0 &&
 
708
                             adfGeoTransform[4] == 0 &&
 
709
                             adfGeoTransform[5] == 1);
 
710
    if (adfGeoTransform[2] != 0 || adfGeoTransform[4] != 0)
 
711
    {
 
712
        CPLError( CE_Failure, CPLE_NotSupported, 
 
713
                  "HF2 driver does not support CreateCopy() from skewed or rotated dataset.\n");
 
714
        return NULL;
 
715
    }
 
716
 
 
717
    GDALDataType eSrcDT = poSrcDS->GetRasterBand(1)->GetRasterDataType();
 
718
    GDALDataType eReqDT;
 
719
    float fVertPres = (float) 0.01;
 
720
    if (eSrcDT == GDT_Byte || eSrcDT == GDT_Int16)
 
721
    {
 
722
        fVertPres = 1;
 
723
        eReqDT = GDT_Int16;
 
724
    }
 
725
    else
 
726
        eReqDT = GDT_Float32;
 
727
 
 
728
/* -------------------------------------------------------------------- */
 
729
/*      Read creation options                                           */
 
730
/* -------------------------------------------------------------------- */
 
731
    const char* pszCompressed = CSLFetchNameValue(papszOptions, "COMPRESS");
 
732
    int bCompress = FALSE;
 
733
    if (pszCompressed)
 
734
        bCompress = CSLTestBoolean(pszCompressed);
 
735
    
 
736
    const char* pszVerticalPrecision = CSLFetchNameValue(papszOptions, "VERTICAL_PRECISION");
 
737
    if (pszVerticalPrecision)
 
738
    {
 
739
        fVertPres = (float) CPLAtofM(pszVerticalPrecision);
 
740
        if (fVertPres <= 0)
 
741
        {
 
742
            CPLError(CE_Warning, CPLE_AppDefined,
 
743
                     "Unsupported value for VERTICAL_PRECISION. Defaulting to 0.01");
 
744
            fVertPres = (float) 0.01;
 
745
        }
 
746
        if (eReqDT == GDT_Int16 && fVertPres > 1)
 
747
            eReqDT = GDT_Float32;
 
748
    }
 
749
 
 
750
    const char* pszBlockSize = CSLFetchNameValue(papszOptions, "BLOCKSIZE");
 
751
    int nTileSize = 256;
 
752
    if (pszBlockSize)
 
753
    {
 
754
        nTileSize = atoi(pszBlockSize);
 
755
        if (nTileSize < 8 || nTileSize > 4096)
 
756
        {
 
757
            CPLError(CE_Warning, CPLE_AppDefined,
 
758
                     "Unsupported value for BLOCKSIZE. Defaulting to 256");
 
759
            nTileSize = 256;
 
760
        }
 
761
    }
 
762
 
 
763
/* -------------------------------------------------------------------- */
 
764
/*      Parse source dataset georeferencing info                        */
 
765
/* -------------------------------------------------------------------- */
 
766
 
 
767
    int nExtendedHeaderLen = 0;
 
768
    if (bHasGeoTransform)
 
769
        nExtendedHeaderLen += 58;
 
770
    const char* pszProjectionRef = poSrcDS->GetProjectionRef();
 
771
    int nDatumCode = -2;
 
772
    int nUTMZone = 0;
 
773
    int bNorth = FALSE;
 
774
    int nEPSGCode = 0;
 
775
    int nExtentUnits = 1;
 
776
    if (pszProjectionRef != NULL && pszProjectionRef[0] != '\0')
 
777
    {
 
778
        OGRSpatialReference oSRS;
 
779
        char* pszTemp = (char*) pszProjectionRef;
 
780
        if (oSRS.importFromWkt(&pszTemp) == OGRERR_NONE)
 
781
        {
 
782
            const char* pszValue = NULL;
 
783
            if( oSRS.GetAuthorityName( "GEOGCS|DATUM" ) != NULL
 
784
                && EQUAL(oSRS.GetAuthorityName( "GEOGCS|DATUM" ),"EPSG") )
 
785
                nDatumCode = atoi(oSRS.GetAuthorityCode( "GEOGCS|DATUM" ));
 
786
            else if ((pszValue = oSRS.GetAttrValue("GEOGCS|DATUM")) != NULL)
 
787
            {
 
788
                if (strstr(pszValue, "WGS") && strstr(pszValue, "84"))
 
789
                    nDatumCode = 6326;
 
790
            }
 
791
 
 
792
            nUTMZone = oSRS.GetUTMZone(&bNorth);
 
793
        }
 
794
        if( oSRS.GetAuthorityName( "PROJCS" ) != NULL
 
795
            && EQUAL(oSRS.GetAuthorityName( "PROJCS" ),"EPSG") )
 
796
            nEPSGCode = atoi(oSRS.GetAuthorityCode( "PROJCS" ));
 
797
 
 
798
        if( oSRS.IsGeographic() )
 
799
        {
 
800
            nExtentUnits = 0;
 
801
        }
 
802
        else
 
803
        {
 
804
            double dfLinear = oSRS.GetLinearUnits();
 
805
 
 
806
            if( ABS(dfLinear - 0.3048) < 0.0000001 )
 
807
                nExtentUnits = 2;
 
808
            else if( ABS(dfLinear - atof(SRS_UL_US_FOOT_CONV)) < 0.00000001 )
 
809
                nExtentUnits = 3;
 
810
            else
 
811
                nExtentUnits = 1;
 
812
        }
 
813
    }
 
814
    if (nDatumCode != -2)
 
815
        nExtendedHeaderLen += 26;
 
816
    if (nUTMZone != 0)
 
817
        nExtendedHeaderLen += 26;
 
818
    if (nEPSGCode)
 
819
        nExtendedHeaderLen += 26;
 
820
 
 
821
/* -------------------------------------------------------------------- */
 
822
/*      Create target file                                              */
 
823
/* -------------------------------------------------------------------- */
 
824
 
 
825
    CPLString osFilename;
 
826
    if (bCompress)
 
827
    {
 
828
        osFilename = "/vsigzip/";
 
829
        osFilename += pszFilename;
 
830
    }
 
831
    else
 
832
        osFilename = pszFilename;
 
833
    VSILFILE* fp = VSIFOpenL(osFilename.c_str(), "wb");
 
834
    if (fp == NULL)
 
835
    {
 
836
        CPLError( CE_Failure, CPLE_AppDefined, 
 
837
                  "Cannot create %s", pszFilename );
 
838
        return NULL;
 
839
    }
 
840
 
 
841
/* -------------------------------------------------------------------- */
 
842
/*      Write header                                                    */
 
843
/* -------------------------------------------------------------------- */
 
844
 
 
845
    VSIFWriteL("HF2\0", 4, 1, fp);
 
846
    WriteShort(fp, 0);
 
847
    WriteInt(fp, nXSize);
 
848
    WriteInt(fp, nYSize);
 
849
    WriteShort(fp, (GInt16) nTileSize);
 
850
    WriteFloat(fp, fVertPres);
 
851
    float fHorizScale = (float) ((fabs(adfGeoTransform[1]) + fabs(adfGeoTransform[5])) / 2);
 
852
    WriteFloat(fp, fHorizScale);
 
853
    WriteInt(fp, nExtendedHeaderLen);
 
854
 
 
855
/* -------------------------------------------------------------------- */
 
856
/*      Write extended header                                           */
 
857
/* -------------------------------------------------------------------- */
 
858
 
 
859
    char szBlockName[16 + 1];
 
860
    if (bHasGeoTransform)
 
861
    {
 
862
        VSIFWriteL("bin\0", 4, 1, fp);
 
863
        memset(szBlockName, 0, 16 + 1);
 
864
        strcpy(szBlockName, "georef-extents");
 
865
        VSIFWriteL(szBlockName, 16, 1, fp);
 
866
        WriteInt(fp, 34);
 
867
        WriteShort(fp, (GInt16) nExtentUnits);
 
868
        WriteDouble(fp, adfGeoTransform[0]);
 
869
        WriteDouble(fp, adfGeoTransform[0] + nXSize * adfGeoTransform[1]);
 
870
        WriteDouble(fp, adfGeoTransform[3] + nYSize * adfGeoTransform[5]);
 
871
        WriteDouble(fp, adfGeoTransform[3]);
 
872
    }
 
873
    if (nUTMZone != 0)
 
874
    {
 
875
        VSIFWriteL("bin\0", 4, 1, fp);
 
876
        memset(szBlockName, 0, 16 + 1);
 
877
        strcpy(szBlockName, "georef-utm");
 
878
        VSIFWriteL(szBlockName, 16, 1, fp);
 
879
        WriteInt(fp, 2);
 
880
        WriteShort(fp, (GInt16) ((bNorth) ? nUTMZone : -nUTMZone));
 
881
    }
 
882
    if (nDatumCode != -2)
 
883
    {
 
884
        VSIFWriteL("bin\0", 4, 1, fp);
 
885
        memset(szBlockName, 0, 16 + 1);
 
886
        strcpy(szBlockName, "georef-datum");
 
887
        VSIFWriteL(szBlockName, 16, 1, fp);
 
888
        WriteInt(fp, 2);
 
889
        WriteShort(fp, (GInt16) nDatumCode);
 
890
    }
 
891
    if (nEPSGCode != 0)
 
892
    {
 
893
        VSIFWriteL("bin\0", 4, 1, fp);
 
894
        memset(szBlockName, 0, 16 + 1);
 
895
        strcpy(szBlockName, "georef-epsg-prj");
 
896
        VSIFWriteL(szBlockName, 16, 1, fp);
 
897
        WriteInt(fp, 2);
 
898
        WriteShort(fp, (GInt16) nEPSGCode);
 
899
    }
 
900
 
 
901
/* -------------------------------------------------------------------- */
 
902
/*      Copy imagery                                                    */
 
903
/* -------------------------------------------------------------------- */
 
904
    int nXBlocks = (nXSize + nTileSize - 1) / nTileSize;
 
905
    int nYBlocks = (nYSize + nTileSize - 1) / nTileSize;
 
906
 
 
907
    void* pTileBuffer = (void*) VSIMalloc(nTileSize * nTileSize * (GDALGetDataTypeSize(eReqDT) / 8));
 
908
    if (pTileBuffer == NULL)
 
909
    {
 
910
        CPLError( CE_Failure, CPLE_OutOfMemory, "Out of memory");
 
911
        VSIFCloseL(fp);
 
912
        return NULL;
 
913
    }
 
914
 
 
915
    int i, j, k, l;
 
916
    CPLErr eErr = CE_None;
 
917
    for(j=0;j<nYBlocks && eErr == CE_None;j++)
 
918
    {
 
919
        for(i=0;i<nXBlocks && eErr == CE_None;i++)
 
920
        {
 
921
            int nReqXSize = MIN(nTileSize, nXSize - i * nTileSize);
 
922
            int nReqYSize = MIN(nTileSize, nYSize - j * nTileSize);
 
923
            eErr = poSrcDS->GetRasterBand(1)->RasterIO(GF_Read,
 
924
                                                i * nTileSize, MAX(0, nYSize - (j + 1) * nTileSize),
 
925
                                                nReqXSize, nReqYSize,
 
926
                                                pTileBuffer, nReqXSize, nReqYSize,
 
927
                                                eReqDT, 0, 0);
 
928
            if (eErr != CE_None)
 
929
                break;
 
930
 
 
931
            if (eReqDT == GDT_Int16)
 
932
            {
 
933
                WriteFloat(fp, 1); /* scale */
 
934
                WriteFloat(fp, 0); /* offset */
 
935
                for(k=0;k<nReqYSize;k++)
 
936
                {
 
937
                    int nLastVal = ((short*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + 0];
 
938
                    GByte nWordSize = 1;
 
939
                    for(l=1;l<nReqXSize;l++)
 
940
                    {
 
941
                        int nVal = ((short*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + l];
 
942
                        int nDiff = nVal - nLastVal;
 
943
                        if (nDiff < -32768 || nDiff > 32767)
 
944
                        {
 
945
                            nWordSize = 4;
 
946
                            break;
 
947
                        }
 
948
                        if (nDiff < -128 || nDiff > 127)
 
949
                            nWordSize = 2;
 
950
                        nLastVal = nVal;
 
951
                    }
 
952
 
 
953
                    VSIFWriteL(&nWordSize, 1, 1, fp);
 
954
                    nLastVal = ((short*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + 0];
 
955
                    WriteInt(fp, nLastVal);
 
956
                    for(l=1;l<nReqXSize;l++)
 
957
                    {
 
958
                        int nVal = ((short*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + l];
 
959
                        int nDiff = nVal - nLastVal;
 
960
                        if (nWordSize == 1)
 
961
                        {
 
962
                            CPLAssert(nDiff >= -128 && nDiff <= 127);
 
963
                            char chDiff = (char)nDiff;
 
964
                            VSIFWriteL(&chDiff, 1, 1, fp);
 
965
                        }
 
966
                        else if (nWordSize == 2)
 
967
                        {
 
968
                            CPLAssert(nDiff >= -32768 && nDiff <= 32767);
 
969
                            WriteShort(fp, (short)nDiff);
 
970
                        }
 
971
                        else
 
972
                        {
 
973
                            WriteInt(fp, nDiff);
 
974
                        }
 
975
                        nLastVal = nVal;
 
976
                    }
 
977
                }
 
978
            }
 
979
            else
 
980
            {
 
981
                float fMinVal = ((float*)pTileBuffer)[0];
 
982
                float fMaxVal = fMinVal;
 
983
                for(k=1;k<nReqYSize*nReqXSize;k++)
 
984
                {
 
985
                    float fVal = ((float*)pTileBuffer)[k];
 
986
                    if (fVal < fMinVal) fMinVal = fVal;
 
987
                    if (fVal > fMaxVal) fMaxVal = fVal;
 
988
                }
 
989
 
 
990
                float fIntRange = (fMaxVal - fMinVal) / fVertPres;
 
991
                float fScale = (fMinVal == fMaxVal) ? 1 : (fMaxVal - fMinVal) / fIntRange;
 
992
                float fOffset = fMinVal;
 
993
                WriteFloat(fp, fScale); /* scale */
 
994
                WriteFloat(fp, fOffset); /* offset */
 
995
                for(k=0;k<nReqYSize;k++)
 
996
                {
 
997
                    float fLastVal = ((float*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + 0];
 
998
                    float fIntLastVal = (fLastVal - fOffset) / fScale;
 
999
                    CPLAssert(fIntLastVal >= -2147483648.0f && fIntLastVal <= 2147483647.0f);
 
1000
                    int nLastVal = (int)fIntLastVal;
 
1001
                    GByte nWordSize = 1;
 
1002
                    for(l=1;l<nReqXSize;l++)
 
1003
                    {
 
1004
                        float fVal = ((float*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + l];
 
1005
                        float fIntVal = (fVal - fOffset) / fScale;
 
1006
                        CPLAssert(fIntVal >= -2147483648.0f && fIntVal <= 2147483647.0f);
 
1007
                        int nVal = (int)fIntVal;
 
1008
                        int nDiff = nVal - nLastVal;
 
1009
                        CPLAssert((int)((GIntBig)nVal - nLastVal) == nDiff);
 
1010
                        if (nDiff < -32768 || nDiff > 32767)
 
1011
                        {
 
1012
                            nWordSize = 4;
 
1013
                            break;
 
1014
                        }
 
1015
                        if (nDiff < -128 || nDiff > 127)
 
1016
                            nWordSize = 2;
 
1017
                        nLastVal = nVal;
 
1018
                    }
 
1019
 
 
1020
                    VSIFWriteL(&nWordSize, 1, 1, fp);
 
1021
                    fLastVal = ((float*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + 0];
 
1022
                    fIntLastVal = (fLastVal - fOffset) / fScale;
 
1023
                    nLastVal = (int)fIntLastVal;
 
1024
                    WriteInt(fp, nLastVal);
 
1025
                    for(l=1;l<nReqXSize;l++)
 
1026
                    {
 
1027
                        float fVal = ((float*)pTileBuffer)[(nReqYSize - k - 1) * nReqXSize + l];
 
1028
                        float fIntVal = (fVal - fOffset) / fScale;
 
1029
                        int nVal = (int)fIntVal;
 
1030
                        int nDiff = nVal - nLastVal;
 
1031
                        CPLAssert((int)((GIntBig)nVal - nLastVal) == nDiff);
 
1032
                        if (nWordSize == 1)
 
1033
                        {
 
1034
                            CPLAssert(nDiff >= -128 && nDiff <= 127);
 
1035
                            char chDiff = (char)nDiff;
 
1036
                            VSIFWriteL(&chDiff, 1, 1, fp);
 
1037
                        }
 
1038
                        else if (nWordSize == 2)
 
1039
                        {
 
1040
                            CPLAssert(nDiff >= -32768 && nDiff <= 32767);
 
1041
                            WriteShort(fp, (short)nDiff);
 
1042
                        }
 
1043
                        else
 
1044
                        {
 
1045
                            WriteInt(fp, nDiff);
 
1046
                        }
 
1047
                        nLastVal = nVal;
 
1048
                    }
 
1049
                }
 
1050
            }
 
1051
 
 
1052
            if( pfnProgress && !pfnProgress( (j * nXBlocks + i + 1) * 1.0 / (nXBlocks * nYBlocks), NULL, pProgressData ) )
 
1053
            {
 
1054
                eErr = CE_Failure;
 
1055
                break;
 
1056
            }
 
1057
        }
 
1058
    }
 
1059
 
 
1060
    CPLFree(pTileBuffer);
 
1061
 
 
1062
    VSIFCloseL(fp);
 
1063
 
 
1064
    return (GDALDataset*) GDALOpen(osFilename.c_str(), GA_ReadOnly);
 
1065
}
 
1066
 
 
1067
/************************************************************************/
 
1068
/*                         GDALRegister_HF2()                           */
 
1069
/************************************************************************/
 
1070
 
 
1071
void GDALRegister_HF2()
 
1072
 
 
1073
{
 
1074
    GDALDriver  *poDriver;
 
1075
 
 
1076
    if( GDALGetDriverByName( "HF2" ) == NULL )
 
1077
    {
 
1078
        poDriver = new GDALDriver();
 
1079
        
 
1080
        poDriver->SetDescription( "HF2" );
 
1081
        poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, 
 
1082
                                   "HF2/HFZ heightfield raster" );
 
1083
        poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, 
 
1084
                                   "frmt_hf2.html" );
 
1085
        poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "hf2" );
 
1086
        poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
 
1087
"<CreationOptionList>"
 
1088
"   <Option name='VERTICAL_PRECISION' type='float' default='0.01' description='Vertical precision.'/>"
 
1089
"   <Option name='COMPRESS' type='boolean' default='false' description='Set to true to produce a GZip compressed file.'/>"
 
1090
"   <Option name='BLOCKSIZE' type='int' default='256' description='Tile size.'/>"
 
1091
"</CreationOptionList>");
 
1092
 
 
1093
        poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
 
1094
 
 
1095
        poDriver->pfnOpen = HF2Dataset::Open;
 
1096
        poDriver->pfnIdentify = HF2Dataset::Identify;
 
1097
        poDriver->pfnCreateCopy = HF2Dataset::CreateCopy;
 
1098
 
 
1099
        GetGDALDriverManager()->RegisterDriver( poDriver );
 
1100
    }
 
1101
}
 
1102