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