Coverage Report

Created: 2025-06-09 07:07

/src/gdal/frmts/leveller/levellerdataset.cpp
Line
Count
Source (jump to first uncovered line)
1
/******************************************************************************
2
 * levellerdataset.cpp,v 1.22
3
 *
4
 * Project:  Leveller TER Driver
5
 * Purpose:  Reader for Leveller TER documents
6
 * Author:   Ray Gardener, Daylon Graphics Ltd.
7
 *
8
 * Portions of this module derived from GDAL drivers by
9
 * Frank Warmerdam, see http://www.gdal.org
10
 *
11
 ******************************************************************************
12
 * Copyright (c) 2005-2007 Daylon Graphics Ltd.
13
 * Copyright (c) 2007-2013, Even Rouault <even dot rouault at spatialys.com>
14
 *
15
 * SPDX-License-Identifier: MIT
16
 ****************************************************************************/
17
18
#include "gdal_frmts.h"
19
#include "gdal_pam.h"
20
#include "ogr_spatialref.h"
21
22
static bool str_equal(const char *_s1, const char *_s2)
23
0
{
24
0
    return 0 == strcmp(_s1, _s2);
25
0
}
26
27
/*GDALDataset *LevellerCreateCopy( const char *, GDALDataset *, int, char **,
28
                                GDALProgressFunc pfnProgress,
29
                                void * pProgressData );
30
*/
31
32
/************************************************************************/
33
/* ==================================================================== */
34
/*                              LevellerDataset                         */
35
/* ==================================================================== */
36
/************************************************************************/
37
38
constexpr size_t kMaxTagNameLen = 63;
39
40
enum
41
{
42
    // Leveller coordsys types.
43
    LEV_COORDSYS_RASTER = 0,
44
    LEV_COORDSYS_LOCAL,
45
    LEV_COORDSYS_GEO
46
};
47
48
enum
49
{
50
    // Leveller digital axis extent styles.
51
    LEV_DA_POSITIONED = 0,
52
    LEV_DA_SIZED,
53
    LEV_DA_PIXEL_SIZED
54
};
55
56
typedef enum
57
{
58
    // Measurement unit IDs, OEM version.
59
    UNITLABEL_UNKNOWN = 0x00000000,
60
    UNITLABEL_PIXEL = 0x70780000,
61
    UNITLABEL_PERCENT = 0x25000000,
62
63
    UNITLABEL_RADIAN = 0x72616400,
64
    UNITLABEL_DEGREE = 0x64656700,
65
    UNITLABEL_ARCMINUTE = 0x6172636D,
66
    UNITLABEL_ARCSECOND = 0x61726373,
67
68
    UNITLABEL_YM = 0x796D0000,
69
    UNITLABEL_ZM = 0x7A6D0000,
70
    UNITLABEL_AM = 0x616D0000,
71
    UNITLABEL_FM = 0x666D0000,
72
    UNITLABEL_PM = 0x706D0000,
73
    UNITLABEL_A = 0x41000000,
74
    UNITLABEL_NM = 0x6E6D0000,
75
    UNITLABEL_U = 0x75000000,
76
    UNITLABEL_UM = 0x756D0000,
77
    UNITLABEL_PPT = 0x70707400,
78
    UNITLABEL_PT = 0x70740000,
79
    UNITLABEL_MM = 0x6D6D0000,
80
    UNITLABEL_P = 0x70000000,
81
    UNITLABEL_CM = 0x636D0000,
82
    UNITLABEL_IN = 0x696E0000,
83
    UNITLABEL_DFT = 0x64667400,
84
    UNITLABEL_DM = 0x646D0000,
85
    UNITLABEL_LI = 0x6C690000,
86
    UNITLABEL_SLI = 0x736C6900,
87
    UNITLABEL_SP = 0x73700000,
88
    UNITLABEL_FT = 0x66740000,
89
    UNITLABEL_SFT = 0x73667400,
90
    UNITLABEL_YD = 0x79640000,
91
    UNITLABEL_SYD = 0x73796400,
92
    UNITLABEL_M = 0x6D000000,
93
    UNITLABEL_FATH = 0x66617468,
94
    UNITLABEL_R = 0x72000000,
95
    UNITLABEL_RD = UNITLABEL_R,
96
    UNITLABEL_DAM = 0x64416D00,
97
    UNITLABEL_DKM = UNITLABEL_DAM,
98
    UNITLABEL_CH = 0x63680000,
99
    UNITLABEL_SCH = 0x73636800,
100
    UNITLABEL_HM = 0x686D0000,
101
    UNITLABEL_F = 0x66000000,
102
    UNITLABEL_KM = 0x6B6D0000,
103
    UNITLABEL_MI = 0x6D690000,
104
    UNITLABEL_SMI = 0x736D6900,
105
    UNITLABEL_NMI = 0x6E6D6900,
106
    UNITLABEL_MEGAM = 0x4D6D0000,
107
    UNITLABEL_LS = 0x6C730000,
108
    UNITLABEL_GM = 0x476D0000,
109
    UNITLABEL_LM = 0x6C6D0000,
110
    UNITLABEL_AU = 0x41550000,
111
    UNITLABEL_TM = 0x546D0000,
112
    UNITLABEL_LHR = 0x6C687200,
113
    UNITLABEL_LD = 0x6C640000,
114
    UNITLABEL_PETAM = 0x506D0000,
115
    UNITLABEL_LY = 0x6C790000,
116
    UNITLABEL_PC = 0x70630000,
117
    UNITLABEL_EXAM = 0x456D0000,
118
    UNITLABEL_KLY = 0x6B6C7900,
119
    UNITLABEL_KPC = 0x6B706300,
120
    UNITLABEL_ZETTAM = 0x5A6D0000,
121
    UNITLABEL_MLY = 0x4D6C7900,
122
    UNITLABEL_MPC = 0x4D706300,
123
    UNITLABEL_YOTTAM = 0x596D0000
124
} UNITLABEL;
125
126
typedef struct
127
{
128
    const char *pszID;
129
    double dScale;
130
    UNITLABEL oemCode;
131
} measurement_unit;
132
133
constexpr double kdays_per_year = 365.25;
134
constexpr double kdLStoM = 299792458.0;
135
constexpr double kdLYtoM = kdLStoM * kdays_per_year * 24 * 60 * 60;
136
constexpr double kdInch = 0.3048 / 12;
137
constexpr double kPI = M_PI;
138
139
constexpr int kFirstLinearMeasureIdx = 9;
140
141
static const measurement_unit kUnits[] = {
142
    {"", 1.0, UNITLABEL_UNKNOWN},
143
    {"px", 1.0, UNITLABEL_PIXEL},
144
    {"%", 1.0, UNITLABEL_PERCENT},  // not actually used
145
146
    {"rad", 1.0, UNITLABEL_RADIAN},
147
    {"\xB0", kPI / 180.0, UNITLABEL_DEGREE},  // \xB0 is Unicode degree symbol
148
    {"d", kPI / 180.0, UNITLABEL_DEGREE},
149
    {"deg", kPI / 180.0, UNITLABEL_DEGREE},
150
    {"'", kPI / (60.0 * 180.0), UNITLABEL_ARCMINUTE},
151
    {"\"", kPI / (3600.0 * 180.0), UNITLABEL_ARCSECOND},
152
153
    {"ym", 1.0e-24, UNITLABEL_YM},
154
    {"zm", 1.0e-21, UNITLABEL_ZM},
155
    {"am", 1.0e-18, UNITLABEL_AM},
156
    {"fm", 1.0e-15, UNITLABEL_FM},
157
    {"pm", 1.0e-12, UNITLABEL_PM},
158
    {"A", 1.0e-10, UNITLABEL_A},
159
    {"nm", 1.0e-9, UNITLABEL_NM},
160
    {"u", 1.0e-6, UNITLABEL_U},
161
    {"um", 1.0e-6, UNITLABEL_UM},
162
    {"ppt", kdInch / 72.27, UNITLABEL_PPT},
163
    {"pt", kdInch / 72.0, UNITLABEL_PT},
164
    {"mm", 1.0e-3, UNITLABEL_MM},
165
    {"p", kdInch / 6.0, UNITLABEL_P},
166
    {"cm", 1.0e-2, UNITLABEL_CM},
167
    {"in", kdInch, UNITLABEL_IN},
168
    {"dft", 0.03048, UNITLABEL_DFT},
169
    {"dm", 0.1, UNITLABEL_DM},
170
    {"li", 0.2011684 /* GDAL 0.20116684023368047 ? */, UNITLABEL_LI},
171
    {"sli", 0.201168402336805, UNITLABEL_SLI},
172
    {"sp", 0.2286, UNITLABEL_SP},
173
    {"ft", 0.3048, UNITLABEL_FT},
174
    {"sft", 1200.0 / 3937.0, UNITLABEL_SFT},
175
    {"yd", 0.9144, UNITLABEL_YD},
176
    {"syd", 0.914401828803658, UNITLABEL_SYD},
177
    {"m", 1.0, UNITLABEL_M},
178
    {"fath", 1.8288, UNITLABEL_FATH},
179
    {"rd", 5.02921, UNITLABEL_RD},
180
    {"dam", 10.0, UNITLABEL_DAM},
181
    {"dkm", 10.0, UNITLABEL_DKM},
182
    {"ch", 20.1168 /* GDAL: 2.0116684023368047 ? */, UNITLABEL_CH},
183
    {"sch", 20.1168402336805, UNITLABEL_SCH},
184
    {"hm", 100.0, UNITLABEL_HM},
185
    {"f", 201.168, UNITLABEL_F},
186
    {"km", 1000.0, UNITLABEL_KM},
187
    {"mi", 1609.344, UNITLABEL_MI},
188
    {"smi", 1609.34721869444, UNITLABEL_SMI},
189
    {"nmi", 1853.0, UNITLABEL_NMI},
190
    {"Mm", 1.0e+6, UNITLABEL_MEGAM},
191
    {"ls", kdLStoM, UNITLABEL_LS},
192
    {"Gm", 1.0e+9, UNITLABEL_GM},
193
    {"lm", kdLStoM * 60, UNITLABEL_LM},
194
    {"AU", 8.317 * kdLStoM * 60, UNITLABEL_AU},
195
    {"Tm", 1.0e+12, UNITLABEL_TM},
196
    {"lhr", 60.0 * 60.0 * kdLStoM, UNITLABEL_LHR},
197
    {"ld", 24 * 60.0 * 60.0 * kdLStoM, UNITLABEL_LD},
198
    {"Pm", 1.0e+15, UNITLABEL_PETAM},
199
    {"ly", kdLYtoM, UNITLABEL_LY},
200
    {"pc", 3.2616 * kdLYtoM, UNITLABEL_PC},
201
    {"Em", 1.0e+18, UNITLABEL_EXAM},
202
    {"kly", 1.0e+3 * kdLYtoM, UNITLABEL_KLY},
203
    {"kpc", 3.2616 * 1.0e+3 * kdLYtoM, UNITLABEL_KPC},
204
    {"Zm", 1.0e+21, UNITLABEL_ZETTAM},
205
    {"Mly", 1.0e+6 * kdLYtoM, UNITLABEL_MLY},
206
    {"Mpc", 3.2616 * 1.0e+6 * kdLYtoM, UNITLABEL_MPC},
207
    {"Ym", 1.0e+24, UNITLABEL_YOTTAM}};
