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