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