1 /******************************************************************************
2 * $Id: btdataset.cpp 29284 2015-06-03 13:26:10Z rouault $
3 *
4 * Project: VTP .bt Driver
5 * Purpose: Implementation of VTP .bt elevation format read/write support.
6 * http://www.vterrain.org/Implementation/Formats/BT.html
7 * Author: Frank Warmerdam, warmerdam@pobox.com
8 *
9 ******************************************************************************
10 * Copyright (c) 2003, Frank Warmerdam
11 * Copyright (c) 2007-2011, Even Rouault <even dot rouault at mines-paris dot org>
12 *
13 * Permission is hereby granted, free of charge, to any person obtaining a
14 * copy of this software and associated documentation files (the "Software"),
15 * to deal in the Software without restriction, including without limitation
16 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 * and/or sell copies of the Software, and to permit persons to whom the
18 * Software is furnished to do so, subject to the following conditions:
19 *
20 * The above copyright notice and this permission notice shall be included
21 * in all copies or substantial portions of the Software.
22 *
23 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 * DEALINGS IN THE SOFTWARE.
30 ****************************************************************************/
31
32 #include "rawdataset.h"
33 #include "ogr_spatialref.h"
34
35 CPL_CVSID("$Id: btdataset.cpp 29284 2015-06-03 13:26:10Z rouault $");
36
37 CPL_C_START
38 void GDALRegister_BT(void);
39 CPL_C_END
40
41 /************************************************************************/
42 /* ==================================================================== */
43 /* BTDataset */
44 /* ==================================================================== */
45 /************************************************************************/
46
47 class BTDataset : public GDALPamDataset
48 {
49 friend class BTRasterBand;
50
51 VSILFILE *fpImage; // image data file.
52
53 int bGeoTransformValid;
54 double adfGeoTransform[6];
55
56 char *pszProjection;
57
58 int nVersionCode; // version times 10.
59
60 int bHeaderModified;
61 unsigned char abyHeader[256];
62
63 float m_fVscale;
64
65 public:
66
67 BTDataset();
68 ~BTDataset();
69
70 virtual const char *GetProjectionRef(void);
71 virtual CPLErr SetProjection( const char * );
72 virtual CPLErr GetGeoTransform( double * );
73 virtual CPLErr SetGeoTransform( double * );
74
75 virtual void FlushCache();
76
77 static GDALDataset *Open( GDALOpenInfo * );
78 static GDALDataset *Create( const char * pszFilename,
79 int nXSize, int nYSize, int nBands,
80 GDALDataType eType, char ** papszOptions );
81 };
82
83 /************************************************************************/
84 /* ==================================================================== */
85 /* BTRasterBand */
86 /* ==================================================================== */
87 /************************************************************************/
88
89 class BTRasterBand : public GDALPamRasterBand
90 {
91 VSILFILE *fpImage;
92
93 public:
94
95 BTRasterBand( GDALDataset * poDS, VSILFILE * fp,
96 GDALDataType eType );
97
98 virtual CPLErr IReadBlock( int, int, void * );
99 virtual CPLErr IWriteBlock( int, int, void * );
100
101 virtual const char* GetUnitType();
102 virtual CPLErr SetUnitType(const char*);
103 virtual double GetNoDataValue( int* = NULL );
104 virtual CPLErr SetNoDataValue( double );
105 };
106
107
108 /************************************************************************/
109 /* BTRasterBand() */
110 /************************************************************************/
111
BTRasterBand(GDALDataset * poDS,VSILFILE * fp,GDALDataType eType)112 BTRasterBand::BTRasterBand( GDALDataset *poDS, VSILFILE *fp, GDALDataType eType )
113
114 {
115 this->poDS = poDS;
116 this->nBand = 1;
117 this->eDataType = eType;
118 this->fpImage = fp;
119
120 nBlockXSize = 1;
121 nBlockYSize = poDS->GetRasterYSize();
122 }
123
124 /************************************************************************/
125 /* IReadBlock() */
126 /************************************************************************/
127
IReadBlock(int nBlockXOff,CPL_UNUSED int nBlockYOff,void * pImage)128 CPLErr BTRasterBand::IReadBlock( int nBlockXOff,
129 CPL_UNUSED int nBlockYOff,
130 void * pImage )
131 {
132 int nDataSize = GDALGetDataTypeSize( eDataType ) / 8;
133 int i;
134
135 CPLAssert( nBlockYOff == 0 );
136
137 /* -------------------------------------------------------------------- */
138 /* Seek to profile. */
139 /* -------------------------------------------------------------------- */
140 if( VSIFSeekL( fpImage,
141 256 + nBlockXOff * nDataSize * (vsi_l_offset) nRasterYSize,
142 SEEK_SET ) != 0 )
143 {
144 CPLError( CE_Failure, CPLE_FileIO,
145 ".bt Seek failed:%s", VSIStrerror( errno ) );
146 return CE_Failure;
147 }
148
149 /* -------------------------------------------------------------------- */
150 /* Read the profile. */
151 /* -------------------------------------------------------------------- */
152 if( VSIFReadL( pImage, nDataSize, nRasterYSize, fpImage ) !=
153 (size_t) nRasterYSize )
154 {
155 CPLError( CE_Failure, CPLE_FileIO,
156 ".bt Read failed:%s", VSIStrerror( errno ) );
157 return CE_Failure;
158 }
159
160 /* -------------------------------------------------------------------- */
161 /* Swap on MSB platforms. */
162 /* -------------------------------------------------------------------- */
163 #ifdef CPL_MSB
164 GDALSwapWords( pImage, nDataSize, nRasterYSize, nDataSize );
165 #endif
166
167 /* -------------------------------------------------------------------- */
168 /* Vertical flip, since GDAL expects values from top to bottom, */
169 /* but in .bt they are bottom to top. */
170 /* -------------------------------------------------------------------- */
171 for( i = 0; i < nRasterYSize / 2; i++ )
172 {
173 GByte abyWrk[8];
174
175 memcpy( abyWrk, ((GByte *) pImage) + i * nDataSize, nDataSize );
176 memcpy( ((GByte *) pImage) + i * nDataSize,
177 ((GByte *) pImage) + (nRasterYSize - i - 1) * nDataSize,
178 nDataSize );
179 memcpy( ((GByte *) pImage) + (nRasterYSize - i - 1) * nDataSize,
180 abyWrk, nDataSize );
181 }
182
183 return CE_None;
184 }
185
186 /************************************************************************/
187 /* IWriteBlock() */
188 /************************************************************************/
189
IWriteBlock(int nBlockXOff,CPL_UNUSED int nBlockYOff,void * pImage)190 CPLErr BTRasterBand::IWriteBlock( int nBlockXOff,
191 CPL_UNUSED int nBlockYOff,
192 void * pImage )
193
194 {
195 int nDataSize = GDALGetDataTypeSize( eDataType ) / 8;
196 GByte *pabyWrkBlock;
197 int i;
198
199 CPLAssert( nBlockYOff == 0 );
200
201 /* -------------------------------------------------------------------- */
202 /* Seek to profile. */
203 /* -------------------------------------------------------------------- */
204 if( VSIFSeekL( fpImage,
205 256 + nBlockXOff * nDataSize * nRasterYSize,
206 SEEK_SET ) != 0 )
207 {
208 CPLError( CE_Failure, CPLE_FileIO,
209 ".bt Seek failed:%s", VSIStrerror( errno ) );
210 return CE_Failure;
211 }
212
213 /* -------------------------------------------------------------------- */
214 /* Allocate working buffer. */
215 /* -------------------------------------------------------------------- */
216 pabyWrkBlock = (GByte *) CPLMalloc(nDataSize * nRasterYSize);
217
218 /* -------------------------------------------------------------------- */
219 /* Vertical flip data into work buffer, since GDAL expects */
220 /* values from top to bottom, but in .bt they are bottom to */
221 /* top. */
222 /* -------------------------------------------------------------------- */
223 for( i = 0; i < nRasterYSize; i++ )
224 {
225 memcpy( pabyWrkBlock + (nRasterYSize - i - 1) * nDataSize,
226 ((GByte *) pImage) + i * nDataSize, nDataSize );
227 }
228
229 /* -------------------------------------------------------------------- */
230 /* Swap on MSB platforms. */
231 /* -------------------------------------------------------------------- */
232 #ifdef CPL_MSB
233 GDALSwapWords( pabyWrkBlock, nDataSize, nRasterYSize, nDataSize );
234 #endif
235
236 /* -------------------------------------------------------------------- */
237 /* Read the profile. */
238 /* -------------------------------------------------------------------- */
239 if( VSIFWriteL( pabyWrkBlock, nDataSize, nRasterYSize, fpImage ) !=
240 (size_t) nRasterYSize )
241 {
242 CPLFree( pabyWrkBlock );
243 CPLError( CE_Failure, CPLE_FileIO,
244 ".bt Write failed:%s", VSIStrerror( errno ) );
245 return CE_Failure;
246 }
247
248 CPLFree( pabyWrkBlock );
249
250 return CE_None;
251 }
252
253
GetNoDataValue(int * pbSuccess)254 double BTRasterBand::GetNoDataValue( int* pbSuccess /*= NULL */ )
255 {
256 if(pbSuccess != NULL)
257 *pbSuccess = TRUE;
258
259 return -32768;
260 }
261
SetNoDataValue(double)262 CPLErr BTRasterBand::SetNoDataValue( double )
263
264 {
265 return CE_None;
266 }
267
268 /************************************************************************/
269 /* GetUnitType() */
270 /************************************************************************/
271
approx_equals(float a,float b)272 static bool approx_equals(float a, float b)
273 {
274 const float epsilon = (float)1e-5;
275 return (fabs(a-b) <= epsilon);
276 }
277
GetUnitType(void)278 const char* BTRasterBand::GetUnitType(void)
279 {
280 const BTDataset& ds = *(BTDataset*)poDS;
281 float f = ds.m_fVscale;
282 if(f == 1.0f)
283 return "m";
284 if(approx_equals(f, 0.3048f))
285 return "ft";
286 if(approx_equals(f, 1200.0f/3937.0f))
287 return "sft";
288
289 // todo: the BT spec allows for any value for
290 // ds.m_fVscale, so rigorous adherence would
291 // require testing for all possible units you
292 // may want to support, such as km, yards, miles, etc.
293 // But m/ft/sft seem to be the top three.
294
295 return "";
296 }
297
298 /************************************************************************/
299 /* SetUnitType() */
300 /************************************************************************/
301
SetUnitType(const char * psz)302 CPLErr BTRasterBand::SetUnitType(const char* psz)
303 {
304 BTDataset& ds = *(BTDataset*)poDS;
305 if(EQUAL(psz, "m"))
306 ds.m_fVscale = 1.0f;
307 else if(EQUAL(psz, "ft"))
308 ds.m_fVscale = 0.3048f;
309 else if(EQUAL(psz, "sft"))
310 ds.m_fVscale = 1200.0f / 3937.0f;
311 else
312 return CE_Failure;
313
314
315 float fScale = ds.m_fVscale;
316
317 CPL_LSBPTR32(&fScale);
318
319 // Update header's elevation scale field.
320 memcpy(ds.abyHeader + 62, &fScale, sizeof(fScale));
321
322 ds.bHeaderModified = TRUE;
323 return CE_None;
324 }
325
326 /************************************************************************/
327 /* ==================================================================== */
328 /* BTDataset */
329 /* ==================================================================== */
330 /************************************************************************/
331
332 /************************************************************************/
333 /* BTDataset() */
334 /************************************************************************/
335
BTDataset()336 BTDataset::BTDataset()
337 {
338 fpImage = NULL;
339 bGeoTransformValid = FALSE;
340 adfGeoTransform[0] = 0.0;
341 adfGeoTransform[1] = 1.0;
342 adfGeoTransform[2] = 0.0;
343 adfGeoTransform[3] = 0.0;
344 adfGeoTransform[4] = 0.0;
345 adfGeoTransform[5] = 1.0;
346 pszProjection = NULL;
347
348 bHeaderModified = FALSE;
349 }
350
351 /************************************************************************/
352 /* ~BTDataset() */
353 /************************************************************************/
354
~BTDataset()355 BTDataset::~BTDataset()
356
357 {
358 FlushCache();
359 if( fpImage != NULL )
360 VSIFCloseL( fpImage );
361 CPLFree( pszProjection );
362 }
363
364 /************************************************************************/
365 /* FlushCache() */
366 /* */
367 /* We override this to include flush out the header block. */
368 /************************************************************************/
369
FlushCache()370 void BTDataset::FlushCache()
371
372 {
373 GDALDataset::FlushCache();
374
375 if( !bHeaderModified )
376 return;
377
378 bHeaderModified = FALSE;
379
380 VSIFSeekL( fpImage, 0, SEEK_SET );
381 VSIFWriteL( abyHeader, 256, 1, fpImage );
382 }
383
384 /************************************************************************/
385 /* GetGeoTransform() */
386 /************************************************************************/
387
GetGeoTransform(double * padfTransform)388 CPLErr BTDataset::GetGeoTransform( double * padfTransform )
389
390 {
391 memcpy( padfTransform, adfGeoTransform, sizeof(double)*6 );
392
393 if( bGeoTransformValid )
394 return CE_None;
395 else
396 return CE_Failure;
397 }
398
399 /************************************************************************/
400 /* SetGeoTransform() */
401 /************************************************************************/
402
SetGeoTransform(double * padfTransform)403 CPLErr BTDataset::SetGeoTransform( double *padfTransform )
404
405 {
406 CPLErr eErr = CE_None;
407
408 memcpy( adfGeoTransform, padfTransform, sizeof(double) * 6 );
409 if( adfGeoTransform[2] != 0.0 || adfGeoTransform[4] != 0.0 )
410 {
411 CPLError( CE_Failure, CPLE_AppDefined,
412 ".bt format does not support rotational coefficients in geotransform, ignoring." );
413 eErr = CE_Failure;
414 }
415
416 /* -------------------------------------------------------------------- */
417 /* Compute bounds, and update header info. */
418 /* -------------------------------------------------------------------- */
419 double dfLeft, dfRight, dfTop, dfBottom;
420
421 dfLeft = adfGeoTransform[0];
422 dfRight = dfLeft + adfGeoTransform[1] * nRasterXSize;
423 dfTop = adfGeoTransform[3];
424 dfBottom = dfTop + adfGeoTransform[5] * nRasterYSize;
425
426 memcpy( abyHeader + 28, &dfLeft, 8 );
427 memcpy( abyHeader + 36, &dfRight, 8 );
428 memcpy( abyHeader + 44, &dfBottom, 8 );
429 memcpy( abyHeader + 52, &dfTop, 8 );
430
431 CPL_LSBPTR64( abyHeader + 28 );
432 CPL_LSBPTR64( abyHeader + 36 );
433 CPL_LSBPTR64( abyHeader + 44 );
434 CPL_LSBPTR64( abyHeader + 52 );
435
436 bHeaderModified = TRUE;
437
438 return eErr;
439 }
440
441 /************************************************************************/
442 /* GetProjectionRef() */
443 /************************************************************************/
444
GetProjectionRef()445 const char *BTDataset::GetProjectionRef()
446
447 {
448 if( pszProjection == NULL )
449 return "";
450 else
451 return pszProjection;
452 }
453
454 /************************************************************************/
455 /* SetProjection() */
456 /************************************************************************/
457
SetProjection(const char * pszNewProjection)458 CPLErr BTDataset::SetProjection( const char *pszNewProjection )
459
460 {
461 CPLErr eErr = CE_None;
462
463 CPLFree( pszProjection );
464 pszProjection = CPLStrdup( pszNewProjection );
465
466 bHeaderModified = TRUE;
467
468 /* -------------------------------------------------------------------- */
469 /* Parse projection. */
470 /* -------------------------------------------------------------------- */
471 OGRSpatialReference oSRS( pszProjection );
472 GInt16 nShortTemp;
473
474 /* -------------------------------------------------------------------- */
475 /* Linear units. */
476 /* -------------------------------------------------------------------- */
477 if( oSRS.IsGeographic() )
478 nShortTemp = 0;
479 else
480 {
481 double dfLinear = oSRS.GetLinearUnits();
482
483 if( ABS(dfLinear - 0.3048) < 0.0000001 )
484 nShortTemp = 2;
485 else if( ABS(dfLinear - CPLAtof(SRS_UL_US_FOOT_CONV)) < 0.00000001 )
486 nShortTemp = 3;
487 else
488 nShortTemp = 1;
489 }
490
491 nShortTemp = CPL_LSBWORD16( 1 );
492 memcpy( abyHeader + 22, &nShortTemp, 2 );
493
494 /* -------------------------------------------------------------------- */
495 /* UTM Zone */
496 /* -------------------------------------------------------------------- */
497 int bNorth;
498
499 nShortTemp = (GInt16) oSRS.GetUTMZone( &bNorth );
500 if( bNorth )
501 nShortTemp = -nShortTemp;
502
503 nShortTemp = CPL_LSBWORD16( nShortTemp );
504 memcpy( abyHeader + 24, &nShortTemp, 2 );
505
506 /* -------------------------------------------------------------------- */
507 /* Datum */
508 /* -------------------------------------------------------------------- */
509 if( oSRS.GetAuthorityName( "GEOGCS|DATUM" ) != NULL
510 && EQUAL(oSRS.GetAuthorityName( "GEOGCS|DATUM" ),"EPSG") )
511 nShortTemp = (GInt16) (atoi(oSRS.GetAuthorityCode( "GEOGCS|DATUM" )) + 2000);
512 else
513 nShortTemp = -2;
514 nShortTemp = CPL_LSBWORD16( nShortTemp ); /* datum unknown */
515 memcpy( abyHeader + 26, &nShortTemp, 2 );
516
517 /* -------------------------------------------------------------------- */
518 /* Write out the projection to a .prj file. */
519 /* -------------------------------------------------------------------- */
520 const char *pszPrjFile = CPLResetExtension( GetDescription(), "prj" );
521 VSILFILE * fp;
522
523 fp = VSIFOpenL( pszPrjFile, "wt" );
524 if( fp != NULL )
525 {
526 VSIFPrintfL( fp, "%s\n", pszProjection );
527 VSIFCloseL( fp );
528 abyHeader[60] = 1;
529 }
530 else
531 {
532 CPLError( CE_Failure, CPLE_AppDefined,
533 "Unable to write out .prj file." );
534 eErr = CE_Failure;
535 }
536
537 return eErr;
538 }
539
540
541 /************************************************************************/
542 /* Open() */
543 /************************************************************************/
544
Open(GDALOpenInfo * poOpenInfo)545 GDALDataset *BTDataset::Open( GDALOpenInfo * poOpenInfo )
546
547 {
548 /* -------------------------------------------------------------------- */
549 /* Verify that this is some form of binterr file. */
550 /* -------------------------------------------------------------------- */
551 if( poOpenInfo->nHeaderBytes < 256)
552 return NULL;
553
554 if( strncmp( (const char *) poOpenInfo->pabyHeader, "binterr", 7 ) != 0 )
555 return NULL;
556
557 /* -------------------------------------------------------------------- */
558 /* Create the dataset. */
559 /* -------------------------------------------------------------------- */
560 BTDataset *poDS;
561
562 poDS = new BTDataset();
563
564 memcpy( poDS->abyHeader, poOpenInfo->pabyHeader, 256 );
565
566 /* -------------------------------------------------------------------- */
567 /* Get the version. */
568 /* -------------------------------------------------------------------- */
569 char szVersion[4];
570
571 strncpy( szVersion, (char *) (poDS->abyHeader + 7), 3 );
572 szVersion[3] = '\0';
573 poDS->nVersionCode = (int) (CPLAtof(szVersion) * 10);
574
575 /* -------------------------------------------------------------------- */
576 /* Extract core header information, being careful about the */
577 /* version. */
578 /* -------------------------------------------------------------------- */
579 GInt32 nIntTemp;
580 GInt16 nDataSize;
581 GDALDataType eType;
582
583 memcpy( &nIntTemp, poDS->abyHeader + 10, 4 );
584 poDS->nRasterXSize = CPL_LSBWORD32( nIntTemp );
585
586 memcpy( &nIntTemp, poDS->abyHeader + 14, 4 );
587 poDS->nRasterYSize = CPL_LSBWORD32( nIntTemp );
588
589 if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
590 {
591 delete poDS;
592 return NULL;
593 }
594
595 memcpy( &nDataSize, poDS->abyHeader+18, 2 );
596 nDataSize = CPL_LSBWORD16( nDataSize );
597
598 if( poDS->abyHeader[20] != 0 && nDataSize == 4 )
599 eType = GDT_Float32;
600 else if( poDS->abyHeader[20] == 0 && nDataSize == 4 )
601 eType = GDT_Int32;
602 else if( poDS->abyHeader[20] == 0 && nDataSize == 2 )
603 eType = GDT_Int16;
604 else
605 {
606 CPLError( CE_Failure, CPLE_AppDefined,
607 ".bt file data type unknown, got datasize=%d.",
608 nDataSize );
609 delete poDS;
610 return NULL;
611 }
612
613 /*
614 rcg, apr 7/06: read offset 62 for vert. units.
615 If zero, assume 1.0 as per spec.
616
617 */
618 memcpy( &poDS->m_fVscale, poDS->abyHeader + 62, 4 );
619 CPL_LSBPTR32(&poDS->m_fVscale);
620 if(poDS->m_fVscale == 0.0f)
621 poDS->m_fVscale = 1.0f;
622
623 /* -------------------------------------------------------------------- */
624 /* Try to read a .prj file if it is indicated. */
625 /* -------------------------------------------------------------------- */
626 OGRSpatialReference oSRS;
627
628 if( poDS->nVersionCode >= 12 && poDS->abyHeader[60] != 0 )
629 {
630 const char *pszPrjFile = CPLResetExtension( poOpenInfo->pszFilename,
631 "prj" );
632 VSILFILE *fp;
633
634 fp = VSIFOpenL( pszPrjFile, "rt" );
635 if( fp != NULL )
636 {
637 char *pszBuffer, *pszBufPtr;
638 int nBufMax = 10000;
639 int nBytes;
640
641 pszBuffer = (char *) CPLMalloc(nBufMax);
642 nBytes = VSIFReadL( pszBuffer, 1, nBufMax-1, fp );
643 VSIFCloseL( fp );
644
645 pszBuffer[nBytes] = '\0';
646
647 pszBufPtr = pszBuffer;
648 if( oSRS.importFromWkt( &pszBufPtr ) != OGRERR_NONE )
649 {
650 CPLError( CE_Warning, CPLE_AppDefined,
651 "Unable to parse .prj file, coordinate system missing." );
652 }
653 CPLFree( pszBuffer );
654 }
655 }
656
657 /* -------------------------------------------------------------------- */
658 /* If we didn't find a .prj file, try to use internal info. */
659 /* -------------------------------------------------------------------- */
660 if( oSRS.GetRoot() == NULL )
661 {
662 GInt16 nUTMZone, nDatum, nHUnits;
663
664 memcpy( &nUTMZone, poDS->abyHeader + 24, 2 );
665 nUTMZone = CPL_LSBWORD16( nUTMZone );
666
667 memcpy( &nDatum, poDS->abyHeader + 26, 2 );
668 nDatum = CPL_LSBWORD16( nDatum );
669
670 memcpy( &nHUnits, poDS->abyHeader + 22, 2 );
671 nHUnits = CPL_LSBWORD16( nHUnits );
672
673 if( nUTMZone != 0 )
674 oSRS.SetUTM( ABS(nUTMZone), nUTMZone > 0 );
675 else if( nHUnits != 0 )
676 oSRS.SetLocalCS( "Unknown" );
677
678 if( nHUnits == 1 )
679 oSRS.SetLinearUnits( SRS_UL_METER, 1.0 );
680 else if( nHUnits == 2 )
681 oSRS.SetLinearUnits( SRS_UL_FOOT, CPLAtof(SRS_UL_FOOT_CONV) );
682 else if( nHUnits == 3 )
683 oSRS.SetLinearUnits( SRS_UL_US_FOOT, CPLAtof(SRS_UL_US_FOOT_CONV) );
684
685 // Translate some of the more obvious old USGS datum codes
686 if( nDatum == 0 )
687 nDatum = 6201;
688 else if( nDatum == 1 )
689 nDatum = 6209;
690 else if( nDatum == 2 )
691 nDatum = 6210;
692 else if( nDatum == 3 )
693 nDatum = 6202;
694 else if( nDatum == 4 )
695 nDatum = 6203;
696 else if( nDatum == 6 )
697 nDatum = 6222;
698 else if( nDatum == 7 )
699 nDatum = 6230;
700 else if( nDatum == 13 )
701 nDatum = 6267;
702 else if( nDatum == 14 )
703 nDatum = 6269;
704 else if( nDatum == 17 )
705 nDatum = 6277;
706 else if( nDatum == 19 )
707 nDatum = 6284;
708 else if( nDatum == 21 )
709 nDatum = 6301;
710 else if( nDatum == 22 )
711 nDatum = 6322;
712 else if( nDatum == 23 )
713 nDatum = 6326;
714
715 if( !oSRS.IsLocal() )
716 {
717 if( nDatum >= 6000 )
718 {
719 char szName[32];
720 sprintf( szName, "EPSG:%d", nDatum-2000 );
721 oSRS.SetWellKnownGeogCS( szName );
722 }
723 else
724 oSRS.SetWellKnownGeogCS( "WGS84" );
725 }
726 }
727
728 /* -------------------------------------------------------------------- */
729 /* Convert coordinate system back to WKT. */
730 /* -------------------------------------------------------------------- */
731 if( oSRS.GetRoot() != NULL )
732 oSRS.exportToWkt( &poDS->pszProjection );
733
734 /* -------------------------------------------------------------------- */
735 /* Get georeferencing bounds. */
736 /* -------------------------------------------------------------------- */
737 if( poDS->nVersionCode >= 11 )
738 {
739 double dfLeft, dfRight, dfTop, dfBottom;
740
741 memcpy( &dfLeft, poDS->abyHeader + 28, 8 );
742 CPL_LSBPTR64( &dfLeft );
743
744 memcpy( &dfRight, poDS->abyHeader + 36, 8 );
745 CPL_LSBPTR64( &dfRight );
746
747 memcpy( &dfBottom, poDS->abyHeader + 44, 8 );
748 CPL_LSBPTR64( &dfBottom );
749
750 memcpy( &dfTop, poDS->abyHeader + 52, 8 );
751 CPL_LSBPTR64( &dfTop );
752
753 poDS->adfGeoTransform[0] = dfLeft;
754 poDS->adfGeoTransform[1] = (dfRight - dfLeft) / poDS->nRasterXSize;
755 poDS->adfGeoTransform[2] = 0.0;
756 poDS->adfGeoTransform[3] = dfTop;
757 poDS->adfGeoTransform[4] = 0.0;
758 poDS->adfGeoTransform[5] = (dfBottom - dfTop) / poDS->nRasterYSize;
759
760 poDS->bGeoTransformValid = TRUE;
761 }
762
763 /* -------------------------------------------------------------------- */
764 /* Re-open the file with the desired access. */
765 /* -------------------------------------------------------------------- */
766
767 if( poOpenInfo->eAccess == GA_Update )
768 poDS->fpImage = VSIFOpenL( poOpenInfo->pszFilename, "rb+" );
769 else
770 poDS->fpImage = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
771
772 if( poDS->fpImage == NULL )
773 {
774 CPLError( CE_Failure, CPLE_OpenFailed,
775 "Failed to re-open %s within BT driver.\n",
776 poOpenInfo->pszFilename );
777 return NULL;
778 }
779 poDS->eAccess = poOpenInfo->eAccess;
780
781 /* -------------------------------------------------------------------- */
782 /* Create band information objects */
783 /* -------------------------------------------------------------------- */
784 poDS->SetBand( 1, new BTRasterBand( poDS, poDS->fpImage, eType ) );
785
786 #ifdef notdef
787 poDS->bGeoTransformValid =
788 GDALReadWorldFile( poOpenInfo->pszFilename, ".wld",
789 poDS->adfGeoTransform );
790 #endif
791
792 /* -------------------------------------------------------------------- */
793 /* Initialize any PAM information. */
794 /* -------------------------------------------------------------------- */
795 poDS->SetDescription( poOpenInfo->pszFilename );
796 poDS->TryLoadXML();
797
798 /* -------------------------------------------------------------------- */
799 /* Check for overviews. */
800 /* -------------------------------------------------------------------- */
801 poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
802
803 return( poDS );
804 }
805
806 /************************************************************************/
807 /* Create() */
808 /************************************************************************/
809
Create(const char * pszFilename,int nXSize,int nYSize,int nBands,GDALDataType eType,CPL_UNUSED char ** papszOptions)810 GDALDataset *BTDataset::Create( const char * pszFilename,
811 int nXSize,
812 int nYSize,
813 int nBands,
814 GDALDataType eType,
815 CPL_UNUSED char ** papszOptions )
816 {
817
818 /* -------------------------------------------------------------------- */
819 /* Verify input options. */
820 /* -------------------------------------------------------------------- */
821 if( eType != GDT_Int16 && eType != GDT_Int32 && eType != GDT_Float32 )
822 {
823 CPLError( CE_Failure, CPLE_AppDefined,
824 "Attempt to create .bt dataset with an illegal\n"
825 "data type (%s), only Int16, Int32 and Float32 supported.\n",
826 GDALGetDataTypeName(eType) );
827
828 return NULL;
829 }
830
831 if( nBands != 1 )
832 {
833 CPLError( CE_Failure, CPLE_AppDefined,
834 "Attempt to create .bt dataset with %d bands, only 1 supported",
835 nBands );
836
837 return NULL;
838 }
839
840 /* -------------------------------------------------------------------- */
841 /* Try to create the file. */
842 /* -------------------------------------------------------------------- */
843 VSILFILE *fp;
844
845 fp = VSIFOpenL( pszFilename, "wb" );
846
847 if( fp == NULL )
848 {
849 CPLError( CE_Failure, CPLE_OpenFailed,
850 "Attempt to create file `%s' failed.\n",
851 pszFilename );
852 return NULL;
853 }
854
855 /* -------------------------------------------------------------------- */
856 /* Setup base header. */
857 /* -------------------------------------------------------------------- */
858 unsigned char abyHeader[256];
859 GInt32 nTemp;
860 GInt16 nShortTemp;
861
862 memset( abyHeader, 0, 256 );
863 memcpy( abyHeader, "binterr1.3", 10 );
864
865 nTemp = CPL_LSBWORD32( nXSize );
866 memcpy( abyHeader+10, &nTemp, 4 );
867
868 nTemp = CPL_LSBWORD32( nYSize );
869 memcpy( abyHeader+14, &nTemp, 4 );
870
871 nShortTemp = (GInt16) (CPL_LSBWORD16( GDALGetDataTypeSize( eType ) / 8 ));
872 memcpy( abyHeader + 18, &nShortTemp, 2 );
873
874 if( eType == GDT_Float32 )
875 abyHeader[20] = 1;
876 else
877 abyHeader[20] = 0;
878
879 nShortTemp = CPL_LSBWORD16( 1 ); /* meters */
880 memcpy( abyHeader + 22, &nShortTemp, 2 );
881
882 nShortTemp = CPL_LSBWORD16( 0 ); /* not utm */
883 memcpy( abyHeader + 24, &nShortTemp, 2 );
884
885 nShortTemp = CPL_LSBWORD16( -2 ); /* datum unknown */
886 memcpy( abyHeader + 26, &nShortTemp, 2 );
887
888 /* -------------------------------------------------------------------- */
889 /* Set dummy extents. */
890 /* -------------------------------------------------------------------- */
891 double dfLeft=0, dfRight=nXSize, dfTop=nYSize, dfBottom=0;
892
893 memcpy( abyHeader + 28, &dfLeft, 8 );
894 memcpy( abyHeader + 36, &dfRight, 8 );
895 memcpy( abyHeader + 44, &dfBottom, 8 );
896 memcpy( abyHeader + 52, &dfTop, 8 );
897
898 CPL_LSBPTR64( abyHeader + 28 );
899 CPL_LSBPTR64( abyHeader + 36 );
900 CPL_LSBPTR64( abyHeader + 44 );
901 CPL_LSBPTR64( abyHeader + 52 );
902
903 /* -------------------------------------------------------------------- */
904 /* Set dummy scale. */
905 /* -------------------------------------------------------------------- */
906 float fScale = 1.0;
907
908 memcpy( abyHeader + 62, &fScale, 4 );
909 CPL_LSBPTR32( abyHeader + 62 );
910
911 /* -------------------------------------------------------------------- */
912 /* Write to disk. */
913 /* -------------------------------------------------------------------- */
914 VSIFWriteL( (void *) abyHeader, 256, 1, fp );
915 if( VSIFSeekL( fp, (GDALGetDataTypeSize(eType)/8) * nXSize * (vsi_l_offset)nYSize - 1,
916 SEEK_CUR ) != 0
917 || VSIFWriteL( abyHeader+255, 1, 1, fp ) != 1 )
918 {
919 CPLError( CE_Failure, CPLE_FileIO,
920 "Failed to extent file to its full size, out of disk space?"
921 );
922
923 VSIFCloseL( fp );
924 VSIUnlink( pszFilename );
925 return NULL;
926 }
927
928 VSIFCloseL( fp );
929
930 return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
931 }
932
933 /************************************************************************/
934 /* GDALRegister_BT() */
935 /************************************************************************/
936
GDALRegister_BT()937 void GDALRegister_BT()
938
939 {
940 GDALDriver *poDriver;
941
942 if( GDALGetDriverByName( "BT" ) == NULL )
943 {
944 poDriver = new GDALDriver();
945
946 poDriver->SetDescription( "BT" );
947 poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
948 poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
949 "VTP .bt (Binary Terrain) 1.3 Format" );
950 poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
951 "frmt_various.html#BT" );
952 poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "bt" );
953 poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
954 "Int16 Int32 Float32" );
955 poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
956
957 poDriver->pfnOpen = BTDataset::Open;
958 poDriver->pfnCreate = BTDataset::Create;
959
960 GetGDALDriverManager()->RegisterDriver( poDriver );
961 }
962 }
963