208
209
// ----------------------------------------------------------------
210
211
static bool approx_equal(double a, double b)
212
0
{
213
0
    const double epsilon = 1e-5;
214
0
    return fabs(a - b) <= epsilon;
215
0
}
216
217
// ----------------------------------------------------------------
218
219
class LevellerRasterBand;
220
221
class LevellerDataset final : public GDALPamDataset
222
{
223
    friend class LevellerRasterBand;
224
    friend class digital_axis;
225
226
    int m_version;
227
228
    char *m_pszFilename;
229
    OGRSpatialReference m_oSRS{};
230
231
    // char              m_szUnits[8];
232
    char m_szElevUnits[8];
233
    double m_dElevScale;  // physical-to-logical scaling.
234
    double m_dElevBase;   // logical offset.
235
    double m_adfTransform[6];
236
    // double            m_dMeasurePerPixel;
237
    double m_dLogSpan[2];
238
239
    VSILFILE *m_fp;
240
    vsi_l_offset m_nDataOffset;
241
242
    bool load_from_file(VSILFILE *, const char *);
243
244
    static bool locate_data(vsi_l_offset &, size_t &, VSILFILE *, const char *);
245
    static bool get(int &, VSILFILE *, const char *);
246
247
    static bool get(size_t &n, VSILFILE *fp, const char *psz)
248
0
    {
249
0
        return get((int &)n, fp, psz);
250
0
    }
251
252
    static bool get(double &, VSILFILE *, const char *);
253
    static bool get(char *, size_t, VSILFILE *, const char *);
254
255
    bool write_header();
256
    bool write_tag(const char *, int);
257
    bool write_tag(const char *, size_t);
258
    bool write_tag(const char *, double);
259
    bool write_tag(const char *, const char *);
260
    bool write_tag_start(const char *, size_t);
261
    bool write(int);
262
    bool write(size_t);
263
    bool write(double);
264
    bool write_byte(size_t);
265
266
    static const measurement_unit *get_uom(const char *);
267
    static const measurement_unit *get_uom(UNITLABEL);
268
    static const measurement_unit *get_uom(double);
269
270
    static bool convert_measure(double, double &, const char *pszUnitsFrom);
271
    bool make_local_coordsys(const char *pszName, const char *pszUnits);
272
    bool make_local_coordsys(const char *pszName, UNITLABEL);
273
    const char *code_to_id(UNITLABEL) const;
274
    UNITLABEL id_to_code(const char *) const;
275
    UNITLABEL meter_measure_to_code(double) const;
276
    bool compute_elev_scaling(const OGRSpatialReference &);
277
    void raw_to_proj(double, double, double &, double &) const;
278
279
  public:
280
    LevellerDataset();
281
    virtual ~LevellerDataset();
282
283
    static GDALDataset *Open(GDALOpenInfo *);
284
    static int Identify(GDALOpenInfo *);
285
    static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
286
                               int nBandsIn, GDALDataType eType,
287
                               char **papszOptions);
288
289
    virtual CPLErr GetGeoTransform(double *) override;
290
291
    virtual CPLErr SetGeoTransform(double *) override;
292
293
    const OGRSpatialReference *GetSpatialRef() const override;
294
    CPLErr SetSpatialRef(const OGRSpatialReference *poSRS) override;
