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