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