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 "ogr_spatialref.h"
21 :
22 1540 : static bool str_equal(const char *_s1, const char *_s2)
23 : {
24 1540 : return 0 == strcmp(_s1, _s2);
25 : }
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 0 : static bool approx_equal(double a, double b)
212 : {
213 0 : const double epsilon = 1e-5;
214 0 : return fabs(a - b) <= epsilon;
215 : }
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 0 : static bool get(size_t &n, VSILFILE *fp, const char *psz)
248 : {
249 0 : return get((int &)n, fp, psz);
250 : }
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 : {
302 0 : m_d[0] = 0.0;
303 0 : m_d[1] = 0.0;
304 0 : }
305 :
306 0 : bool get(LevellerDataset &ds, VSILFILE *fp, int n)
307 : {
308 : 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 : }
323 :
324 0 : double origin(size_t pixels) const
325 : {
326 0 : if (m_fixedEnd == 1)
327 : {
328 0 : switch (m_eStyle)
329 : {
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 : }
336 : }
337 0 : return m_d[0];
338 : }
339 :
340 0 : double scaling(size_t pixels) const
341 : {
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 : }
348 :
349 0 : double length(int pixels) const
350 : {
351 : // Return the signed length of the axis.
352 :
353 0 : switch (m_eStyle)
354 : {
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 : }
364 0 : CPLAssert(false);
365 : return 0.0;
366 : }
367 :
368 : protected:
369 : int m_eStyle;
370 : size_t 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 2 : LevellerRasterBand::LevellerRasterBand(LevellerDataset *poDSIn)
408 2 : : m_pLine(nullptr), m_bFirstTime(true)
409 : {
410 2 : poDS = poDSIn;
411 2 : nBand = 1;
412 :
413 2 : eDataType = GDT_Float32;
414 :
415 2 : nBlockXSize = poDS->GetRasterXSize();
416 2 : nBlockYSize = 1; // poDS->GetRasterYSize();
417 2 : }
418 :
419 : /************************************************************************/
420 : /* Init() */
421 : /************************************************************************/
422 :
423 2 : bool LevellerRasterBand::Init()
424 : {
425 2 : m_pLine = reinterpret_cast<float *>(
426 2 : VSI_MALLOC2_VERBOSE(sizeof(float), nBlockXSize));
427 2 : return m_pLine != nullptr;
428 : }
429 :
430 4 : LevellerRasterBand::~LevellerRasterBand()
431 : {
432 2 : CPLFree(m_pLine);
433 4 : }
434 :
435 : /************************************************************************/
436 : /* IWriteBlock() */
437 : /************************************************************************/
438 :
439 0 : CPLErr LevellerRasterBand::IWriteBlock(CPL_UNUSED int nBlockXOff,
440 : int nBlockYOff, void *pImage)
441 : {
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 : {
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 : }
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 : {
466 0 : for (size_t x = 0; x < (size_t)nBlockXSize; x++)
467 : {
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 : }
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 : }
479 :
480 0 : return CE_Failure;
481 : }
482 :
483 0 : CPLErr LevellerRasterBand::SetUnitType(const char *psz)
484 : {
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 : }
494 :
495 : /************************************************************************/
496 : /* IReadBlock() */
497 : /************************************************************************/
498 :
499 96 : CPLErr LevellerRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
500 : void *pImage)
501 :
502 : {
503 : CPLAssert(sizeof(float) == sizeof(GInt32));
504 96 : CPLAssert(nBlockXOff == 0);
505 96 : CPLAssert(pImage != nullptr);
506 :
507 96 : LevellerDataset *poGDS = reinterpret_cast<LevellerDataset *>(poDS);
508 :
509 : /* -------------------------------------------------------------------- */
510 : /* Seek to scanline. */
511 : /* -------------------------------------------------------------------- */
512 96 : const size_t rowbytes = nBlockXSize * sizeof(float);
513 :
514 96 : if (0 != VSIFSeekL(poGDS->m_fp,
515 96 : poGDS->m_nDataOffset + nBlockYOff * rowbytes, SEEK_SET))
516 : {
517 0 : CPLError(CE_Failure, CPLE_FileIO, "Leveller seek failed: %s",
518 0 : VSIStrerror(errno));
519 0 : return CE_Failure;
520 : }
521 :
522 : /* -------------------------------------------------------------------- */
523 : /* Read the scanline into the image buffer. */
524 : /* -------------------------------------------------------------------- */
525 :
526 96 : if (VSIFReadL(pImage, rowbytes, 1, poGDS->m_fp) != 1)
527 : {
528 0 : CPLError(CE_Failure, CPLE_FileIO, "Leveller read failed: %s",
529 0 : VSIStrerror(errno));
530 0 : return CE_Failure;
531 : }
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 96 : float *pf = reinterpret_cast<float *>(pImage);
544 :
545 96 : if (poGDS->m_version < 6)
546 : {
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 : }
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 96 : return CE_None;
561 : }
562 :
563 : /************************************************************************/
564 : /* GetUnitType() */
565 : /************************************************************************/
566 0 : const char *LevellerRasterBand::GetUnitType()
567 : {
568 : // Return elevation units.
569 :
570 0 : LevellerDataset *poGDS = reinterpret_cast<LevellerDataset *>(poDS);
571 :
572 0 : return poGDS->m_szElevUnits;
573 : }
574 :
575 : /************************************************************************/
576 : /* GetScale() */
577 : /************************************************************************/
578 :
579 0 : double LevellerRasterBand::GetScale(int *pbSuccess)
580 : {
581 0 : LevellerDataset *poGDS = reinterpret_cast<LevellerDataset *>(poDS);
582 0 : if (pbSuccess != nullptr)
583 0 : *pbSuccess = TRUE;
584 0 : return poGDS->m_dElevScale;
585 : }
586 :
587 : /************************************************************************/
588 : /* GetOffset() */
589 : /************************************************************************/
590 :
591 0 : double LevellerRasterBand::GetOffset(int *pbSuccess)
592 : {
593 0 : LevellerDataset *poGDS = reinterpret_cast<LevellerDataset *>(poDS);
594 0 : if (pbSuccess != nullptr)
595 0 : *pbSuccess = TRUE;
596 0 : return poGDS->m_dElevBase;
597 : }
598 :
599 : /************************************************************************/
600 : /* ==================================================================== */
601 : /* LevellerDataset */
602 : /* ==================================================================== */
603 : /************************************************************************/
604 :
605 : /************************************************************************/
606 : /* LevellerDataset() */
607 : /************************************************************************/
608 :
609 6 : LevellerDataset::LevellerDataset()
610 : : m_version(0), m_pszFilename(nullptr), m_dElevScale(), m_dElevBase(),
611 6 : m_fp(nullptr), m_nDataOffset()
612 : {
613 6 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
614 6 : memset(m_szElevUnits, 0, sizeof(m_szElevUnits));
615 6 : memset(m_adfTransform, 0, sizeof(m_adfTransform));
616 6 : memset(m_dLogSpan, 0, sizeof(m_dLogSpan));
617 6 : }
618 :
619 : /************************************************************************/
620 : /* ~LevellerDataset() */
621 : /************************************************************************/
622 :
623 12 : LevellerDataset::~LevellerDataset()
624 : {
625 6 : FlushCache(true);
626 :
627 6 : CPLFree(m_pszFilename);
628 :
629 6 : if (m_fp != nullptr)
630 4 : VSIFCloseL(m_fp);
631 12 : }
632 :
633 0 : static double degrees_to_radians(double d)
634 : {
635 0 : return d * (M_PI / 180);
636 : }
637 :
638 0 : static double average(double a, double b)
639 : {
640 0 : return 0.5 * (a + b);
641 : }
642 :
643 0 : void LevellerDataset::raw_to_proj(double x, double y, double &xp,
644 : double &yp) const
645 : {
646 0 : xp = x * m_adfTransform[1] + m_adfTransform[0];
647 0 : yp = y * m_adfTransform[5] + m_adfTransform[3];
648 0 : }
649 :
650 0 : bool LevellerDataset::compute_elev_scaling(const OGRSpatialReference &sr)
651 : {
652 0 : const char *pszGroundUnits = nullptr;
653 :
654 0 : if (!sr.IsGeographic())
655 : {
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 : }
667 : else
668 : {
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 : 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 : 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 : }
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 : }
708 :
709 0 : bool LevellerDataset::write_header()
710 : {
711 : 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 : {
719 0 : CPLError(CE_Failure, CPLE_FileIO, "Could not write header");
720 0 : return false;
721 : }
722 :
723 0 : m_dElevBase = 0.0;
724 0 : m_dElevScale = 1.0;
725 :
726 0 : if (m_oSRS.IsEmpty())
727 : {
728 0 : write_tag("csclass", LEV_COORDSYS_RASTER);
729 : }
730 : else
731 : {
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 : {
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 : }
756 :
757 0 : if (m_oSRS.IsLocal())
758 : {
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 : }
765 : else
766 : {
767 0 : write_tag("csclass", LEV_COORDSYS_GEO);
768 : }
769 :
770 0 : if (m_adfTransform[2] != 0.0 || m_adfTransform[4] != 0.0)
771 : {
772 0 : CPLError(CE_Failure, CPLE_IllegalArg,
773 : "Cannot handle rotated geotransform");
774 0 : return false;
775 : }
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 : }
793 :
794 0 : this->write_tag_start("hf_data",
795 0 : sizeof(float) * nRasterXSize * nRasterYSize);
796 :
797 0 : return true;
798 : }
799 :
800 : /************************************************************************/
801 : /* SetGeoTransform() */
802 : /************************************************************************/
803 :
804 0 : CPLErr LevellerDataset::SetGeoTransform(double *padfGeoTransform)
805 : {
806 0 : memcpy(m_adfTransform, padfGeoTransform, sizeof(m_adfTransform));
807 :
808 0 : return CE_None;
809 : }
810 :
811 : /************************************************************************/
812 : /* SetSpatialRef() */
813 : /************************************************************************/
814 :
815 0 : CPLErr LevellerDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
816 : {
817 0 : m_oSRS.Clear();
818 0 : if (poSRS)
819 0 : m_oSRS = *poSRS;
820 :
821 0 : return CE_None;
822 : }
823 :
824 : /************************************************************************/
825 : /* Create() */
826 : /************************************************************************/
827 49 : GDALDataset *LevellerDataset::Create(const char *pszFilename, int nXSize,
828 : int nYSize, int nBandsIn,
829 : GDALDataType eType, char **papszOptions)
830 : {
831 49 : if (nBandsIn != 1)
832 : {
833 22 : CPLError(CE_Failure, CPLE_IllegalArg, "Band count must be 1");
834 22 : return nullptr;
835 : }
836 :
837 27 : if (eType != GDT_Float32)
838 : {
839 23 : CPLError(CE_Failure, CPLE_IllegalArg, "Pixel type must be Float32");
840 23 : return nullptr;
841 : }
842 :
843 4 : if (nXSize < 2 || nYSize < 2)
844 : {
845 0 : CPLError(CE_Failure, CPLE_IllegalArg,
846 : "One or more raster dimensions too small");
847 0 : return nullptr;
848 : }
849 :
850 4 : LevellerDataset *poDS = new LevellerDataset;
851 :
852 4 : poDS->eAccess = GA_Update;
853 :
854 4 : poDS->m_pszFilename = CPLStrdup(pszFilename);
855 :
856 4 : poDS->m_fp = VSIFOpenL(pszFilename, "wb+");
857 :
858 4 : if (poDS->m_fp == nullptr)
859 : {
860 2 : CPLError(CE_Failure, CPLE_OpenFailed,
861 : "Attempt to create file `%s' failed.", pszFilename);
862 2 : delete poDS;
863 2 : return nullptr;
864 : }
865 :
866 : // Header will be written the first time IWriteBlock
867 : // is called.
868 :
869 2 : poDS->nRasterXSize = nXSize;
870 2 : poDS->nRasterYSize = nYSize;
871 :
872 2 : const char *pszValue = CSLFetchNameValue(papszOptions, "MINUSERPIXELVALUE");
873 2 : if (pszValue != nullptr)
874 0 : poDS->m_dLogSpan[0] = CPLAtof(pszValue);
875 : else
876 : {
877 2 : delete poDS;
878 2 : CPLError(CE_Failure, CPLE_IllegalArg,
879 : "MINUSERPIXELVALUE must be specified.");
880 2 : return nullptr;
881 : }
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 : {
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 : }
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 : {
902 0 : delete poDS;
903 0 : return nullptr;
904 : }
905 :
906 0 : return poDS;
907 : }
908 :
909 0 : bool LevellerDataset::write_byte(size_t n)
910 : {
911 0 : unsigned char uch = static_cast<unsigned char>(n);
912 0 : return 1 == VSIFWriteL(&uch, 1, 1, m_fp);
913 : }
914 :
915 0 : bool LevellerDataset::write(int n)
916 : {
917 0 : CPL_LSBPTR32(&n);
918 0 : return 1 == VSIFWriteL(&n, sizeof(n), 1, m_fp);
919 : }
920 :
921 0 : bool LevellerDataset::write(size_t n)
922 : {
923 0 : GUInt32 n32 = (GUInt32)n;
924 0 : CPL_LSBPTR32(&n32);
925 0 : return 1 == VSIFWriteL(&n32, sizeof(n32), 1, m_fp);
926 : }
927 :
928 0 : bool LevellerDataset::write(double d)
929 : {
930 0 : CPL_LSBPTR64(&d);
931 0 : return 1 == VSIFWriteL(&d, sizeof(d), 1, m_fp);
932 : }
933 :
934 0 : bool LevellerDataset::write_tag_start(const char *pszTag, size_t n)
935 : {
936 0 : if (this->write_byte(strlen(pszTag)))
937 : {
938 0 : return (1 == VSIFWriteL(pszTag, strlen(pszTag), 1, m_fp) &&
939 0 : this->write(n));
940 : }
941 :
942 0 : return false;
943 : }
944 :
945 0 : bool LevellerDataset::write_tag(const char *pszTag, int n)
946 : {
947 0 : return (this->write_tag_start(pszTag, sizeof(n)) && this->write(n));
948 : }
949 :
950 0 : bool LevellerDataset::write_tag(const char *pszTag, size_t n)
951 : {
952 0 : return (this->write_tag_start(pszTag, sizeof(n)) && this->write(n));
953 : }
954 :
955 0 : bool LevellerDataset::write_tag(const char *pszTag, double d)
956 : {
957 0 : return (this->write_tag_start(pszTag, sizeof(d)) && this->write(d));
958 : }
959 :
960 0 : bool LevellerDataset::write_tag(const char *pszTag, const char *psz)
961 : {
962 0 : CPLAssert(strlen(pszTag) <= kMaxTagNameLen);
963 :
964 : 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 : {
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 : }
974 0 : return false;
975 : }
976 :
977 10 : bool LevellerDataset::locate_data(vsi_l_offset &offset, size_t &len,
978 : VSILFILE *fp, const char *pszTag)
979 : {
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 10 : if (0 != VSIFSeekL(fp, 5, SEEK_SET))
986 0 : return false;
987 :
988 10 : const int kMaxDescLen = 64;
989 : for (;;)
990 : {
991 : unsigned char c;
992 1498 : if (1 != VSIFReadL(&c, sizeof(c), 1, fp))
993 10 : return false;
994 :
995 1498 : const size_t descriptorLen = c;
996 1498 : if (descriptorLen == 0 || descriptorLen > (size_t)kMaxDescLen)
997 0 : return false;
998 :
999 : char descriptor[kMaxDescLen + 1];
1000 1498 : if (1 != VSIFReadL(descriptor, descriptorLen, 1, fp))
1001 0 : return false;
1002 :
1003 : GUInt32 datalen;
1004 1498 : if (1 != VSIFReadL(&datalen, sizeof(datalen), 1, fp))
1005 0 : return false;
1006 :
1007 1498 : CPL_LSBPTR32(&datalen);
1008 1498 : descriptor[descriptorLen] = 0;
1009 1498 : if (str_equal(descriptor, pszTag))
1010 : {
1011 10 : len = (size_t)datalen;
1012 10 : offset = VSIFTellL(fp);
1013 10 : return true;
1014 : }
1015 : else
1016 : {
1017 : // Seek to next tag.
1018 1488 : if (0 != VSIFSeekL(fp, (vsi_l_offset)datalen, SEEK_CUR))
1019 0 : return false;
1020 : }
1021 1488 : }
1022 : }
1023 :
1024 : /************************************************************************/
1025 : /* get() */
1026 : /************************************************************************/
1027 :
1028 4 : bool LevellerDataset::get(int &n, VSILFILE *fp, const char *psz)
1029 : {
1030 : vsi_l_offset offset;
1031 : size_t len;
1032 :
1033 4 : if (locate_data(offset, len, fp, psz))
1034 : {
1035 : GInt32 value;
1036 4 : if (1 == VSIFReadL(&value, sizeof(value), 1, fp))
1037 : {
1038 4 : CPL_LSBPTR32(&value);
1039 4 : n = static_cast<int>(value);
1040 4 : return true;
1041 : }
1042 : }
1043 0 : return false;
1044 : }
1045 :
1046 : /************************************************************************/
1047 : /* get() */
1048 : /************************************************************************/
1049 :
1050 2 : bool LevellerDataset::get(double &d, VSILFILE *fp, const char *pszTag)
1051 : {
1052 : vsi_l_offset offset;
1053 : size_t len;
1054 :
1055 2 : if (locate_data(offset, len, fp, pszTag))
1056 : {
1057 2 : if (1 == VSIFReadL(&d, sizeof(d), 1, fp))
1058 : {
1059 2 : CPL_LSBPTR64(&d);
1060 2 : return true;
1061 : }
1062 : }
1063 0 : return false;
1064 : }
1065 :
1066 : /************************************************************************/
1067 : /* get() */
1068 : /************************************************************************/
1069 2 : bool LevellerDataset::get(char *pszValue, size_t maxchars, VSILFILE *fp,
1070 : const char *pszTag)
1071 : {
1072 : char szTag[65];
1073 :
1074 : // We can assume 8-bit encoding, so just go straight
1075 : // to the *_d tag.
1076 2 : snprintf(szTag, sizeof(szTag), "%sd", pszTag);
1077 :
1078 : vsi_l_offset offset;
1079 : size_t len;
1080 :
1081 2 : if (locate_data(offset, len, fp, szTag))
1082 : {
1083 2 : if (len > maxchars)
1084 0 : return false;
1085 :
1086 2 : if (1 == VSIFReadL(pszValue, len, 1, fp))
1087 : {
1088 2 : pszValue[len] = 0; // terminate C-string
1089 2 : return true;
1090 : }
1091 : }
1092 :
1093 0 : return false;
1094 : }
1095 :
1096 0 : UNITLABEL LevellerDataset::meter_measure_to_code(double dM) const
1097 : {
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 : }
1105 :
1106 0 : UNITLABEL LevellerDataset::id_to_code(const char *pszUnits) const
1107 : {
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 : }
1113 :
1114 0 : const char *LevellerDataset::code_to_id(UNITLABEL code) const
1115 : {
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 : }
1121 :
1122 0 : const measurement_unit *LevellerDataset::get_uom(const char *pszUnits)
1123 : {
1124 0 : for (size_t i = 0; i < CPL_ARRAYSIZE(kUnits); i++)
1125 : {
1126 0 : if (strcmp(pszUnits, kUnits[i].pszID) == 0)
1127 0 : return &kUnits[i];
1128 : }
1129 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unknown measurement units: %s",
1130 : pszUnits);
1131 0 : return nullptr;
1132 : }
1133 :
1134 0 : const measurement_unit *LevellerDataset::get_uom(UNITLABEL code)
1135 : {
1136 0 : for (size_t i = 0; i < CPL_ARRAYSIZE(kUnits); i++)
1137 : {
1138 0 : if (kUnits[i].oemCode == code)
1139 0 : return &kUnits[i];
1140 : }
1141 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unknown measurement unit code: %08x",
1142 : code);
1143 0 : return nullptr;
1144 : }
1145 :
1146 0 : const measurement_unit *LevellerDataset::get_uom(double dM)
1147 : {
1148 0 : for (size_t i = kFirstLinearMeasureIdx; i < CPL_ARRAYSIZE(kUnits); i++)
1149 : {
1150 0 : if (dM >= 1.0e-4)
1151 : {
1152 0 : if (approx_equal(dM, kUnits[i].dScale))
1153 0 : return &kUnits[i];
1154 : }
1155 0 : else if (dM == kUnits[i].dScale)
1156 0 : return &kUnits[i];
1157 : }
1158 0 : CPLError(CE_Failure, CPLE_AppDefined,
1159 : "Unknown measurement conversion factor: %f", dM);
1160 0 : return nullptr;
1161 : }
1162 :
1163 : /************************************************************************/
1164 : /* convert_measure() */
1165 : /************************************************************************/
1166 :
1167 2 : bool LevellerDataset::convert_measure(double d, double &dResult,
1168 : const char *pszSpace)
1169 : {
1170 : // Convert a measure to meters.
1171 :
1172 42 : for (size_t i = kFirstLinearMeasureIdx; i < CPL_ARRAYSIZE(kUnits); i++)
1173 : {
1174 42 : if (str_equal(pszSpace, kUnits[i].pszID))
1175 : {
1176 2 : dResult = d * kUnits[i].dScale;
1177 2 : return true;
1178 : }
1179 : }
1180 0 : CPLError(CE_Failure, CPLE_FileIO, "Unknown linear measurement unit: '%s'",
1181 : pszSpace);
1182 0 : return false;
1183 : }
1184 :
1185 2 : bool LevellerDataset::make_local_coordsys(const char *pszName,
1186 : const char *pszUnits)
1187 : {
1188 2 : m_oSRS.SetLocalCS(pszName);
1189 : double d;
1190 4 : return (convert_measure(1.0, d, pszUnits) &&
1191 4 : OGRERR_NONE == m_oSRS.SetLinearUnits(pszUnits, d));
1192 : }
1193 :
1194 0 : bool LevellerDataset::make_local_coordsys(const char *pszName, UNITLABEL code)
1195 : {
1196 0 : const char *pszUnitID = code_to_id(code);
1197 0 : return pszUnitID != nullptr && make_local_coordsys(pszName, pszUnitID);
1198 : }
1199 :
1200 : /************************************************************************/
1201 : /* load_from_file() */
1202 : /************************************************************************/
1203 :
1204 2 : bool LevellerDataset::load_from_file(VSILFILE *file, const char *pszFilename)
1205 : {
1206 : // get hf dimensions
1207 2 : if (!get(nRasterXSize, file, "hf_w"))
1208 : {
1209 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1210 : "Cannot determine heightfield width.");
1211 0 : return false;
1212 : }
1213 :
1214 2 : if (!get(nRasterYSize, file, "hf_b"))
1215 : {
1216 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1217 : "Cannot determine heightfield breadth.");
1218 0 : return false;
1219 : }
1220 :
1221 2 : if (nRasterXSize < 2 || nRasterYSize < 2)
1222 : {
1223 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1224 : "Heightfield raster dimensions too small.");
1225 0 : return false;
1226 : }
1227 :
1228 : // Record start of pixel data
1229 : size_t datalen;
1230 2 : if (!locate_data(m_nDataOffset, datalen, file, "hf_data"))
1231 : {
1232 0 : CPLError(CE_Failure, CPLE_OpenFailed, "Cannot locate elevation data.");
1233 0 : return false;
1234 : }
1235 :
1236 : // Sanity check: do we have enough pixels?
1237 2 : if (static_cast<GUIntBig>(datalen) !=
1238 2 : static_cast<GUIntBig>(nRasterXSize) *
1239 2 : static_cast<GUIntBig>(nRasterYSize) * sizeof(float))
1240 : {
1241 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1242 : "File does not have enough data.");
1243 0 : return false;
1244 : }
1245 :
1246 : // Defaults for raster coordsys.
1247 2 : m_adfTransform[0] = 0.0;
1248 2 : m_adfTransform[1] = 1.0;
1249 2 : m_adfTransform[2] = 0.0;
1250 2 : m_adfTransform[3] = 0.0;
1251 2 : m_adfTransform[4] = 0.0;
1252 2 : m_adfTransform[5] = 1.0;
1253 :
1254 2 : m_dElevScale = 1.0;
1255 2 : m_dElevBase = 0.0;
1256 2 : strcpy(m_szElevUnits, "");
1257 :
1258 2 : if (m_version >= 7)
1259 : {
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 : {
1266 : // Get projection details and units.
1267 0 : if (csclass == LEV_COORDSYS_LOCAL)
1268 : {
1269 : UNITLABEL unitcode;
1270 : // char szLocalUnits[8];
1271 : 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 : {
1278 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1279 : "Cannot define local coordinate system.");
1280 0 : return false;
1281 : }
1282 : }
1283 0 : else if (csclass == LEV_COORDSYS_GEO)
1284 : {
1285 : char szWKT[1024];
1286 0 : if (!get(szWKT, 1023, file, "coordsys_wkt"))
1287 0 : return false;
1288 :
1289 0 : m_oSRS.importFromWkt(szWKT);
1290 : }
1291 : else
1292 : {
1293 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1294 : "Unknown coordinate system type in %s.", pszFilename);
1295 0 : return false;
1296 : }
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 : {
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 : }
1311 : }
1312 :
1313 : // Get vertical (elev) coordsys.
1314 0 : int bHasVertCS = FALSE;
1315 0 : if (get(bHasVertCS, file, "coordsys_haselevm") && bHasVertCS)
1316 : {
1317 0 : get(m_dElevScale, file, "coordsys_em_scale");
1318 0 : get(m_dElevBase, file, "coordsys_em_base");
1319 : UNITLABEL unitcode;
1320 : int unitcode_int;
1321 0 : if (get(unitcode_int, file, "coordsys_em_units"))
1322 : {
1323 0 : unitcode = static_cast<UNITLABEL>(unitcode_int);
1324 0 : const char *pszUnitID = code_to_id(unitcode);
1325 0 : if (pszUnitID != nullptr)
1326 : {
1327 0 : strncpy(m_szElevUnits, pszUnitID, sizeof(m_szElevUnits));
1328 0 : m_szElevUnits[sizeof(m_szElevUnits) - 1] = '\0';
1329 : }
1330 : else
1331 : {
1332 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1333 : "Unknown OEM elevation unit of measure (%d)",
1334 : unitcode);
1335 0 : return false;
1336 : }
1337 : }
1338 : // datum and localcs are currently unused.
1339 : }
1340 : }
1341 : else
1342 : {
1343 : // Legacy files use world units.
1344 : char szWorldUnits[32];
1345 2 : strcpy(szWorldUnits, "m");
1346 :
1347 2 : double dWorldscale = 1.0;
1348 :
1349 2 : if (get(dWorldscale, file, "hf_worldspacing"))
1350 : {
1351 : // m_bHasWorldscale = true;
1352 2 : if (get(szWorldUnits, sizeof(szWorldUnits) - 1, file,
1353 : "hf_worldspacinglabel"))
1354 : {
1355 : // Drop long name, if present.
1356 2 : char *p = strchr(szWorldUnits, ' ');
1357 2 : if (p != nullptr)
1358 2 : *p = 0;
1359 : }
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 2 : m_adfTransform[0] = -0.5 * dWorldscale * (nRasterXSize - 1);
1377 2 : m_adfTransform[3] = -0.5 * dWorldscale * (nRasterYSize - 1);
1378 2 : m_adfTransform[1] = dWorldscale;
1379 2 : m_adfTransform[5] = dWorldscale;
1380 : }
1381 2 : 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 2 : if (!make_local_coordsys("Leveller world space", szWorldUnits))
1387 : {
1388 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1389 : "Cannot define local coordinate system.");
1390 0 : return false;
1391 : }
1392 : }
1393 :
1394 2 : return true;
1395 : }
1396 :
1397 : /************************************************************************/
1398 : /* GetSpatialRef() */
1399 : /************************************************************************/
1400 :
1401 0 : const OGRSpatialReference *LevellerDataset::GetSpatialRef() const
1402 : {
1403 0 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
1404 : }
1405 :
1406 : /************************************************************************/
1407 : /* GetGeoTransform() */
1408 : /************************************************************************/
1409 :
1410 0 : CPLErr LevellerDataset::GetGeoTransform(double *padfTransform)
1411 :
1412 : {
1413 0 : memcpy(padfTransform, m_adfTransform, sizeof(m_adfTransform));
1414 0 : return CE_None;
1415 : }
1416 :
1417 : /************************************************************************/
1418 : /* Identify() */
1419 : /************************************************************************/
1420 :
1421 55407 : int LevellerDataset::Identify(GDALOpenInfo *poOpenInfo)
1422 : {
1423 55407 : if (poOpenInfo->nHeaderBytes < 4)
1424 49196 : return FALSE;
1425 :
1426 6211 : return STARTS_WITH_CI(
1427 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader), "trrn");
1428 : }
1429 :
1430 : /************************************************************************/
1431 : /* Open() */
1432 : /************************************************************************/
1433 :
1434 2 : GDALDataset *LevellerDataset::Open(GDALOpenInfo *poOpenInfo)
1435 : {
1436 : // The file should have at least 5 header bytes
1437 : // and hf_w, hf_b, and hf_data tags.
1438 :
1439 2 : if (poOpenInfo->nHeaderBytes < 5 + 13 + 13 + 16 ||
1440 2 : poOpenInfo->fpL == nullptr)
1441 0 : return nullptr;
1442 :
1443 2 : if (!LevellerDataset::Identify(poOpenInfo))
1444 0 : return nullptr;
1445 :
1446 2 : const int version = poOpenInfo->pabyHeader[4];
1447 2 : if (version < 4 || version > 9)
1448 0 : return nullptr;
1449 :
1450 : /* -------------------------------------------------------------------- */
1451 : /* Create a corresponding GDALDataset. */
1452 : /* -------------------------------------------------------------------- */
1453 :
1454 2 : LevellerDataset *poDS = new LevellerDataset();
1455 :
1456 2 : poDS->m_version = version;
1457 :
1458 2 : poDS->m_fp = poOpenInfo->fpL;
1459 2 : poOpenInfo->fpL = nullptr;
1460 2 : poDS->eAccess = poOpenInfo->eAccess;
1461 :
1462 : /* -------------------------------------------------------------------- */
1463 : /* Read the file. */
1464 : /* -------------------------------------------------------------------- */
1465 2 : if (!poDS->load_from_file(poDS->m_fp, poOpenInfo->pszFilename))
1466 : {
1467 0 : delete poDS;
1468 0 : return nullptr;
1469 : }
1470 :
1471 : /* -------------------------------------------------------------------- */
1472 : /* Create band information objects. */
1473 : /* -------------------------------------------------------------------- */
1474 2 : LevellerRasterBand *poBand = new LevellerRasterBand(poDS);
1475 2 : poDS->SetBand(1, poBand);
1476 2 : if (!poBand->Init())
1477 : {
1478 0 : delete poDS;
1479 0 : return nullptr;
1480 : }
1481 :
1482 2 : poDS->SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
1483 :
1484 : /* -------------------------------------------------------------------- */
1485 : /* Initialize any PAM information. */
1486 : /* -------------------------------------------------------------------- */
1487 2 : poDS->SetDescription(poOpenInfo->pszFilename);
1488 2 : poDS->TryLoadXML();
1489 :
1490 : /* -------------------------------------------------------------------- */
1491 : /* Check for external overviews. */
1492 : /* -------------------------------------------------------------------- */
1493 4 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename,
1494 2 : poOpenInfo->GetSiblingFiles());
1495 :
1496 2 : return poDS;
1497 : }
1498 :
1499 : /************************************************************************/
1500 : /* GDALRegister_Leveller() */
1501 : /************************************************************************/
1502 :
1503 1595 : void GDALRegister_Leveller()
1504 :
1505 : {
1506 1595 : if (GDALGetDriverByName("Leveller") != nullptr)
1507 302 : return;
1508 :
1509 1293 : GDALDriver *poDriver = new GDALDriver();
1510 :
1511 1293 : poDriver->SetDescription("Leveller");
1512 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1513 1293 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "ter");
1514 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Leveller heightfield");
1515 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
1516 1293 : "drivers/raster/leveller.html");
1517 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1518 :
1519 1293 : poDriver->pfnIdentify = LevellerDataset::Identify;
1520 1293 : poDriver->pfnOpen = LevellerDataset::Open;
1521 1293 : poDriver->pfnCreate = LevellerDataset::Create;
1522 :
1523 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1524 : }
|
