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 : const std::string osPrjFile =
536 20 : CPLResetExtensionSafe(GetDescription(), "prj");
537 10 : VSILFILE *fp = VSIFOpenL(osPrjFile.c_str(), "wt");
538 10 : if (fp != nullptr)
539 : {
540 10 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fp, "%s\n", pszProjection));
541 10 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
542 10 : abyHeader[60] = 1;
543 : }
544 : else
545 : {
546 0 : CPLError(CE_Failure, CPLE_AppDefined,
547 : "Unable to write out .prj file.");
548 0 : eErr = CE_Failure;
549 : }
550 10 : CPLFree(pszProjection);
551 : }
552 :
553 10 : return eErr;
554 : }
555 :
556 : /************************************************************************/
557 : /* Open() */
558 : /************************************************************************/
559 :
560 30991 : GDALDataset *BTDataset::Open(GDALOpenInfo *poOpenInfo)
561 :
562 : {
563 : /* -------------------------------------------------------------------- */
564 : /* Verify that this is some form of binterr file. */
565 : /* -------------------------------------------------------------------- */
566 30991 : if (poOpenInfo->nHeaderBytes < 256 || poOpenInfo->fpL == nullptr)
567 28148 : return nullptr;
568 :
569 2843 : if (!STARTS_WITH((const char *)poOpenInfo->pabyHeader, "binterr"))
570 2816 : return nullptr;
571 :
572 : /* -------------------------------------------------------------------- */
573 : /* Create the dataset. */
574 : /* -------------------------------------------------------------------- */
575 27 : BTDataset *poDS = new BTDataset();
576 :
577 27 : memcpy(poDS->abyHeader, poOpenInfo->pabyHeader, 256);
578 :
579 : /* -------------------------------------------------------------------- */
580 : /* Get the version. */
581 : /* -------------------------------------------------------------------- */
582 27 : char szVersion[4] = {};
583 :
584 27 : strncpy(szVersion, reinterpret_cast<char *>(poDS->abyHeader + 7), 3);
585 27 : szVersion[3] = '\0';
586 27 : poDS->nVersionCode = static_cast<int>(CPLAtof(szVersion) * 10);
587 :
588 : /* -------------------------------------------------------------------- */
589 : /* Extract core header information, being careful about the */
590 : /* version. */
591 : /* -------------------------------------------------------------------- */
592 :
593 27 : GInt32 nIntTemp = 0;
594 27 : memcpy(&nIntTemp, poDS->abyHeader + 10, 4);
595 27 : poDS->nRasterXSize = CPL_LSBWORD32(nIntTemp);
596 :
597 27 : memcpy(&nIntTemp, poDS->abyHeader + 14, 4);
598 27 : poDS->nRasterYSize = CPL_LSBWORD32(nIntTemp);
599 :
600 27 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
601 : {
602 0 : delete poDS;
603 0 : return nullptr;
604 : }
605 :
606 27 : GInt16 nDataSize = 0;
607 27 : memcpy(&nDataSize, poDS->abyHeader + 18, 2);
608 27 : CPL_LSBPTR16(&nDataSize);
609 :
610 27 : GDALDataType eType = GDT_Unknown;
611 27 : if (poDS->abyHeader[20] != 0 && nDataSize == 4)
612 15 : eType = GDT_Float32;
613 12 : else if (poDS->abyHeader[20] == 0 && nDataSize == 4)
614 6 : eType = GDT_Int32;
615 6 : else if (poDS->abyHeader[20] == 0 && nDataSize == 2)
616 6 : eType = GDT_Int16;
617 : else
618 : {
619 0 : CPLError(CE_Failure, CPLE_AppDefined,
620 : ".bt file data type unknown, got datasize=%d.", nDataSize);
621 0 : delete poDS;
622 0 : return nullptr;
623 : }
624 :
625 : /*
626 : rcg, apr 7/06: read offset 62 for vert. units.
627 : If zero, assume 1.0 as per spec.
628 :
629 : */
630 27 : memcpy(&poDS->m_fVscale, poDS->abyHeader + 62, 4);
631 27 : CPL_LSBPTR32(&poDS->m_fVscale);
632 27 : if (poDS->m_fVscale == 0.0f)
633 0 : poDS->m_fVscale = 1.0f;
634 :
635 : /* -------------------------------------------------------------------- */
636 : /* Try to read a .prj file if it is indicated. */
637 : /* -------------------------------------------------------------------- */
638 27 : OGRSpatialReference oSRS;
639 27 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
640 :
641 27 : if (poDS->nVersionCode >= 12 && poDS->abyHeader[60] != 0)
642 : {
643 : const std::string osPrjFile =
644 22 : CPLResetExtensionSafe(poOpenInfo->pszFilename, "prj");
645 11 : VSILFILE *fp = VSIFOpenL(osPrjFile.c_str(), "rt");
646 11 : if (fp != nullptr)
647 : {
648 11 : const int nBufMax = 10000;
649 :
650 11 : char *pszBuffer = static_cast<char *>(CPLMalloc(nBufMax));
651 : const int nBytes =
652 11 : static_cast<int>(VSIFReadL(pszBuffer, 1, nBufMax - 1, fp));
653 11 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
654 :
655 11 : pszBuffer[nBytes] = '\0';
656 :
657 11 : if (oSRS.importFromWkt(pszBuffer) != OGRERR_NONE)
658 : {
659 0 : CPLError(CE_Warning, CPLE_AppDefined,
660 : "Unable to parse .prj file, "
661 : "coordinate system missing.");
662 : }
663 11 : CPLFree(pszBuffer);
664 : }
665 : }
666 :
667 : /* -------------------------------------------------------------------- */
668 : /* If we didn't find a .prj file, try to use internal info. */
669 : /* -------------------------------------------------------------------- */
670 27 : if (oSRS.GetRoot() == nullptr)
671 : {
672 16 : GInt16 nUTMZone = 0;
673 16 : memcpy(&nUTMZone, poDS->abyHeader + 24, 2);
674 16 : CPL_LSBPTR16(&nUTMZone);
675 :
676 16 : GInt16 nDatum = 0;
677 16 : memcpy(&nDatum, poDS->abyHeader + 26, 2);
678 16 : CPL_LSBPTR16(&nDatum);
679 :
680 16 : GInt16 nHUnits = 0;
681 16 : memcpy(&nHUnits, poDS->abyHeader + 22, 2);
682 16 : CPL_LSBPTR16(&nHUnits);
683 :
684 16 : if (nUTMZone != 0)
685 0 : oSRS.SetUTM(std::abs(static_cast<int>(nUTMZone)), nUTMZone > 0);
686 16 : else if (nHUnits != 0)
687 16 : oSRS.SetLocalCS("Unknown");
688 :
689 16 : if (nHUnits == 1)
690 16 : oSRS.SetLinearUnits(SRS_UL_METER, 1.0);
691 0 : else if (nHUnits == 2)
692 0 : oSRS.SetLinearUnits(SRS_UL_FOOT, CPLAtof(SRS_UL_FOOT_CONV));
693 0 : else if (nHUnits == 3)
694 0 : oSRS.SetLinearUnits(SRS_UL_US_FOOT, CPLAtof(SRS_UL_US_FOOT_CONV));
695 :
696 : // Translate some of the more obvious old USGS datum codes
697 16 : if (nDatum == 0)
698 0 : nDatum = 6201;
699 16 : else if (nDatum == 1)
700 0 : nDatum = 6209;
701 16 : else if (nDatum == 2)
702 0 : nDatum = 6210;
703 16 : else if (nDatum == 3)
704 0 : nDatum = 6202;
705 16 : else if (nDatum == 4)
706 0 : nDatum = 6203;
707 16 : else if (nDatum == 6)
708 0 : nDatum = 6222;
709 16 : else if (nDatum == 7)
710 0 : nDatum = 6230;
711 16 : else if (nDatum == 13)
712 0 : nDatum = 6267;
713 16 : else if (nDatum == 14)
714 0 : nDatum = 6269;
715 16 : else if (nDatum == 17)
716 0 : nDatum = 6277;
717 16 : else if (nDatum == 19)
718 0 : nDatum = 6284;
719 16 : else if (nDatum == 21)
720 0 : nDatum = 6301;
721 16 : else if (nDatum == 22)
722 0 : nDatum = 6322;
723 16 : else if (nDatum == 23)
724 0 : nDatum = 6326;
725 :
726 16 : if (!oSRS.IsLocal())
727 : {
728 0 : if (nDatum >= 6000)
729 : {
730 : char szName[32];
731 0 : snprintf(szName, sizeof(szName), "EPSG:%d", nDatum - 2000);
732 0 : oSRS.SetWellKnownGeogCS(szName);
733 : }
734 : else
735 0 : oSRS.SetWellKnownGeogCS("WGS84");
736 : }
737 : }
738 :
739 : /* -------------------------------------------------------------------- */
740 : /* Convert coordinate system back to WKT. */
741 : /* -------------------------------------------------------------------- */
742 27 : if (oSRS.GetRoot() != nullptr)
743 27 : poDS->m_oSRS = std::move(oSRS);
744 :
745 : /* -------------------------------------------------------------------- */
746 : /* Get georeferencing bounds. */
747 : /* -------------------------------------------------------------------- */
748 27 : if (poDS->nVersionCode >= 11)
749 : {
750 27 : double dfLeft = 0.0;
751 27 : memcpy(&dfLeft, poDS->abyHeader + 28, 8);
752 27 : CPL_LSBPTR64(&dfLeft);
753 :
754 27 : double dfRight = 0.0;
755 27 : memcpy(&dfRight, poDS->abyHeader + 36, 8);
756 27 : CPL_LSBPTR64(&dfRight);
757 :
758 27 : double dfBottom = 0.0;
759 27 : memcpy(&dfBottom, poDS->abyHeader + 44, 8);
760 27 : CPL_LSBPTR64(&dfBottom);
761 :
762 27 : double dfTop = 0.0;
763 27 : memcpy(&dfTop, poDS->abyHeader + 52, 8);
764 27 : CPL_LSBPTR64(&dfTop);
765 :
766 27 : poDS->adfGeoTransform[0] = dfLeft;
767 27 : poDS->adfGeoTransform[1] = (dfRight - dfLeft) / poDS->nRasterXSize;
768 27 : poDS->adfGeoTransform[2] = 0.0;
769 27 : poDS->adfGeoTransform[3] = dfTop;
770 27 : poDS->adfGeoTransform[4] = 0.0;
771 27 : poDS->adfGeoTransform[5] = (dfBottom - dfTop) / poDS->nRasterYSize;
772 :
773 27 : poDS->bGeoTransformValid = TRUE;
774 : }
775 :
776 27 : poDS->eAccess = poOpenInfo->eAccess;
777 27 : poDS->fpImage = poOpenInfo->fpL;
778 27 : poOpenInfo->fpL = nullptr;
779 :
780 : /* -------------------------------------------------------------------- */
781 : /* Create band information objects */
782 : /* -------------------------------------------------------------------- */
783 27 : poDS->SetBand(1, new BTRasterBand(poDS, poDS->fpImage, eType));
784 :
785 : #ifdef notdef
786 : poDS->bGeoTransformValid = GDALReadWorldFile(poOpenInfo->pszFilename,
787 : ".wld", poDS->adfGeoTransform);
788 : #endif
789 :
790 : /* -------------------------------------------------------------------- */
791 : /* Initialize any PAM information. */
792 : /* -------------------------------------------------------------------- */
793 27 : poDS->SetDescription(poOpenInfo->pszFilename);
794 27 : poDS->TryLoadXML();
795 :
796 : /* -------------------------------------------------------------------- */
797 : /* Check for overviews. */
798 : /* -------------------------------------------------------------------- */
799 27 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
800 :
801 27 : return poDS;
802 : }
803 :
804 : /************************************************************************/
805 : /* Create() */
806 : /************************************************************************/
807 :
808 66 : GDALDataset *BTDataset::Create(const char *pszFilename, int nXSize, int nYSize,
809 : int nBandsIn, GDALDataType eType,
810 : CPL_UNUSED char **papszOptions)
811 : {
812 :
813 : /* -------------------------------------------------------------------- */
814 : /* Verify input options. */
815 : /* -------------------------------------------------------------------- */
816 66 : if (eType != GDT_Int16 && eType != GDT_Int32 && eType != GDT_Float32)
817 : {
818 38 : CPLError(CE_Failure, CPLE_AppDefined,
819 : "Attempt to create .bt dataset with an illegal "
820 : "data type (%s), only Int16, Int32 and Float32 supported.",
821 : GDALGetDataTypeName(eType));
822 :
823 38 : return nullptr;
824 : }
825 :
826 28 : if (nBandsIn != 1)
827 : {
828 3 : CPLError(
829 : CE_Failure, CPLE_AppDefined,
830 : "Attempt to create .bt dataset with %d bands, only 1 supported",
831 : nBandsIn);
832 :
833 3 : return nullptr;
834 : }
835 :
836 : /* -------------------------------------------------------------------- */
837 : /* Try to create the file. */
838 : /* -------------------------------------------------------------------- */
839 25 : VSILFILE *fp = VSIFOpenL(pszFilename, "wb");
840 :
841 25 : if (fp == nullptr)
842 : {
843 3 : CPLError(CE_Failure, CPLE_OpenFailed,
844 : "Attempt to create file `%s' failed.", pszFilename);
845 3 : return nullptr;
846 : }
847 :
848 : /* -------------------------------------------------------------------- */
849 : /* Setup base header. */
850 : /* -------------------------------------------------------------------- */
851 22 : unsigned char abyHeader[256] = {};
852 :
853 22 : memcpy(abyHeader, "binterr1.3", 10);
854 :
855 22 : GInt32 nTemp = CPL_LSBWORD32(nXSize);
856 22 : memcpy(abyHeader + 10, &nTemp, 4);
857 :
858 22 : nTemp = CPL_LSBWORD32(nYSize);
859 22 : memcpy(abyHeader + 14, &nTemp, 4);
860 :
861 : GInt16 nShortTemp = static_cast<GInt16>(
862 22 : CPL_LSBWORD16((GInt16)(GDALGetDataTypeSize(eType) / 8)));
863 22 : memcpy(abyHeader + 18, &nShortTemp, 2);
864 :
865 22 : if (eType == GDT_Float32)
866 6 : abyHeader[20] = 1;
867 : else
868 16 : abyHeader[20] = 0;
869 :
870 22 : nShortTemp = CPL_LSBWORD16(1); /* meters */
871 22 : memcpy(abyHeader + 22, &nShortTemp, 2);
872 :
873 22 : nShortTemp = CPL_LSBWORD16(0); /* not utm */
874 22 : memcpy(abyHeader + 24, &nShortTemp, 2);
875 :
876 22 : nShortTemp = CPL_LSBWORD16(-2); /* datum unknown */
877 22 : memcpy(abyHeader + 26, &nShortTemp, 2);
878 :
879 : /* -------------------------------------------------------------------- */
880 : /* Set dummy extents. */
881 : /* -------------------------------------------------------------------- */
882 22 : double dfLeft = 0.0;
883 22 : double dfRight = nXSize;
884 22 : double dfTop = nYSize;
885 22 : double dfBottom = 0.0;
886 :
887 22 : memcpy(abyHeader + 28, &dfLeft, 8);
888 22 : memcpy(abyHeader + 36, &dfRight, 8);
889 22 : memcpy(abyHeader + 44, &dfBottom, 8);
890 22 : memcpy(abyHeader + 52, &dfTop, 8);
891 :
892 22 : CPL_LSBPTR64(abyHeader + 28);
893 22 : CPL_LSBPTR64(abyHeader + 36);
894 22 : CPL_LSBPTR64(abyHeader + 44);
895 22 : CPL_LSBPTR64(abyHeader + 52);
896 :
897 : /* -------------------------------------------------------------------- */
898 : /* Set dummy scale. */
899 : /* -------------------------------------------------------------------- */
900 22 : float fScale = 1.0;
901 :
902 22 : memcpy(abyHeader + 62, &fScale, 4);
903 22 : CPL_LSBPTR32(abyHeader + 62);
904 :
905 : /* -------------------------------------------------------------------- */
906 : /* Write to disk. */
907 : /* -------------------------------------------------------------------- */
908 22 : if (VSIFWriteL(abyHeader, 256, 1, fp) != 1 ||
909 16 : VSIFSeekL(fp,
910 16 : static_cast<vsi_l_offset>(GDALGetDataTypeSizeBytes(eType)) *
911 16 : nXSize * nYSize -
912 : 1,
913 38 : SEEK_CUR) != 0 ||
914 16 : VSIFWriteL(abyHeader + 255, 1, 1, fp) != 1)
915 : {
916 10 : CPLError(CE_Failure, CPLE_FileIO,
917 : "Failed to extent file to its full size, out of disk space?");
918 :
919 10 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
920 10 : VSIUnlink(pszFilename);
921 10 : return nullptr;
922 : }
923 :
924 12 : if (VSIFCloseL(fp) != 0)
925 : {
926 0 : CPLError(CE_Failure, CPLE_FileIO,
927 : "Failed to extent file to its full size, out of disk space?");
928 0 : VSIUnlink(pszFilename);
929 0 : return nullptr;
930 : }
931 :
932 12 : return GDALDataset::Open(pszFilename, GDAL_OF_RASTER | GDAL_OF_UPDATE);
933 : }
934 :
935 : /************************************************************************/
936 : /* GDALRegister_BT() */
937 : /************************************************************************/
938 :
939 1682 : void GDALRegister_BT()
940 :
941 : {
942 1682 : if (GDALGetDriverByName("BT") != nullptr)
943 301 : return;
944 :
945 1381 : GDALDriver *poDriver = new GDALDriver();
946 :
947 1381 : poDriver->SetDescription("BT");
948 1381 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
949 1381 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
950 1381 : "VTP .bt (Binary Terrain) 1.3 Format");
951 1381 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/bt.html");
952 1381 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "bt");
953 1381 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
954 1381 : "Int16 Int32 Float32");
955 1381 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
956 :
957 1381 : poDriver->pfnOpen = BTDataset::Open;
958 1381 : poDriver->pfnCreate = BTDataset::Create;
959 :
960 1381 : GetGDALDriverManager()->RegisterDriver(poDriver);
961 : }
|
