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