295
};
296
297
class digital_axis
298
{
299
  public:
300
0
    digital_axis() : m_eStyle(LEV_DA_PIXEL_SIZED), m_fixedEnd(0)
301
0
    {
302
0
        m_d[0] = 0.0;
303
0
        m_d[1] = 0.0;
304
0
    }
305
306
    bool get(LevellerDataset &ds, VSILFILE *fp, int n)
307
0
    {
308
0
        char szTag[32];
309
0
        snprintf(szTag, sizeof(szTag), "coordsys_da%d_style", n);
310
0
        if (!ds.get(m_eStyle, fp, szTag))
311
0
            return false;
312
0
        snprintf(szTag, sizeof(szTag), "coordsys_da%d_fixedend", n);
313
0
        if (!ds.get(m_fixedEnd, fp, szTag))
314
0
            return false;
315
0
        snprintf(szTag, sizeof(szTag), "coordsys_da%d_v0", n);
316
0
        if (!ds.get(m_d[0], fp, szTag))
317
0
            return false;
318
0
        snprintf(szTag, sizeof(szTag), "coordsys_da%d_v1", n);
319
0
        if (!ds.get(m_d[1], fp, szTag))
320
0
            return false;
321
0
        return true;
322
0
    }
323
324
    double origin(size_t pixels) const
325
0
    {
326
0
        if (m_fixedEnd == 1)
327
0
        {
328
0
            switch (m_eStyle)
329
0
            {
330
0
                case LEV_DA_SIZED:
331
0
                    return m_d[1] + m_d[0];
332
333
0
                case LEV_DA_PIXEL_SIZED:
334
0
                    return m_d[1] + (m_d[0] * (pixels - 1));
335
0
            }
336
0
        }
337
0
        return m_d[0];
338
0
    }
339
340
    double scaling(size_t pixels) const
341
0
    {
342
0
        CPLAssert(pixels > 1);
343
0
        if (m_eStyle == LEV_DA_PIXEL_SIZED)
344
0
            return m_d[1 - m_fixedEnd];
345
346
0
        return this->length(static_cast<int>(pixels)) / (pixels - 1);
347
0
    }
348
349
    double length(int pixels) const
350
0
    {
351
        // Return the signed length of the axis.
352
353
0
        switch (m_eStyle)
354
0
        {
355
0
            case LEV_DA_POSITIONED:
356
0
                return m_d[1] - m_d[0];
357
358
0
            case LEV_DA_SIZED:
359
0
                return m_d[1 - m_fixedEnd];
360
361
0
            case LEV_DA_PIXEL_SIZED:
362
0
                return m_d[1 - m_fixedEnd] * (pixels - 1);
363
0
        }
364
0
        CPLAssert(false);
365
0
        return 0.0;
366
0
    }
367
368
  protected:
369
    int m_eStyle;
370
    int m_fixedEnd;
371
    double m_d[2];
372
};
373
374
/************************************************************************/
375
/* ==================================================================== */
376
/*                            LevellerRasterBand                        */
377
/* ==================================================================== */
378
/************************************************************************/
379
380
class LevellerRasterBand final : public GDALPamRasterBand
381
{
382
    friend class LevellerDataset;
383
384
    float *m_pLine;
385
    bool m_bFirstTime;
386
387
  public:
388
    explicit LevellerRasterBand(LevellerDataset *);
389
    virtual ~LevellerRasterBand();
390
391
    bool Init();
392
393
    // Geomeasure support.
394
    virtual const char *GetUnitType() override;
395
    virtual double GetScale(int *pbSuccess = nullptr) override;
396
    virtual double GetOffset(int *pbSuccess = nullptr) override;
397
398
    virtual CPLErr IReadBlock(int, int, void *) override;
399
    virtual CPLErr IWriteBlock(int, int, void *) override;
400
    virtual CPLErr SetUnitType(const char *) override;
401
};
402
403
/************************************************************************/
404
/*                         LevellerRasterBand()                         */
405
/************************************************************************/
406
407
LevellerRasterBand::LevellerRasterBand(LevellerDataset *poDSIn)
408
0
    : m_pLine(nullptr), m_bFirstTime(true)
409
0
{
410
0
    poDS = poDSIn;
411
0
    nBand = 1;
412
413
0
    eDataType = GDT_Float32;
414
415
0
    nBlockXSize = poDS->GetRasterXSize();
416
0
    nBlockYSize = 1;  // poDS->GetRasterYSize();
417
0
}
418
419
/************************************************************************/
420
/*                           Init()                                     */
421
/************************************************************************/
422
423
bool LevellerRasterBand::Init()
424
0
{
425
0
    m_pLine = reinterpret_cast<float *>(
426
0
        VSI_MALLOC2_VERBOSE(sizeof(float), nBlockXSize));
427
0
    return m_pLine != nullptr;
428
0
}
429
430
LevellerRasterBand::~LevellerRasterBand()
431
0
{
432
0
    CPLFree(m_pLine);
433
0
}
434
435
/************************************************************************/
436
/*                             IWriteBlock()                            */
437
/************************************************************************/
438
439
CPLErr LevellerRasterBand::IWriteBlock(CPL_UNUSED int nBlockXOff,
440
                                       int nBlockYOff, void *pImage)
441
0
{
442
0
    CPLAssert(nBlockXOff == 0);
443
0
    CPLAssert(pImage != nullptr);
444
0
    CPLAssert(m_pLine != nullptr);
445
446
    /*  #define sgn(_n) ((_n) < 0 ? -1 : ((_n) > 0 ? 1 : 0) )
447
        #define sround(_f)                          \
448
        (int)((_f) + (0.5 * sgn(_f)))
449
    */
450
0
    const size_t pixelsize = sizeof(float);
451
452
0
    LevellerDataset &ds = *reinterpret_cast<LevellerDataset *>(poDS);
453
0
    if (m_bFirstTime)
454
0
    {
455
0
        m_bFirstTime = false;
456
0
        if (!ds.write_header())
457
0
            return CE_Failure;
458
0
        ds.m_nDataOffset = VSIFTellL(ds.m_fp);
459
0
    }
460
0
    const size_t rowbytes = nBlockXSize * pixelsize;
461
0
    const float *pfImage = reinterpret_cast<float *>(pImage);
462
463
0
    if (0 ==
464
0
        VSIFSeekL(ds.m_fp, ds.m_nDataOffset + nBlockYOff * rowbytes, SEEK_SET))
465
0
    {
466
0
        for (size_t x = 0; x < (size_t)nBlockXSize; x++)
467
0
        {
468
            // Convert logical elevations to physical.
469
0
            m_pLine[x] = static_cast<float>((pfImage[x] - ds.m_dElevBase) /
470
0
                                            ds.m_dElevScale);
471
0
        }
472
473
#ifdef CPL_MSB
474
        GDALSwapWords(m_pLine, pixelsize, nBlockXSize, pixelsize);
475
#endif
476
0
        if (1 == VSIFWriteL(m_pLine, rowbytes, 1, ds.m_fp))
477
0
            return CE_None;
478
0
    }
479
480
0
    return CE_Failure;
481
0
}
482
483
CPLErr LevellerRasterBand::SetUnitType(const char *psz)
484
0
{
485
0
    LevellerDataset &ds = *reinterpret_cast<LevellerDataset *>(poDS);
486
487
0
    if (strlen(psz) >= sizeof(ds.m_szElevUnits))
488
0
        return CE_Failure;
489
490
0
    strcpy(ds.m_szElevUnits, psz);
491
492
0
    return CE_None;
493
0
}
494
495
/************************************************************************/
496
/*                             IReadBlock()                             */
497
/************************************************************************/
498
499
CPLErr LevellerRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
500
                                      void *pImage)
