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