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