501
502
0
{
503
0
    CPLAssert(sizeof(float) == sizeof(GInt32));
504
0
    CPLAssert(nBlockXOff == 0);
505
0
    CPLAssert(pImage != nullptr);
506
507
0
    LevellerDataset *poGDS = reinterpret_cast<LevellerDataset *>(poDS);
508
509
    /* -------------------------------------------------------------------- */
510
    /*      Seek to scanline.                                               */
511
    /* -------------------------------------------------------------------- */
512
0
    const size_t rowbytes = nBlockXSize * sizeof(float);
513
514
0
    if (0 != VSIFSeekL(poGDS->m_fp,
515
0
                       poGDS->m_nDataOffset + nBlockYOff * rowbytes, SEEK_SET))
516
0
    {
517
0
        CPLError(CE_Failure, CPLE_FileIO, "Leveller seek failed: %s",
518
0
                 VSIStrerror(errno));
519
0
        return CE_Failure;
520
0
    }
521
522
    /* -------------------------------------------------------------------- */
523
    /*      Read the scanline into the image buffer.                        */
524
    /* -------------------------------------------------------------------- */
525
526
0
    if (VSIFReadL(pImage, rowbytes, 1, poGDS->m_fp) != 1)
527
0
    {
528
0
        CPLError(CE_Failure, CPLE_FileIO, "Leveller read failed: %s",
529
0
                 VSIStrerror(errno));
530
0
        return CE_Failure;
531
0
    }
532
533
/* -------------------------------------------------------------------- */
534
/*      Swap on MSB platforms.                                          */
535
/* -------------------------------------------------------------------- */
536
#ifdef CPL_MSB
537
    GDALSwapWords(pImage, 4, nRasterXSize, 4);
538
#endif
539
540
    /* -------------------------------------------------------------------- */
541
    /*      Convert from legacy-format fixed-point if necessary.            */
542
    /* -------------------------------------------------------------------- */
543
0
    float *pf = reinterpret_cast<float *>(pImage);
544
545
0
    if (poGDS->m_version < 6)
546
0
    {
547
0
        GInt32 *pi = reinterpret_cast<int *>(pImage);
548
0
        for (size_t i = 0; i < (size_t)nBlockXSize; i++)
549
0
            pf[i] = static_cast<float>(pi[i]) / 65536;
550
0
    }
551
552
/* -------------------------------------------------------------------- */
553
/*      Convert raw elevations to realworld elevs.                      */
554
/* -------------------------------------------------------------------- */
555
#if 0
556
    for(size_t i = 0; i < nBlockXSize; i++)
557
        pf[i] *= poGDS->m_dWorldscale; //this->GetScale();
558
#endif
559
560
0
    return CE_None;
561
0
}
562
563
/************************************************************************/
564
/*                            GetUnitType()                             */
565
/************************************************************************/
566
const char *LevellerRasterBand::GetUnitType()
567
0
{
568
    // Return elevation units.
569
570
0
    LevellerDataset *poGDS = reinterpret_cast<LevellerDataset *>(poDS);
571
572
0
    return poGDS->m_szElevUnits;
573
0
}
574
575
/************************************************************************/
576
/*                              GetScale()                              */
577
/************************************************************************/
578
579
double LevellerRasterBand::GetScale(int *pbSuccess)
580
0
{
581
0
    LevellerDataset *poGDS = reinterpret_cast<LevellerDataset *>(poDS);
582
0
    if (pbSuccess != nullptr)
583
0
        *pbSuccess = TRUE;
584
0
    return poGDS->m_dElevScale;
585
0
}
586
587
/************************************************************************/
588
/*                             GetOffset()                              */
589
/************************************************************************/
590
591
double LevellerRasterBand::GetOffset(int *pbSuccess)
592
0
{
593
0
    LevellerDataset *poGDS = reinterpret_cast<LevellerDataset *>(poDS);
594
0
    if (pbSuccess != nullptr)
595
0
        *pbSuccess = TRUE;
596
0
    return poGDS->m_dElevBase;
597
0
}
598
599
/************************************************************************/
600
/* ==================================================================== */
601
/*                              LevellerDataset                         */
602
/* ==================================================================== */
603
/************************************************************************/
604
605
/************************************************************************/
606
/*                          LevellerDataset()                           */
607
/************************************************************************/
608
609
LevellerDataset::LevellerDataset()
610
0
    : m_version(0), m_pszFilename(nullptr), m_dElevScale(), m_dElevBase(),
611
0
      m_fp(nullptr), m_nDataOffset()
612
0
{
613
0
    m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
614
0
    memset(m_szElevUnits, 0, sizeof(m_szElevUnits));
615
0
    memset(m_adfTransform, 0, sizeof(m_adfTransform));
616
0
    memset(m_dLogSpan, 0, sizeof(m_dLogSpan));
617
0
}
618
619
/************************************************************************/
620
/*                          ~LevellerDataset()                          */
621
/************************************************************************/
622
623
LevellerDataset::~LevellerDataset()
624
0
{
625
0
    FlushCache(true);
626
627
0
    CPLFree(m_pszFilename);
628
629
0
    if (m_fp != nullptr)
630
0
        VSIFCloseL(m_fp);
631
0
}
632
633
static double degrees_to_radians(double d)
634
0
{
635
0
    return d * (M_PI / 180);
636
0
}
637
638
static double average(double a, double b)
639
0
{
640
0
    return 0.5 * (a + b);
641
0
}
642
643
void LevellerDataset::raw_to_proj(double x, double y, double &xp,
644
                                  double &yp) const
645
0
{
646
0
    xp = x * m_adfTransform[1] + m_adfTransform[0];
647
0
    yp = y * m_adfTransform[5] + m_adfTransform[3];
648
0
}
649
650
bool LevellerDataset::compute_elev_scaling(const OGRSpatialReference &sr)
651
0
{
652
0
    const char *pszGroundUnits = nullptr;
653
654
0
    if (!sr.IsGeographic())
655
0
    {
656
        // For projected or local CS, the elev scale is
657
        // the average ground scale.
658
0
        m_dElevScale = average(m_adfTransform[1], m_adfTransform[5]);
659
660
0
        const double dfLinear = sr.GetLinearUnits();
661
0
        const measurement_unit *pu = this->get_uom(dfLinear);
662
0
        if (pu == nullptr)
663
0
            return false;
664
665
0
        pszGroundUnits = pu->pszID;
666
0
    }
667
0
    else
668
0
    {
669
0
        pszGroundUnits = "m";
670
671
0
        const double kdEarthCircumPolar = 40007849;
672
0
        const double kdEarthCircumEquat = 40075004;
673
674
0
        const double xr = 0.5 * nRasterXSize;
675
0
        const double yr = 0.5 * nRasterYSize;
676
677
0
        double xg[2], yg[2];
678
0
        raw_to_proj(xr, yr, xg[0], yg[0]);
679
0
        raw_to_proj(xr + 1, yr + 1, xg[1], yg[1]);
680
681
        // The earths' circumference shrinks using a sin()
682
        // curve as we go up in latitude.
683
0
        const double dLatCircum =
684
0
            kdEarthCircumEquat * sin(degrees_to_radians(90.0 - yg[0]));
685
686
        // Derive meter distance between geolongitudes
687
        // in xg[0] and xg[1].
688
0
        const double dx = fabs(xg[1] - xg[0]) / 360.0 * dLatCircum;
689
0
        const double dy = fabs(yg[1] - yg[0]) / 360.0 * kdEarthCircumPolar;
690
691
0
        m_dElevScale = average(dx, dy);
692
0
    }
693
694
0
    m_dElevBase = m_dLogSpan[0];
695
696
    // Convert from ground units to elev units.
697
0
    const measurement_unit *puG = this->get_uom(pszGroundUnits);
698
0
    const measurement_unit *puE = this->get_uom(m_szElevUnits);
699
700
0
    if (puG == nullptr || puE == nullptr)
701
0
        return false;
702
703
0
    const double g_to_e = puG->dScale / puE->dScale;
704
705
0
    m_dElevScale *= g_to_e;
706
0
    return true;
707
0
}
708
709
bool LevellerDataset::write_header()
710
0
{
711
0
    char szHeader[5];
712
0
    strcpy(szHeader, "trrn");
713
0
    szHeader[4] = 7;  // TER v7 introduced w/ Lev 2.6.
714
715
0
    if (1 != VSIFWriteL(szHeader, 5, 1, m_fp) ||
716
0
        !this->write_tag("hf_w", (size_t)nRasterXSize) ||
717
0
        !this->write_tag("hf_b", (size_t)nRasterYSize))
718
0
    {
719
0
        CPLError(CE_Failure, CPLE_FileIO, "Could not write header");
720
0
        return false;
721
0
    }
722
723
0
    m_dElevBase = 0.0;
724
0
    m_dElevScale = 1.0;
725
726
0
    if (m_oSRS.IsEmpty())
727
0
    {
728
0
        write_tag("csclass", LEV_COORDSYS_RASTER);
729
0
    }
730
0
    else
731
0
    {
732
0
        char *pszWkt = nullptr;
733
0
        m_oSRS.exportToWkt(&pszWkt);
734
0
        if (pszWkt)
735
0
            write_tag("coordsys_wkt", pszWkt);
736
0
        CPLFree(pszWkt);
737
0
        const UNITLABEL units_elev = this->id_to_code(m_szElevUnits);
738
739
0
        const int bHasECS =
740
0
            (units_elev != UNITLABEL_PIXEL && units_elev != UNITLABEL_UNKNOWN);
741
742
0
        write_tag("coordsys_haselevm", bHasECS);
743
744
0
        if (bHasECS)
745
0
        {
746
0
            if (!this->compute_elev_scaling(m_oSRS))
747
0
                return false;
748
749
            // Raw-to-real units scaling.
750
0
            write_tag("coordsys_em_scale", m_dElevScale);
751
752
            // Elev offset, in real units.
753
0
            write_tag("coordsys_em_base", m_dElevBase);
754
0
            write_tag("coordsys_em_units", units_elev);
755
0
        }
756
757
0
        if (m_oSRS.IsLocal())
758
0
        {
759
0
            write_tag("csclass", LEV_COORDSYS_LOCAL);
760
761
0
            const double dfLinear = m_oSRS.GetLinearUnits();
762
0
            const int n = this->meter_measure_to_code(dfLinear);
763
0
            write_tag("coordsys_units", n);
764
0
        }
765
0
        else
766
0
        {
767
0
            write_tag("csclass", LEV_COORDSYS_GEO);
768
0
        }
769
770
0
        if (m_adfTransform[2] != 0.0 || m_adfTransform[4] != 0.0)
771
0
        {
772
0
            CPLError(CE_Failure, CPLE_IllegalArg,
773
0
                     "Cannot handle rotated geotransform");
774
0
            return false;
775
0
        }
776
777
        // todo: GDAL gridpost spacing is based on extent / rastersize
778
        // instead of extent / (rastersize-1) like Leveller.
779
        // We need to look into this and adjust accordingly.
780
781
        // Write north-south digital axis.
782
0
        write_tag("coordsys_da0_style", LEV_DA_PIXEL_SIZED);
783
0
        write_tag("coordsys_da0_fixedend", 0);
784
0
        write_tag("coordsys_da0_v0", m_adfTransform[3]);
785
0
        write_tag("coordsys_da0_v1", m_adfTransform[5]);
786
787
        // Write east-west digital axis.
788
0
        write_tag("coordsys_da1_style", LEV_DA_PIXEL_SIZED);
789
0
        write_tag("coordsys_da1_fixedend", 0);
790
0
        write_tag("coordsys_da1_v0", m_adfTransform[0]);
791
0
        write_tag("coordsys_da1_v1", m_adfTransform[1]);
792
0
    }
793
794
0
    this->write_tag_start("hf_data",
795
0
                          sizeof(float) * nRasterXSize * nRasterYSize);
796
797
0
    return true;
798
0
}
799
800
/************************************************************************/
801
/*                          SetGeoTransform()                           */
802
/************************************************************************/
803
804
CPLErr LevellerDataset::SetGeoTransform(double *padfGeoTransform)
805
0
{
806
0
    memcpy(m_adfTransform, padfGeoTransform, sizeof(m_adfTransform));
807
808
0
    return CE_None;
809
0
}
810
811
/************************************************************************/
812
/*                           SetSpatialRef()                            */
813
/************************************************************************/
814
815
CPLErr LevellerDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
816
0
{
817
0
    m_oSRS.Clear();
818
0
    if (poSRS)
819
0
        m_oSRS = *poSRS;
820
821
0
    return CE_None;
822
0
}
823
824
/************************************************************************/
825
/*                           Create()                                   */
826
/************************************************************************/
827
GDALDataset *LevellerDataset::Create(const char *pszFilename, int nXSize,
828
                                     int nYSize, int nBandsIn,
829
                                     GDALDataType eType, char **papszOptions)
