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