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