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 : GDALGeoTransform m_gt{};
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(GDALGeoTransform &) const override;
58 : CPLErr SetGeoTransform(const GDALGeoTransform &) 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 20 : double BTRasterBand::GetNoDataValue(int *pbSuccess /*= NULL */)
242 : {
243 : // First check in PAM
244 20 : int bSuccess = FALSE;
245 20 : double dfRet = GDALPamRasterBand::GetNoDataValue(&bSuccess);
246 20 : if (bSuccess)
247 : {
248 0 : if (pbSuccess != nullptr)
249 0 : *pbSuccess = TRUE;
250 0 : return dfRet;
251 : }
252 :
253 : // Otherwise defaults to -32768
254 20 : if (pbSuccess != nullptr)
255 16 : *pbSuccess = TRUE;
256 :
257 20 : 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 : memset(abyHeader, 0, sizeof(abyHeader));
354 27 : }
355 :
356 : /************************************************************************/
357 : /* ~BTDataset() */
358 : /************************************************************************/
359 :
360 54 : BTDataset::~BTDataset()
361 :
362 : {
363 27 : BTDataset::FlushCache(true);
364 27 : if (fpImage != nullptr)
365 : {
366 27 : if (VSIFCloseL(fpImage) != 0)
367 : {
368 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
369 : }
370 : }
371 54 : }
372 :
373 : /************************************************************************/
374 : /* FlushCache() */
375 : /* */
376 : /* We override this to include flush out the header block. */
377 : /************************************************************************/
378 :
379 28 : CPLErr BTDataset::FlushCache(bool bAtClosing)
380 :
381 : {
382 28 : CPLErr eErr = GDALDataset::FlushCache(bAtClosing);
383 :
384 28 : if (!bHeaderModified)
385 17 : return eErr;
386 :
387 11 : bHeaderModified = FALSE;
388 :
389 22 : if (VSIFSeekL(fpImage, 0, SEEK_SET) != 0 ||
390 11 : VSIFWriteL(abyHeader, 256, 1, fpImage) != 1)
391 : {
392 0 : eErr = CE_Failure;
393 : }
394 11 : return eErr;
395 : }
396 :
397 : /************************************************************************/
398 : /* GetGeoTransform() */
399 : /************************************************************************/
400 :
401 7 : CPLErr BTDataset::GetGeoTransform(GDALGeoTransform >) const
402 :
403 : {
404 7 : gt = m_gt;
405 :
406 7 : if (bGeoTransformValid)
407 7 : return CE_None;
408 : else
409 0 : return CE_Failure;
410 : }
411 :
412 : /************************************************************************/
413 : /* SetGeoTransform() */
414 : /************************************************************************/
415 :
416 11 : CPLErr BTDataset::SetGeoTransform(const GDALGeoTransform >)
417 :
418 : {
419 11 : CPLErr eErr = CE_None;
420 :
421 11 : m_gt = gt;
422 11 : if (gt[2] != 0.0 || gt[4] != 0.0)
423 : {
424 0 : CPLError(CE_Failure, CPLE_AppDefined,
425 : ".bt format does not support rotational coefficients "
426 : "in geotransform, ignoring.");
427 0 : eErr = CE_Failure;
428 : }
429 :
430 : /* -------------------------------------------------------------------- */
431 : /* Compute bounds, and update header info. */
432 : /* -------------------------------------------------------------------- */
433 11 : const double dfLeft = m_gt[0];
434 11 : const double dfRight = dfLeft + m_gt[1] * nRasterXSize;
435 11 : const double dfTop = m_gt[3];
436 11 : const double dfBottom = dfTop + m_gt[5] * nRasterYSize;
437 :
438 11 : memcpy(abyHeader + 28, &dfLeft, 8);
439 11 : memcpy(abyHeader + 36, &dfRight, 8);
440 11 : memcpy(abyHeader + 44, &dfBottom, 8);
441 11 : memcpy(abyHeader + 52, &dfTop, 8);
442 :
443 11 : CPL_LSBPTR64(abyHeader + 28);
444 11 : CPL_LSBPTR64(abyHeader + 36);
445 11 : CPL_LSBPTR64(abyHeader + 44);
446 11 : CPL_LSBPTR64(abyHeader + 52);
447 :
448 11 : bHeaderModified = TRUE;
449 :
450 11 : return eErr;
451 : }
452 :
453 : /************************************************************************/
454 : /* SetSpatialRef() */
455 : /************************************************************************/
456 :
457 10 : CPLErr BTDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
458 :
459 : {
460 10 : CPLErr eErr = CE_None;
461 :
462 10 : if (poSRS)
463 10 : m_oSRS = *poSRS;
464 : else
465 0 : m_oSRS.Clear();
466 :
467 10 : bHeaderModified = TRUE;
468 :
469 : /* -------------------------------------------------------------------- */
470 : /* Parse projection. */
471 : /* -------------------------------------------------------------------- */
472 :
473 : /* -------------------------------------------------------------------- */
474 : /* Linear units. */
475 : /* -------------------------------------------------------------------- */
476 : #if 0
477 : if( m_oSRS.IsGeographic() )
478 : {
479 : // nShortTemp = 0;
480 : }
481 : else
482 : {
483 : const double dfLinear = m_oSRS.GetLinearUnits();
484 :
485 : if( std::abs(dfLinear - 0.3048) < 0.0000001 )
486 : nShortTemp = 2;
487 : else if( std::abs(dfLinear - CPLAtof(SRS_UL_US_FOOT_CONV))
488 : < 0.00000001 )
489 : nShortTemp = 3;
490 : else
491 : nShortTemp = 1;
492 : }
493 : #endif
494 10 : GInt16 nShortTemp = CPL_LSBWORD16(1);
495 10 : memcpy(abyHeader + 22, &nShortTemp, 2);
496 :
497 : /* -------------------------------------------------------------------- */
498 : /* UTM Zone */
499 : /* -------------------------------------------------------------------- */
500 10 : int bNorth = FALSE;
501 :
502 10 : nShortTemp = static_cast<GInt16>(m_oSRS.GetUTMZone(&bNorth));
503 10 : if (bNorth)
504 3 : nShortTemp = -nShortTemp;
505 :
506 10 : CPL_LSBPTR16(&nShortTemp);
507 10 : memcpy(abyHeader + 24, &nShortTemp, 2);
508 :
509 : /* -------------------------------------------------------------------- */
510 : /* Datum */
511 : /* -------------------------------------------------------------------- */
512 14 : if (m_oSRS.GetAuthorityName("GEOGCS|DATUM") != nullptr &&
513 4 : EQUAL(m_oSRS.GetAuthorityName("GEOGCS|DATUM"), "EPSG"))
514 4 : nShortTemp = static_cast<GInt16>(
515 4 : atoi(m_oSRS.GetAuthorityCode("GEOGCS|DATUM")) + 2000);
516 : else
517 6 : nShortTemp = -2;
518 10 : CPL_LSBPTR16(&nShortTemp); /* datum unknown */
519 10 : memcpy(abyHeader + 26, &nShortTemp, 2);
520 :
521 : /* -------------------------------------------------------------------- */
522 : /* Write out the projection to a .prj file. */
523 : /* -------------------------------------------------------------------- */
524 10 : char *pszProjection = nullptr;
525 10 : const char *const apszOptions[] = {"FORMAT=WKT1", nullptr};
526 10 : m_oSRS.exportToWkt(&pszProjection, apszOptions);
527 10 : if (pszProjection)
528 : {
529 : const std::string osPrjFile =
530 20 : CPLResetExtensionSafe(GetDescription(), "prj");
531 10 : VSILFILE *fp = VSIFOpenL(osPrjFile.c_str(), "wt");
532 10 : if (fp != nullptr)
533 : {
534 10 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fp, "%s\n", pszProjection));
535 10 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
536 10 : abyHeader[60] = 1;
537 : }
538 : else
539 : {
540 0 : CPLError(CE_Failure, CPLE_AppDefined,
541 : "Unable to write out .prj file.");
542 0 : eErr = CE_Failure;
543 : }
544 10 : CPLFree(pszProjection);
545 : }
546 :
547 10 : return eErr;
548 : }
549 :
550 : /************************************************************************/
551 : /* Open() */
552 : /************************************************************************/
553 :
554 33915 : GDALDataset *BTDataset::Open(GDALOpenInfo *poOpenInfo)
555 :
556 : {
557 : /* -------------------------------------------------------------------- */
558 : /* Verify that this is some form of binterr file. */
559 : /* -------------------------------------------------------------------- */
560 33915 : if (poOpenInfo->nHeaderBytes < 256 || poOpenInfo->fpL == nullptr)
561 31013 : return nullptr;
562 :
563 2902 : if (!STARTS_WITH(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
564 : "binterr"))
565 2875 : return nullptr;
566 :
567 : /* -------------------------------------------------------------------- */
568 : /* Create the dataset. */
569 : /* -------------------------------------------------------------------- */
570 27 : BTDataset *poDS = new BTDataset();
571 :
572 27 : memcpy(poDS->abyHeader, poOpenInfo->pabyHeader, 256);
573 :
574 : /* -------------------------------------------------------------------- */
575 : /* Get the version. */
576 : /* -------------------------------------------------------------------- */
577 27 : char szVersion[4] = {};
578 :
579 27 : strncpy(szVersion, reinterpret_cast<char *>(poDS->abyHeader + 7), 3);
580 27 : szVersion[3] = '\0';
581 27 : poDS->nVersionCode = static_cast<int>(CPLAtof(szVersion) * 10);
582 :
583 : /* -------------------------------------------------------------------- */
584 : /* Extract core header information, being careful about the */
585 : /* version. */
586 : /* -------------------------------------------------------------------- */
587 :
588 27 : GInt32 nIntTemp = 0;
589 27 : memcpy(&nIntTemp, poDS->abyHeader + 10, 4);
590 27 : poDS->nRasterXSize = CPL_LSBWORD32(nIntTemp);
591 :
592 27 : memcpy(&nIntTemp, poDS->abyHeader + 14, 4);
593 27 : poDS->nRasterYSize = CPL_LSBWORD32(nIntTemp);
594 :
595 27 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
596 : {
597 0 : delete poDS;
598 0 : return nullptr;
599 : }
600 :
601 27 : GInt16 nDataSize = 0;
602 27 : memcpy(&nDataSize, poDS->abyHeader + 18, 2);
603 27 : CPL_LSBPTR16(&nDataSize);
604 :
605 27 : GDALDataType eType = GDT_Unknown;
606 27 : if (poDS->abyHeader[20] != 0 && nDataSize == 4)
607 15 : eType = GDT_Float32;
608 12 : else if (poDS->abyHeader[20] == 0 && nDataSize == 4)
609 6 : eType = GDT_Int32;
610 6 : else if (poDS->abyHeader[20] == 0 && nDataSize == 2)
611 6 : eType = GDT_Int16;
612 : else
613 : {
614 0 : CPLError(CE_Failure, CPLE_AppDefined,
615 : ".bt file data type unknown, got datasize=%d.", nDataSize);
616 0 : delete poDS;
617 0 : return nullptr;
618 : }
619 :
620 : /*
621 : rcg, apr 7/06: read offset 62 for vert. units.
622 : If zero, assume 1.0 as per spec.
623 :
624 : */
625 27 : memcpy(&poDS->m_fVscale, poDS->abyHeader + 62, 4);
626 27 : CPL_LSBPTR32(&poDS->m_fVscale);
627 27 : if (poDS->m_fVscale == 0.0f)
628 0 : poDS->m_fVscale = 1.0f;
629 :
630 : /* -------------------------------------------------------------------- */
631 : /* Try to read a .prj file if it is indicated. */
632 : /* -------------------------------------------------------------------- */
633 27 : OGRSpatialReference oSRS;
634 27 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
635 :
636 27 : if (poDS->nVersionCode >= 12 && poDS->abyHeader[60] != 0)
637 : {
638 : const std::string osPrjFile =
639 22 : CPLResetExtensionSafe(poOpenInfo->pszFilename, "prj");
640 11 : VSILFILE *fp = VSIFOpenL(osPrjFile.c_str(), "rt");
641 11 : if (fp != nullptr)
642 : {
643 11 : const int nBufMax = 10000;
644 :
645 11 : char *pszBuffer = static_cast<char *>(CPLMalloc(nBufMax));
646 : const int nBytes =
647 11 : static_cast<int>(VSIFReadL(pszBuffer, 1, nBufMax - 1, fp));
648 11 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
649 :
650 11 : pszBuffer[nBytes] = '\0';
651 :
652 11 : if (oSRS.importFromWkt(pszBuffer) != OGRERR_NONE)
653 : {
654 0 : CPLError(CE_Warning, CPLE_AppDefined,
655 : "Unable to parse .prj file, "
656 : "coordinate system missing.");
657 : }
658 11 : CPLFree(pszBuffer);
659 : }
660 : }
661 :
662 : /* -------------------------------------------------------------------- */
663 : /* If we didn't find a .prj file, try to use internal info. */
664 : /* -------------------------------------------------------------------- */
665 27 : if (oSRS.GetRoot() == nullptr)
666 : {
667 16 : GInt16 nUTMZone = 0;
668 16 : memcpy(&nUTMZone, poDS->abyHeader + 24, 2);
669 16 : CPL_LSBPTR16(&nUTMZone);
670 :
671 16 : GInt16 nDatum = 0;
672 16 : memcpy(&nDatum, poDS->abyHeader + 26, 2);
673 16 : CPL_LSBPTR16(&nDatum);
674 :
675 16 : GInt16 nHUnits = 0;
676 16 : memcpy(&nHUnits, poDS->abyHeader + 22, 2);
677 16 : CPL_LSBPTR16(&nHUnits);
678 :
679 16 : if (nUTMZone != 0)
680 0 : oSRS.SetUTM(std::abs(static_cast<int>(nUTMZone)), nUTMZone > 0);
681 16 : else if (nHUnits != 0)
682 16 : oSRS.SetLocalCS("Unknown");
683 :
684 16 : if (nHUnits == 1)
685 16 : oSRS.SetLinearUnits(SRS_UL_METER, 1.0);
686 0 : else if (nHUnits == 2)
687 0 : oSRS.SetLinearUnits(SRS_UL_FOOT, CPLAtof(SRS_UL_FOOT_CONV));
688 0 : else if (nHUnits == 3)
689 0 : oSRS.SetLinearUnits(SRS_UL_US_FOOT, CPLAtof(SRS_UL_US_FOOT_CONV));
690 :
691 : // Translate some of the more obvious old USGS datum codes
692 16 : if (nDatum == 0)
693 0 : nDatum = 6201;
694 16 : else if (nDatum == 1)
695 0 : nDatum = 6209;
696 16 : else if (nDatum == 2)
697 0 : nDatum = 6210;
698 16 : else if (nDatum == 3)
699 0 : nDatum = 6202;
700 16 : else if (nDatum == 4)
701 0 : nDatum = 6203;
702 16 : else if (nDatum == 6)
703 0 : nDatum = 6222;
704 16 : else if (nDatum == 7)
705 0 : nDatum = 6230;
706 16 : else if (nDatum == 13)
707 0 : nDatum = 6267;
708 16 : else if (nDatum == 14)
709 0 : nDatum = 6269;
710 16 : else if (nDatum == 17)
711 0 : nDatum = 6277;
712 16 : else if (nDatum == 19)
713 0 : nDatum = 6284;
714 16 : else if (nDatum == 21)
715 0 : nDatum = 6301;
716 16 : else if (nDatum == 22)
717 0 : nDatum = 6322;
718 16 : else if (nDatum == 23)
719 0 : nDatum = 6326;
720 :
721 16 : if (!oSRS.IsLocal())
722 : {
723 0 : if (nDatum >= 6000)
724 : {
725 : char szName[32];
726 0 : snprintf(szName, sizeof(szName), "EPSG:%d", nDatum - 2000);
727 0 : oSRS.SetWellKnownGeogCS(szName);
728 : }
729 : else
730 0 : oSRS.SetWellKnownGeogCS("WGS84");
731 : }
732 : }
733 :
734 : /* -------------------------------------------------------------------- */
735 : /* Convert coordinate system back to WKT. */
736 : /* -------------------------------------------------------------------- */
737 27 : if (oSRS.GetRoot() != nullptr)
738 27 : poDS->m_oSRS = std::move(oSRS);
739 :
740 : /* -------------------------------------------------------------------- */
741 : /* Get georeferencing bounds. */
742 : /* -------------------------------------------------------------------- */
743 27 : if (poDS->nVersionCode >= 11)
744 : {
745 27 : double dfLeft = 0.0;
746 27 : memcpy(&dfLeft, poDS->abyHeader + 28, 8);
747 27 : CPL_LSBPTR64(&dfLeft);
748 :
749 27 : double dfRight = 0.0;
750 27 : memcpy(&dfRight, poDS->abyHeader + 36, 8);
751 27 : CPL_LSBPTR64(&dfRight);
752 :
753 27 : double dfBottom = 0.0;
754 27 : memcpy(&dfBottom, poDS->abyHeader + 44, 8);
755 27 : CPL_LSBPTR64(&dfBottom);
756 :
757 27 : double dfTop = 0.0;
758 27 : memcpy(&dfTop, poDS->abyHeader + 52, 8);
759 27 : CPL_LSBPTR64(&dfTop);
760 :
761 27 : poDS->m_gt[0] = dfLeft;
762 27 : poDS->m_gt[1] = (dfRight - dfLeft) / poDS->nRasterXSize;
763 27 : poDS->m_gt[2] = 0.0;
764 27 : poDS->m_gt[3] = dfTop;
765 27 : poDS->m_gt[4] = 0.0;
766 27 : poDS->m_gt[5] = (dfBottom - dfTop) / poDS->nRasterYSize;
767 :
768 27 : poDS->bGeoTransformValid = TRUE;
769 : }
770 :
771 27 : poDS->eAccess = poOpenInfo->eAccess;
772 27 : poDS->fpImage = poOpenInfo->fpL;
773 27 : poOpenInfo->fpL = nullptr;
774 :
775 : /* -------------------------------------------------------------------- */
776 : /* Create band information objects */
777 : /* -------------------------------------------------------------------- */
778 27 : poDS->SetBand(1, new BTRasterBand(poDS, poDS->fpImage, eType));
779 :
780 : #ifdef notdef
781 : poDS->bGeoTransformValid =
782 : GDALReadWorldFile(poOpenInfo->pszFilename, ".wld", m_gt.data());
783 : #endif
784 :
785 : /* -------------------------------------------------------------------- */
786 : /* Initialize any PAM information. */
787 : /* -------------------------------------------------------------------- */
788 27 : poDS->SetDescription(poOpenInfo->pszFilename);
789 27 : poDS->TryLoadXML();
790 :
791 : /* -------------------------------------------------------------------- */
792 : /* Check for overviews. */
793 : /* -------------------------------------------------------------------- */
794 27 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
795 :
796 27 : return poDS;
797 : }
798 :
799 : /************************************************************************/
800 : /* Create() */
801 : /************************************************************************/
802 :
803 66 : GDALDataset *BTDataset::Create(const char *pszFilename, int nXSize, int nYSize,
804 : int nBandsIn, GDALDataType eType,
805 : CPL_UNUSED char **papszOptions)
806 : {
807 :
808 : /* -------------------------------------------------------------------- */
809 : /* Verify input options. */
810 : /* -------------------------------------------------------------------- */
811 66 : if (eType != GDT_Int16 && eType != GDT_Int32 && eType != GDT_Float32)
812 : {
813 38 : CPLError(CE_Failure, CPLE_AppDefined,
814 : "Attempt to create .bt dataset with an illegal "
815 : "data type (%s), only Int16, Int32 and Float32 supported.",
816 : GDALGetDataTypeName(eType));
817 :
818 38 : return nullptr;
819 : }
820 :
821 28 : if (nBandsIn != 1)
822 : {
823 3 : CPLError(
824 : CE_Failure, CPLE_AppDefined,
825 : "Attempt to create .bt dataset with %d bands, only 1 supported",
826 : nBandsIn);
827 :
828 3 : return nullptr;
829 : }
830 :
831 : /* -------------------------------------------------------------------- */
832 : /* Try to create the file. */
833 : /* -------------------------------------------------------------------- */
834 25 : VSILFILE *fp = VSIFOpenL(pszFilename, "wb");
835 :
836 25 : if (fp == nullptr)
837 : {
838 3 : CPLError(CE_Failure, CPLE_OpenFailed,
839 : "Attempt to create file `%s' failed.", pszFilename);
840 3 : return nullptr;
841 : }
842 :
843 : /* -------------------------------------------------------------------- */
844 : /* Setup base header. */
845 : /* -------------------------------------------------------------------- */
846 22 : unsigned char abyHeader[256] = {};
847 :
848 22 : memcpy(abyHeader, "binterr1.3", 10);
849 :
850 22 : GInt32 nTemp = CPL_LSBWORD32(nXSize);
851 22 : memcpy(abyHeader + 10, &nTemp, 4);
852 :
853 22 : nTemp = CPL_LSBWORD32(nYSize);
854 22 : memcpy(abyHeader + 14, &nTemp, 4);
855 :
856 : GInt16 nShortTemp = static_cast<GInt16>(
857 22 : CPL_LSBWORD16(static_cast<GInt16>(GDALGetDataTypeSizeBytes(eType))));
858 22 : memcpy(abyHeader + 18, &nShortTemp, 2);
859 :
860 22 : if (eType == GDT_Float32)
861 6 : abyHeader[20] = 1;
862 : else
863 16 : abyHeader[20] = 0;
864 :
865 22 : nShortTemp = CPL_LSBWORD16(1); /* meters */
866 22 : memcpy(abyHeader + 22, &nShortTemp, 2);
867 :
868 22 : nShortTemp = CPL_LSBWORD16(0); /* not utm */
869 22 : memcpy(abyHeader + 24, &nShortTemp, 2);
870 :
871 22 : nShortTemp = CPL_LSBWORD16(-2); /* datum unknown */
872 22 : memcpy(abyHeader + 26, &nShortTemp, 2);
873 :
874 : /* -------------------------------------------------------------------- */
875 : /* Set dummy extents. */
876 : /* -------------------------------------------------------------------- */
877 22 : double dfLeft = 0.0;
878 22 : double dfRight = nXSize;
879 22 : double dfTop = nYSize;
880 22 : double dfBottom = 0.0;
881 :
882 22 : memcpy(abyHeader + 28, &dfLeft, 8);
883 22 : memcpy(abyHeader + 36, &dfRight, 8);
884 22 : memcpy(abyHeader + 44, &dfBottom, 8);
885 22 : memcpy(abyHeader + 52, &dfTop, 8);
886 :
887 22 : CPL_LSBPTR64(abyHeader + 28);
888 22 : CPL_LSBPTR64(abyHeader + 36);
889 22 : CPL_LSBPTR64(abyHeader + 44);
890 22 : CPL_LSBPTR64(abyHeader + 52);
891 :
892 : /* -------------------------------------------------------------------- */
893 : /* Set dummy scale. */
894 : /* -------------------------------------------------------------------- */
895 22 : float fScale = 1.0;
896 :
897 22 : memcpy(abyHeader + 62, &fScale, 4);
898 22 : CPL_LSBPTR32(abyHeader + 62);
899 :
900 : /* -------------------------------------------------------------------- */
901 : /* Write to disk. */
902 : /* -------------------------------------------------------------------- */
903 22 : if (VSIFWriteL(abyHeader, 256, 1, fp) != 1 ||
904 16 : VSIFSeekL(fp,
905 16 : static_cast<vsi_l_offset>(GDALGetDataTypeSizeBytes(eType)) *
906 16 : nXSize * nYSize -
907 : 1,
908 38 : SEEK_CUR) != 0 ||
909 16 : VSIFWriteL(abyHeader + 255, 1, 1, fp) != 1)
910 : {
911 10 : CPLError(CE_Failure, CPLE_FileIO,
912 : "Failed to extent file to its full size, out of disk space?");
913 :
914 10 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
915 10 : VSIUnlink(pszFilename);
916 10 : return nullptr;
917 : }
918 :
919 12 : if (VSIFCloseL(fp) != 0)
920 : {
921 0 : CPLError(CE_Failure, CPLE_FileIO,
922 : "Failed to extent file to its full size, out of disk space?");
923 0 : VSIUnlink(pszFilename);
924 0 : return nullptr;
925 : }
926 :
927 12 : return GDALDataset::Open(pszFilename, GDAL_OF_RASTER | GDAL_OF_UPDATE);
928 : }
929 :
930 : /************************************************************************/
931 : /* GDALRegister_BT() */
932 : /************************************************************************/
933 :
934 2024 : void GDALRegister_BT()
935 :
936 : {
937 2024 : if (GDALGetDriverByName("BT") != nullptr)
938 283 : return;
939 :
940 1741 : GDALDriver *poDriver = new GDALDriver();
941 :
942 1741 : poDriver->SetDescription("BT");
943 1741 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
944 1741 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
945 1741 : "VTP .bt (Binary Terrain) 1.3 Format");
946 1741 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/bt.html");
947 1741 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "bt");
948 1741 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
949 1741 : "Int16 Int32 Float32");
950 1741 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
951 :
952 1741 : poDriver->pfnOpen = BTDataset::Open;
953 1741 : poDriver->pfnCreate = BTDataset::Create;
954 :
955 1741 : GetGDALDriverManager()->RegisterDriver(poDriver);
956 : }
|