830
0
{
831
0
    if (nBandsIn != 1)
832
0
    {
833
0
        CPLError(CE_Failure, CPLE_IllegalArg, "Band count must be 1");
834
0
        return nullptr;
835
0
    }
836
837
0
    if (eType != GDT_Float32)
838
0
    {
839
0
        CPLError(CE_Failure, CPLE_IllegalArg, "Pixel type must be Float32");
840
0
        return nullptr;
841
0
    }
842
843
0
    if (nXSize < 2 || nYSize < 2)
844
0
    {
845
0
        CPLError(CE_Failure, CPLE_IllegalArg,
846
0
                 "One or more raster dimensions too small");
847
0
        return nullptr;
848
0
    }
849
850
0
    LevellerDataset *poDS = new LevellerDataset;
851
852
0
    poDS->eAccess = GA_Update;
853
854
0
    poDS->m_pszFilename = CPLStrdup(pszFilename);
855
856
0
    poDS->m_fp = VSIFOpenL(pszFilename, "wb+");
857
858
0
    if (poDS->m_fp == nullptr)
859
0
    {
860
0
        CPLError(CE_Failure, CPLE_OpenFailed,
861
0
                 "Attempt to create file `%s' failed.", pszFilename);
862
0
        delete poDS;
863
0
        return nullptr;
864
0
    }
865
866
    // Header will be written the first time IWriteBlock
867
    // is called.
868
869
0
    poDS->nRasterXSize = nXSize;
870
0
    poDS->nRasterYSize = nYSize;
871
872
0
    const char *pszValue = CSLFetchNameValue(papszOptions, "MINUSERPIXELVALUE");
873
0
    if (pszValue != nullptr)
874
0
        poDS->m_dLogSpan[0] = CPLAtof(pszValue);
875
0
    else
876
0
    {
877
0
        delete poDS;
878
0
        CPLError(CE_Failure, CPLE_IllegalArg,
879
0
                 "MINUSERPIXELVALUE must be specified.");
880
0
        return nullptr;
881
0
    }
882
883
0
    pszValue = CSLFetchNameValue(papszOptions, "MAXUSERPIXELVALUE");
884
0
    if (pszValue != nullptr)
885
0
        poDS->m_dLogSpan[1] = CPLAtof(pszValue);
886
887
0
    if (poDS->m_dLogSpan[1] < poDS->m_dLogSpan[0])
888
0
    {
889
0
        double t = poDS->m_dLogSpan[0];
890
0
        poDS->m_dLogSpan[0] = poDS->m_dLogSpan[1];
891
0
        poDS->m_dLogSpan[1] = t;
892
0
    }
893
894
    // --------------------------------------------------------------------
895
    //      Instance a band.
896
    // --------------------------------------------------------------------
897
0
    LevellerRasterBand *poBand = new LevellerRasterBand(poDS);
898
0
    poDS->SetBand(1, poBand);
899
900
0
    if (!poBand->Init())
901
0
    {
902
0
        delete poDS;
903
0
        return nullptr;
904
0
    }
905
906
0
    return poDS;
907
0
}
908
909
bool LevellerDataset::write_byte(size_t n)
910
0
{
911
0
    unsigned char uch = static_cast<unsigned char>(n);
912
0
    return 1 == VSIFWriteL(&uch, 1, 1, m_fp);
913
0
}
914
915
bool LevellerDataset::write(int n)
916
0
{
917
0
    CPL_LSBPTR32(&n);
918
0
    return 1 == VSIFWriteL(&n, sizeof(n), 1, m_fp);
919
0
}
920
921
bool LevellerDataset::write(size_t n)
922
0
{
923
0
    GUInt32 n32 = (GUInt32)n;
924
0
    CPL_LSBPTR32(&n32);
925
0
    return 1 == VSIFWriteL(&n32, sizeof(n32), 1, m_fp);
926
0
}
927
928
bool LevellerDataset::write(double d)
929
0
{
930
0
    CPL_LSBPTR64(&d);
931
0
    return 1 == VSIFWriteL(&d, sizeof(d), 1, m_fp);
932
0
}
933
934
bool LevellerDataset::write_tag_start(const char *pszTag, size_t n)
935
0
{
936
0
    if (this->write_byte(strlen(pszTag)))
937
0
    {
938
0
        return (1 == VSIFWriteL(pszTag, strlen(pszTag), 1, m_fp) &&
939
0
                this->write(n));
940
0
    }
941
942
0
    return false;
943
0
}
944
945
bool LevellerDataset::write_tag(const char *pszTag, int n)
946
0
{
947
0
    return (this->write_tag_start(pszTag, sizeof(n)) && this->write(n));
948
0
}
949
950
bool LevellerDataset::write_tag(const char *pszTag, size_t n)
951
0
{
952
0
    return (this->write_tag_start(pszTag, sizeof(n)) && this->write(n));
953
0
}
954
955
bool LevellerDataset::write_tag(const char *pszTag, double d)
956
0
{
957
0
    return (this->write_tag_start(pszTag, sizeof(d)) && this->write(d));
958
0
}
959
960
bool LevellerDataset::write_tag(const char *pszTag, const char *psz)
961
0
{
962
0
    CPLAssert(strlen(pszTag) <= kMaxTagNameLen);
963
964
0
    char sz[kMaxTagNameLen + 1];
965
0
    snprintf(sz, sizeof(sz), "%sl", pszTag);
966
0
    const size_t len = strlen(psz);
967
968
0
    if (len > 0 && this->write_tag(sz, len))
969
0
    {
970
0
        snprintf(sz, sizeof(sz), "%sd", pszTag);
971
0
        this->write_tag_start(sz, len);
972
0
        return 1 == VSIFWriteL(psz, len, 1, m_fp);
973
0
    }
974
0
    return false;
975
0
}
976
977
bool LevellerDataset::locate_data(vsi_l_offset &offset, size_t &len,
978
                                  VSILFILE *fp, const char *pszTag)
979
0
{
980
    // Locate the file offset of the desired tag's data.
981
    // If it is not available, return false.
982
    // If the tag is found, leave the filemark at the
983
    // start of its data.
984
985
0
    if (0 != VSIFSeekL(fp, 5, SEEK_SET))
986
0
        return false;
987
988
0
    const int kMaxDescLen = 64;
989
0
    for (;;)
990
0
    {
991
0
        unsigned char c;
992
0
        if (1 != VSIFReadL(&c, sizeof(c), 1, fp))
993
0
            return false;
994
995
0
        const size_t descriptorLen = c;
996
0
        if (descriptorLen == 0 || descriptorLen > (size_t)kMaxDescLen)
997
0
            return false;
998
999
0
        char descriptor[kMaxDescLen + 1];
1000
0
        if (1 != VSIFReadL(descriptor, descriptorLen, 1, fp))
1001
0
            return false;
1002
1003
0
        GUInt32 datalen;
1004
0
        if (1 != VSIFReadL(&datalen, sizeof(datalen), 1, fp))
1005
0
            return false;
1006
1007
0
        CPL_LSBPTR32(&datalen);
1008
0
        descriptor[descriptorLen] = 0;
1009
0
        if (str_equal(descriptor, pszTag))
1010
0
        {
1011
0
            len = (size_t)datalen;
1012
0
            offset = VSIFTellL(fp);
1013
0
            return true;
1014
0
        }
1015
0
        else
1016
0
        {
1017
            // Seek to next tag.
1018
0
            if (0 != VSIFSeekL(fp, (vsi_l_offset)datalen, SEEK_CUR))
1019
0
                return false;
1020
0
        }
1021
0
    }
1022
0
}
1023
1024
/************************************************************************/
1025
/*                                get()                                 */
1026
/************************************************************************/
1027
1028
bool LevellerDataset::get(int &n, VSILFILE *fp, const char *psz)
1029
0
{
1030
0
    vsi_l_offset offset;
1031
0
    size_t len;
1032
1033
0
    if (locate_data(offset, len, fp, psz))
1034
0
    {
1035
0
        GInt32 value;
1036
0
        if (1 == VSIFReadL(&value, sizeof(value), 1, fp))
1037
0
        {
1038
0
            CPL_LSBPTR32(&value);
1039
0
            n = static_cast<int>(value);
1040
0
            return true;
1041
0
        }
1042
0
    }
1043
0
    return false;
1044
0
}
1045
1046
/************************************************************************/
1047
/*                                get()                                 */
1048
/************************************************************************/
1049
1050
bool LevellerDataset::get(double &d, VSILFILE *fp, const char *pszTag)
1051
0
{
1052
0
    vsi_l_offset offset;
1053
0
    size_t len;
1054
1055
0
    if (locate_data(offset, len, fp, pszTag))
1056
0
    {
1057
0
        if (1 == VSIFReadL(&d, sizeof(d), 1, fp))
1058
0
        {
1059
0
            CPL_LSBPTR64(&d);
1060
0
            return true;
1061
0
        }
1062
0
    }
1063
0
    return false;
1064
0
}
1065
1066
/************************************************************************/
1067
/*                                get()                                 */
1068
/************************************************************************/
1069
bool LevellerDataset::get(char *pszValue, size_t maxchars, VSILFILE *fp,
1070
                          const char *pszTag)
1071
0
{
1072
0
    char szTag[65];
1073
1074
    // We can assume 8-bit encoding, so just go straight
1075
    // to the *_d tag.
1076
0
    snprintf(szTag, sizeof(szTag), "%sd", pszTag);
1077
1078
0
    vsi_l_offset offset;
1079
0
    size_t len;
1080
1081
0
    if (locate_data(offset, len, fp, szTag))
1082
0
    {
1083
0
        if (len > maxchars)
1084
0
            return false;
1085
1086
0
        if (1 == VSIFReadL(pszValue, len, 1, fp))
1087
0
        {
1088
0
            pszValue[len] = 0;  // terminate C-string
1089
0
            return true;
1090
0
        }
1091
0
    }
1092
1093
0
    return false;
1094
0
}
1095
1096
UNITLABEL LevellerDataset::meter_measure_to_code(double dM) const
1097
0
{
1098
    // Convert a meter conversion factor to its UOM OEM code.
1099
    // If the factor is close to the approximation margin, then
1100
    // require exact equality, otherwise be loose.
1101
1102
0
    const measurement_unit *pu = this->get_uom(dM);
1103
0
    return pu != nullptr ? pu->oemCode : UNITLABEL_UNKNOWN;
1104
0
}
1105
1106
UNITLABEL LevellerDataset::id_to_code(const char *pszUnits) const
1107
0
{
1108
    // Convert a readable UOM to its OEM code.
1109
1110
0
    const measurement_unit *pu = this->get_uom(pszUnits);
1111
0
    return pu != nullptr ? pu->oemCode : UNITLABEL_UNKNOWN;
1112
0
}
1113
1114
const char *LevellerDataset::code_to_id(UNITLABEL code) const
1115
0
{
1116
    // Convert a measurement unit's OEM ID to its readable ID.
1117
1118
0
    const measurement_unit *pu = this->get_uom(code);
1119
0
    return pu != nullptr ? pu->pszID : nullptr;
1120
0
}
1121
1122
const measurement_unit *LevellerDataset::get_uom(const char *pszUnits)
1123
0
{
1124
0
    for (size_t i = 0; i < CPL_ARRAYSIZE(kUnits); i++)
1125
0
    {
1126
0
        if (strcmp(pszUnits, kUnits[i].pszID) == 0)
1127
0
            return &kUnits[i];
1128
0
    }
1129
0
    CPLError(CE_Failure, CPLE_AppDefined, "Unknown measurement units: %s",
1130
0
             pszUnits);
1131
0
    return nullptr;
1132
0
}
1133
1134
const measurement_unit *LevellerDataset::get_uom(UNITLABEL code)
1135
0
{
1136
0
    for (size_t i = 0; i < CPL_ARRAYSIZE(kUnits); i++)
1137
0
    {
1138
0
        if (kUnits[i].oemCode == code)
1139
0
            return &kUnits[i];
1140
0
    }
1141
0
    CPLError(CE_Failure, CPLE_AppDefined, "Unknown measurement unit code: %08x",
1142
0
             code);
1143
0
    return nullptr;
1144
0
}
1145
1146
const measurement_unit *LevellerDataset::get_uom(double dM)
1147
0
{
1148
0
    for (size_t i = kFirstLinearMeasureIdx; i < CPL_ARRAYSIZE(kUnits); i++)
1149
0
    {
1150
0
        if (dM >= 1.0e-4)
1151
0
        {
1152
0
            if (approx_equal(dM, kUnits[i].dScale))
1153
0
                return &kUnits[i];
1154
0
        }
1155
0
        else if (dM == kUnits[i].dScale)
1156
0
            return &kUnits[i];
1157
0
    }
1158
0
    CPLError(CE_Failure, CPLE_AppDefined,
1159
0
             "Unknown measurement conversion factor: %f", dM);
1160
0
    return nullptr;
1161
0
}
1162
1163
/************************************************************************/
1164
/*                          convert_measure()                           */
1165
/************************************************************************/
1166
1167
bool LevellerDataset::convert_measure(double d, double &dResult,
1168
                                      const char *pszSpace)
1169
0
{
1170
    // Convert a measure to meters.
1171
1172
0
    for (size_t i = kFirstLinearMeasureIdx; i < CPL_ARRAYSIZE(kUnits); i++)
1173
0
    {
1174
0
        if (str_equal(pszSpace, kUnits[i].pszID))
1175
0
        {
1176
0
            dResult = d * kUnits[i].dScale;
1177
0
            return true;
1178
0
        }
1179
0
    }
1180
0
    CPLError(CE_Failure, CPLE_FileIO, "Unknown linear measurement unit: '%s'",
1181
0
             pszSpace);
1182
0
    return false;
1183
0
}
1184
1185
bool LevellerDataset::make_local_coordsys(const char *pszName,
1186
                                          const char *pszUnits)
1187
0
{
1188
0
    m_oSRS.SetLocalCS(pszName);
1189
0
    double d;
1190
0
    return (convert_measure(1.0, d, pszUnits) &&
1191
0
            OGRERR_NONE == m_oSRS.SetLinearUnits(pszUnits, d));
1192
0
}
1193
1194
bool LevellerDataset::make_local_coordsys(const char *pszName, UNITLABEL code)
1195
0
{
1196
0
    const char *pszUnitID = code_to_id(code);
1197
0
    return pszUnitID != nullptr && make_local_coordsys(pszName, pszUnitID);
1198
0
}
1199
1200
/************************************************************************/
1201
/*                            load_from_file()                            */
1202
/************************************************************************/
1203
1204
bool LevellerDataset::load_from_file(VSILFILE *file, const char *pszFilename)
1205
0
{
1206
    // get hf dimensions
1207
0
    if (!get(nRasterXSize, file, "hf_w"))
1208
0
    {
1209
0
        CPLError(CE_Failure, CPLE_OpenFailed,
1210
0
                 "Cannot determine heightfield width.");
1211
0
        return false;
1212
0
    }
1213
1214
0
    if (!get(nRasterYSize, file, "hf_b"))
1215
0
    {
1216
0
        CPLError(CE_Failure, CPLE_OpenFailed,
1217
0
                 "Cannot determine heightfield breadth.");
1218
0
        return false;
1219
0
    }
1220
1221
0
    if (nRasterXSize < 2 || nRasterYSize < 2)
1222
0
    {
1223
0
        CPLError(CE_Failure, CPLE_OpenFailed,
1224
0
                 "Heightfield raster dimensions too small.");
1225
0
        return false;
1226
0
    }
1227
1228
    // Record start of pixel data
1229
0
    size_t datalen;
1230
0
    if (!locate_data(m_nDataOffset, datalen, file, "hf_data"))
1231
0
    {
1232
0
        CPLError(CE_Failure, CPLE_OpenFailed, "Cannot locate elevation data.");
1233
0
        return false;
1234
0
    }
1235
1236
    // Sanity check: do we have enough pixels?
1237
0
    if (static_cast<GUIntBig>(datalen) !=
1238
0
        static_cast<GUIntBig>(nRasterXSize) *
1239
0
            static_cast<GUIntBig>(nRasterYSize) * sizeof(float))
1240
0
    {
1241
0
        CPLError(CE_Failure, CPLE_OpenFailed,
1242
0
                 "File does not have enough data.");
1243
0
        return false;
1244
0
    }
1245
1246
    // Defaults for raster coordsys.
1247
0
    m_adfTransform[0] = 0.0;
1248
0
    m_adfTransform[1] = 1.0;
1249
0
    m_adfTransform[2] = 0.0;
1250
0
    m_adfTransform[3] = 0.0;
1251
0
    m_adfTransform[4] = 0.0;
1252
0
    m_adfTransform[5] = 1.0;
1253
1254
0
    m_dElevScale = 1.0;
1255
0
    m_dElevBase = 0.0;
1256
0
    strcpy(m_szElevUnits, "");
1257
1258
0
    if (m_version >= 7)
1259
0
    {
1260
        // Read coordsys info.
1261
0
        int csclass = LEV_COORDSYS_RASTER;
1262
0
        /* (void) */ get(csclass, file, "csclass");
1263
1264
0
        if (csclass != LEV_COORDSYS_RASTER)
1265
0
        {
1266
            // Get projection details and units.
1267
0
            if (csclass == LEV_COORDSYS_LOCAL)
1268
0
            {
1269
0
                UNITLABEL unitcode;
1270
                // char szLocalUnits[8];
1271
0
                int unitcode_int;
1272
0
                if (!get(unitcode_int, file, "coordsys_units"))
1273
0
                    unitcode_int = UNITLABEL_M;
1274
0
                unitcode = static_cast<UNITLABEL>(unitcode_int);
1275
1276
0
                if (!make_local_coordsys("Leveller", unitcode))
1277
0
                {
1278
0
                    CPLError(CE_Failure, CPLE_OpenFailed,
1279
0
                             "Cannot define local coordinate system.");
1280
0
                    return false;
1281
0
                }
1282
0
            }
1283
0
            else if (csclass == LEV_COORDSYS_GEO)
1284
0
            {
1285
0
                char szWKT[1024];
1286
0
                if (!get(szWKT, 1023, file, "coordsys_wkt"))
1287
0
                    return false;
1288
1289
0
                m_oSRS.importFromWkt(szWKT);
1290
0
            }
1291
0
            else
1292
0
            {
1293
0
                CPLError(CE_Failure, CPLE_OpenFailed,
1294
0
                         "Unknown coordinate system type in %s.", pszFilename);
1295
0
                return false;
1296
0
            }
1297
1298
            // Get ground extents.
1299
0
            digital_axis axis_ns, axis_ew;
1300
1301
0
            if (axis_ns.get(*this, file, 0) && axis_ew.get(*this, file, 1))
1302
0
            {
1303
0
                m_adfTransform[0] = axis_ew.origin(nRasterXSize);
1304
0
                m_adfTransform[1] = axis_ew.scaling(nRasterXSize);
1305
0
                m_adfTransform[2] = 0.0;
1306
1307
0
                m_adfTransform[3] = axis_ns.origin(nRasterYSize);
1308
0
                m_adfTransform[4] = 0.0;
1309
0
                m_adfTransform[5] = axis_ns.scaling(nRasterYSize);
1310
0
            }
1311
0
        }
1312
1313
        // Get vertical (elev) coordsys.
1314
0
        int bHasVertCS = FALSE;
1315
0
        if (get(bHasVertCS, file, "coordsys_haselevm") && bHasVertCS)
1316
0
        {
1317
0
            get(m_dElevScale, file, "coordsys_em_scale");
1318
0
            get(m_dElevBase, file, "coordsys_em_base");
1319
0
            UNITLABEL unitcode;
1320
0
            int unitcode_int;
1321
0
            if (get(unitcode_int, file, "coordsys_em_units"))
1322
0
            {
1323
0
                unitcode = static_cast<UNITLABEL>(unitcode_int);
1324
0
                const char *pszUnitID = code_to_id(unitcode);
1325
0
                if (pszUnitID != nullptr)
1326
0
                {
1327
0
                    strncpy(m_szElevUnits, pszUnitID, sizeof(m_szElevUnits));
1328
0
                    m_szElevUnits[sizeof(m_szElevUnits) - 1] = '\0';
1329
0
                }
1330
0
                else
1331
0
                {
1332
0
                    CPLError(CE_Failure, CPLE_OpenFailed,
1333
0
                             "Unknown OEM elevation unit of measure (%d)",
1334
0
                             unitcode);
1335
0
                    return false;
1336
0
                }
1337
0
            }
1338
            // datum and localcs are currently unused.
1339
0
        }
1340
0
    }
1341
0
    else
1342
0
    {
1343
        // Legacy files use world units.
1344
0
        char szWorldUnits[32];
1345
0
        strcpy(szWorldUnits, "m");
1346
1347
0
        double dWorldscale = 1.0;
1348
1349
0
        if (get(dWorldscale, file, "hf_worldspacing"))
1350
0
        {
1351
            // m_bHasWorldscale = true;
1352
0
            if (get(szWorldUnits, sizeof(szWorldUnits) - 1, file,
1353
0
                    "hf_worldspacinglabel"))
1354
0
            {
1355
                // Drop long name, if present.
1356
0
                char *p = strchr(szWorldUnits, ' ');
1357
0
                if (p != nullptr)
1358
0
                    *p = 0;
1359
0
            }
1360
1361
#if 0
1362
          // If the units are something besides m/ft/sft,
1363
          // then convert them to meters.
1364
1365
          if(!str_equal("m", szWorldUnits)
1366
             && !str_equal("ft", szWorldUnits)
1367
             && !str_equal("sft", szWorldUnits))
1368
          {
1369
              dWorldscale = this->convert_measure(dWorldscale, szWorldUnits);
1370
              strcpy(szWorldUnits, "m");
1371
          }
1372
#endif
1373
1374
            // Our extents are such that the origin is at the
1375
            // center of the heightfield.
1376
0
            m_adfTransform[0] = -0.5 * dWorldscale * (nRasterXSize - 1);
1377
0
            m_adfTransform[3] = -0.5 * dWorldscale * (nRasterYSize - 1);
1378
0
            m_adfTransform[1] = dWorldscale;
1379
0
            m_adfTransform[5] = dWorldscale;
1380
0
        }
1381
0
        m_dElevScale = dWorldscale;  // this was 1.0 before because
1382
        // we were converting to real elevs ourselves, but
1383
        // some callers may want both the raw pixels and the
1384
        // transform to get real elevs.
1385
1386
0
        if (!make_local_coordsys("Leveller world space", szWorldUnits))
1387
0
        {
1388
0
            CPLError(CE_Failure, CPLE_OpenFailed,
1389
0
                     "Cannot define local coordinate system.");
1390
0
            return false;
1391
0
        }
1392
0
    }
1393
1394
0
    return true;
1395
0
}
1396
1397
/************************************************************************/
1398
/*                          GetSpatialRef()                             */
1399
/************************************************************************/
1400
1401
const OGRSpatialReference *LevellerDataset::GetSpatialRef() const
1402
0
{
1403
0
    return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
1404
0
}
1405
1406
/************************************************************************/
1407
/*                          GetGeoTransform()                           */
1408
/************************************************************************/
1409
1410
CPLErr LevellerDataset::GetGeoTransform(double *padfTransform)
1411
1412
0
{
1413
0
    memcpy(padfTransform, m_adfTransform, sizeof(m_adfTransform));
1414
0
    return CE_None;
1415
0
}
1416
1417
/************************************************************************/
1418
/*                              Identify()                              */
1419
/************************************************************************/
1420
1421
int LevellerDataset::Identify(GDALOpenInfo *poOpenInfo)
1422
16.2k
{
1423
16.2k
    if (poOpenInfo->nHeaderBytes < 4)
1424
476
        return FALSE;
1425
1426
15.7k
    return STARTS_WITH_CI(
1427
16.2k
        reinterpret_cast<const char *>(poOpenInfo->pabyHeader), "trrn");
1428
16.2k
}
1429
1430
/************************************************************************/
1431
/*                                Open()                                */
1432
/************************************************************************/
1433
1434
GDALDataset *LevellerDataset::Open(GDALOpenInfo *poOpenInfo)
1435
0
{
1436
    // The file should have at least 5 header bytes
1437
    // and hf_w, hf_b, and hf_data tags.
1438
1439
0
    if (poOpenInfo->nHeaderBytes < 5 + 13 + 13 + 16 ||
1440
0
        poOpenInfo->fpL == nullptr)
1441
0
        return nullptr;
1442
1443
0
    if (!LevellerDataset::Identify(poOpenInfo))
1444
0
        return nullptr;
1445
1446
0
    const int version = poOpenInfo->pabyHeader[4];
1447
0
    if (version < 4 || version > 12)
1448
0
        return nullptr;
1449
1450
    /* -------------------------------------------------------------------- */
1451
    /*      Create a corresponding GDALDataset.                             */
1452
    /* -------------------------------------------------------------------- */
1453
1454
0
    LevellerDataset *poDS = new LevellerDataset();
1455
1456
0
    poDS->m_version = version;
1457
1458
0
    poDS->m_fp = poOpenInfo->fpL;
1459
0
    poOpenInfo->fpL = nullptr;
1460
0
    poDS->eAccess = poOpenInfo->eAccess;
1461
1462
    /* -------------------------------------------------------------------- */
1463
    /*      Read the file.                                                  */
1464
    /* -------------------------------------------------------------------- */
1465
0
    if (!poDS->load_from_file(poDS->m_fp, poOpenInfo->pszFilename))
1466
0
    {
1467
0
        delete poDS;
1468
0
        return nullptr;
1469
0
    }
1470
1471
    /* -------------------------------------------------------------------- */
1472
    /*      Create band information objects.                                */
1473
    /* -------------------------------------------------------------------- */
1474
0
    LevellerRasterBand *poBand = new LevellerRasterBand(poDS);
1475
0
    poDS->SetBand(1, poBand);
1476
0
    if (!poBand->Init())
1477
0
    {
1478
0
        delete poDS;
1479
0
        return nullptr;
1480
0
    }
1481
1482
0
    poDS->SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
1483
1484
    /* -------------------------------------------------------------------- */
1485
    /*      Initialize any PAM information.                                 */
1486
    /* -------------------------------------------------------------------- */
1487
0
    poDS->SetDescription(poOpenInfo->pszFilename);
1488
0
    poDS->TryLoadXML();
1489
1490
    /* -------------------------------------------------------------------- */
1491
    /*      Check for external overviews.                                   */
1492
    /* -------------------------------------------------------------------- */
1493
0
    poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename,
1494
0
                                poOpenInfo->GetSiblingFiles());
1495
1496
0
    return poDS;
1497
0
}
1498
1499
/************************************************************************/
1500
/*                        GDALRegister_Leveller()                       */
1501
/************************************************************************/
1502
1503
void GDALRegister_Leveller()
1504
1505
2
{
1506
2
    if (GDALGetDriverByName("Leveller") != nullptr)
1507
0
        return;
1508
1509
2
    GDALDriver *poDriver = new GDALDriver();
1510
1511
2
    poDriver->SetDescription("Leveller");
1512
2
    poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1513
2
    poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "ter");
1514
2
    poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Leveller heightfield");
1515
2
    poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
1516
2
                              "drivers/raster/leveller.html");
1517
2
    poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1518
1519
2
    poDriver->pfnIdentify = LevellerDataset::Identify;
1520
2
    poDriver->pfnOpen = LevellerDataset::Open;
1521
2
    poDriver->pfnCreate = LevellerDataset::Create;
1522
1523
2
    GetGDALDriverManager()->RegisterDriver(poDriver);
1524
2
}