Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: USGS DEM Driver
4 : * Purpose: All reader for USGS DEM Reader
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : * Portions of this module derived from the VTP USGS DEM driver by Ben
8 : * Discoe, see http://www.vterrain.org
9 : *
10 : ******************************************************************************
11 : * Copyright (c) 2001, Frank Warmerdam <warmerdam@pobox.com>
12 : * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
13 : *
14 : * SPDX-License-Identifier: MIT
15 : ****************************************************************************/
16 :
17 : #include "gdal_frmts.h"
18 : #include "gdal_pam.h"
19 : #include "gdal_driver.h"
20 : #include "gdal_drivermanager.h"
21 : #include "gdal_openinfo.h"
22 : #include "gdal_cpp_functions.h"
23 : #include "ogr_spatialref.h"
24 :
25 : #include <algorithm>
26 : #include <cmath>
27 :
28 : typedef struct
29 : {
30 : double x;
31 : double y;
32 : } DPoint2;
33 :
34 : constexpr int USGSDEM_NODATA = -32767;
35 :
36 : GDALDataset *USGSDEMCreateCopy(const char *, GDALDataset *, int, char **,
37 : GDALProgressFunc pfnProgress,
38 : void *pProgressData);
39 :
40 : /************************************************************************/
41 : /* ReadInt() */
42 : /************************************************************************/
43 :
44 222 : static int ReadInt(VSILFILE *fp)
45 : {
46 : char c;
47 222 : int nRead = 0;
48 : char szBuffer[12];
49 222 : bool bInProlog = true;
50 :
51 : while (true)
52 : {
53 2266 : if (VSIFReadL(&c, 1, 1, fp) != 1)
54 : {
55 4 : return 0;
56 : }
57 2262 : if (bInProlog)
58 : {
59 1972 : if (!isspace(static_cast<unsigned char>(c)))
60 : {
61 218 : bInProlog = false;
62 : }
63 : }
64 2262 : if (!bInProlog)
65 : {
66 508 : if (c != '-' && c != '+' && !(c >= '0' && c <= '9'))
67 : {
68 218 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, VSIFTellL(fp) - 1, SEEK_SET));
69 218 : break;
70 : }
71 290 : if (nRead < 11)
72 290 : szBuffer[nRead] = c;
73 290 : nRead++;
74 : }
75 : }
76 218 : szBuffer[std::min(nRead, 11)] = 0;
77 218 : return atoi(szBuffer);
78 : }
79 :
80 : typedef struct
81 : {
82 : VSILFILE *fp;
83 : int max_size;
84 : char *buffer;
85 : int buffer_size;
86 : int cur_index;
87 : } Buffer;
88 :
89 : /************************************************************************/
90 : /* USGSDEMRefillBuffer() */
91 : /************************************************************************/
92 :
93 42 : static void USGSDEMRefillBuffer(Buffer *psBuffer)
94 : {
95 42 : memmove(psBuffer->buffer, psBuffer->buffer + psBuffer->cur_index,
96 42 : psBuffer->buffer_size - psBuffer->cur_index);
97 :
98 42 : psBuffer->buffer_size -= psBuffer->cur_index;
99 42 : psBuffer->buffer_size += static_cast<int>(
100 84 : VSIFReadL(psBuffer->buffer + psBuffer->buffer_size, 1,
101 42 : psBuffer->max_size - psBuffer->buffer_size, psBuffer->fp));
102 42 : psBuffer->cur_index = 0;
103 42 : }
104 :
105 : /************************************************************************/
106 : /* USGSDEMGetCurrentFilePos() */
107 : /************************************************************************/
108 :
109 10 : static vsi_l_offset USGSDEMGetCurrentFilePos(const Buffer *psBuffer)
110 : {
111 10 : return VSIFTellL(psBuffer->fp) - psBuffer->buffer_size +
112 10 : psBuffer->cur_index;
113 : }
114 :
115 : /************************************************************************/
116 : /* USGSDEMSetCurrentFilePos() */
117 : /************************************************************************/
118 :
119 10 : static void USGSDEMSetCurrentFilePos(Buffer *psBuffer, vsi_l_offset nNewPos)
120 : {
121 10 : vsi_l_offset nCurPosFP = VSIFTellL(psBuffer->fp);
122 10 : if (nNewPos >= nCurPosFP - psBuffer->buffer_size && nNewPos < nCurPosFP)
123 : {
124 5 : psBuffer->cur_index =
125 5 : static_cast<int>(nNewPos - (nCurPosFP - psBuffer->buffer_size));
126 : }
127 : else
128 : {
129 5 : CPL_IGNORE_RET_VAL(VSIFSeekL(psBuffer->fp, nNewPos, SEEK_SET));
130 5 : psBuffer->buffer_size = 0;
131 5 : psBuffer->cur_index = 0;
132 : }
133 10 : }
134 :
135 : /************************************************************************/
136 : /* USGSDEMReadIntFromBuffer() */
137 : /************************************************************************/
138 :
139 942758 : static int USGSDEMReadIntFromBuffer(Buffer *psBuffer, int *pbSuccess = nullptr)
140 : {
141 : char c;
142 :
143 : while (true)
144 : {
145 942758 : if (psBuffer->cur_index >= psBuffer->buffer_size)
146 : {
147 38 : USGSDEMRefillBuffer(psBuffer);
148 38 : if (psBuffer->cur_index >= psBuffer->buffer_size)
149 : {
150 0 : if (pbSuccess)
151 0 : *pbSuccess = FALSE;
152 0 : return 0;
153 : }
154 : }
155 :
156 942758 : c = psBuffer->buffer[psBuffer->cur_index];
157 942758 : psBuffer->cur_index++;
158 942758 : if (!isspace(static_cast<unsigned char>(c)))
159 14451 : break;
160 : }
161 :
162 14451 : GIntBig nVal = 0;
163 14451 : int nSign = 1;
164 14451 : if (c == '-')
165 4944 : nSign = -1;
166 9507 : else if (c == '+')
167 0 : nSign = 1;
168 9507 : else if (c >= '0' && c <= '9')
169 9507 : nVal = c - '0';
170 : else
171 : {
172 0 : if (pbSuccess)
173 0 : *pbSuccess = FALSE;
174 0 : return 0;
175 : }
176 :
177 : while (true)
178 : {
179 47517 : if (psBuffer->cur_index >= psBuffer->buffer_size)
180 : {
181 0 : USGSDEMRefillBuffer(psBuffer);
182 0 : if (psBuffer->cur_index >= psBuffer->buffer_size)
183 : {
184 0 : if (pbSuccess)
185 0 : *pbSuccess = TRUE;
186 0 : return static_cast<int>(nSign * nVal);
187 : }
188 : }
189 :
190 47517 : c = psBuffer->buffer[psBuffer->cur_index];
191 47517 : if (c >= '0' && c <= '9')
192 : {
193 33066 : psBuffer->cur_index++;
194 33066 : if (nVal * nSign < INT_MAX && nVal * nSign > INT_MIN)
195 : {
196 33066 : nVal = nVal * 10 + (c - '0');
197 33066 : if (nVal * nSign > INT_MAX)
198 : {
199 0 : nVal = INT_MAX;
200 0 : nSign = 1;
201 : }
202 33066 : else if (nVal * nSign < INT_MIN)
203 : {
204 0 : nVal = INT_MIN;
205 0 : nSign = 1;
206 : }
207 : }
208 : }
209 : else
210 : {
211 14451 : if (pbSuccess)
212 14451 : *pbSuccess = TRUE;
213 14451 : return static_cast<int>(nSign * nVal);
214 : }
215 : }
216 : }
217 :
218 : /************************************************************************/
219 : /* USGSDEMReadDoubleFromBuffer() */
220 : /************************************************************************/
221 :
222 5184 : static double USGSDEMReadDoubleFromBuffer(Buffer *psBuffer, int nCharCount,
223 : int *pbSuccess = nullptr)
224 :
225 : {
226 5184 : if (psBuffer->cur_index + nCharCount > psBuffer->buffer_size)
227 : {
228 4 : USGSDEMRefillBuffer(psBuffer);
229 4 : if (psBuffer->cur_index + nCharCount > psBuffer->buffer_size)
230 : {
231 1 : if (pbSuccess)
232 1 : *pbSuccess = FALSE;
233 1 : return 0;
234 : }
235 : }
236 :
237 5183 : char *szPtr = psBuffer->buffer + psBuffer->cur_index;
238 5183 : char backupC = szPtr[nCharCount];
239 5183 : szPtr[nCharCount] = 0;
240 129575 : for (int i = 0; i < nCharCount; i++)
241 : {
242 124392 : if (szPtr[i] == 'D')
243 5164 : szPtr[i] = 'E';
244 : }
245 :
246 5183 : double dfVal = CPLAtof(szPtr);
247 5183 : szPtr[nCharCount] = backupC;
248 5183 : psBuffer->cur_index += nCharCount;
249 :
250 5183 : if (pbSuccess)
251 5183 : *pbSuccess = TRUE;
252 5183 : return dfVal;
253 : }
254 :
255 : /************************************************************************/
256 : /* DConvert() */
257 : /************************************************************************/
258 :
259 246 : static double DConvert(VSILFILE *fp, int nCharCount)
260 :
261 : {
262 : char szBuffer[100];
263 :
264 246 : CPL_IGNORE_RET_VAL(VSIFReadL(szBuffer, nCharCount, 1, fp));
265 246 : szBuffer[nCharCount] = '\0';
266 :
267 6366 : for (int i = 0; i < nCharCount; i++)
268 : {
269 6120 : if (szBuffer[i] == 'D')
270 178 : szBuffer[i] = 'E';
271 : }
272 :
273 492 : return CPLAtof(szBuffer);
274 : }
275 :
276 : /************************************************************************/
277 : /* ==================================================================== */
278 : /* USGSDEMDataset */
279 : /* ==================================================================== */
280 : /************************************************************************/
281 :
282 : class USGSDEMRasterBand;
283 :
284 : class USGSDEMDataset final : public GDALPamDataset
285 : {
286 : friend class USGSDEMRasterBand;
287 :
288 : int nDataStartOffset;
289 : GDALDataType eNaturalDataFormat;
290 :
291 : GDALGeoTransform m_gt{};
292 : OGRSpatialReference m_oSRS{};
293 :
294 : double fVRes;
295 :
296 : const char *pszUnits;
297 :
298 : int LoadFromFile(VSILFILE *);
299 :
300 : VSILFILE *fp;
301 :
302 : public:
303 : USGSDEMDataset();
304 : ~USGSDEMDataset() override;
305 :
306 : static int Identify(GDALOpenInfo *);
307 : static GDALDataset *Open(GDALOpenInfo *);
308 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
309 : const OGRSpatialReference *GetSpatialRef() const override;
310 : };
311 :
312 : /************************************************************************/
313 : /* ==================================================================== */
314 : /* USGSDEMRasterBand */
315 : /* ==================================================================== */
316 : /************************************************************************/
317 :
318 : class USGSDEMRasterBand final : public GDALPamRasterBand
319 : {
320 : friend class USGSDEMDataset;
321 :
322 : public:
323 : explicit USGSDEMRasterBand(USGSDEMDataset *);
324 :
325 : const char *GetUnitType() override;
326 : double GetNoDataValue(int *pbSuccess = nullptr) override;
327 : CPLErr IReadBlock(int, int, void *) override;
328 : };
329 :
330 : /************************************************************************/
331 : /* USGSDEMRasterBand() */
332 : /************************************************************************/
333 :
334 18 : USGSDEMRasterBand::USGSDEMRasterBand(USGSDEMDataset *poDSIn)
335 :
336 : {
337 18 : this->poDS = poDSIn;
338 18 : this->nBand = 1;
339 :
340 18 : eDataType = poDSIn->eNaturalDataFormat;
341 :
342 18 : nBlockXSize = poDSIn->GetRasterXSize();
343 18 : nBlockYSize = poDSIn->GetRasterYSize();
344 18 : }
345 :
346 : /************************************************************************/
347 : /* IReadBlock() */
348 : /************************************************************************/
349 :
350 9 : CPLErr USGSDEMRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
351 : CPL_UNUSED int nBlockYOff, void *pImage)
352 :
353 : {
354 : /* int bad = FALSE; */
355 9 : USGSDEMDataset *poGDS = cpl::down_cast<USGSDEMDataset *>(poDS);
356 :
357 : /* -------------------------------------------------------------------- */
358 : /* Initialize image buffer to nodata value. */
359 : /* -------------------------------------------------------------------- */
360 9 : GDALCopyWords(&USGSDEM_NODATA, GDT_Int32, 0, pImage, GetRasterDataType(),
361 : GDALGetDataTypeSizeBytes(GetRasterDataType()),
362 9 : GetXSize() * GetYSize());
363 :
364 : /* -------------------------------------------------------------------- */
365 : /* Seek to data. */
366 : /* -------------------------------------------------------------------- */
367 9 : CPL_IGNORE_RET_VAL(VSIFSeekL(poGDS->fp, poGDS->nDataStartOffset, 0));
368 :
369 9 : double dfYMin = poGDS->m_gt[3] + (GetYSize() - 0.5) * poGDS->m_gt[5];
370 :
371 : /* -------------------------------------------------------------------- */
372 : /* Read all the profiles into the image buffer. */
373 : /* -------------------------------------------------------------------- */
374 :
375 : Buffer sBuffer;
376 9 : sBuffer.max_size = 32768;
377 9 : sBuffer.buffer = static_cast<char *>(CPLMalloc(sBuffer.max_size + 1));
378 9 : sBuffer.fp = poGDS->fp;
379 9 : sBuffer.buffer_size = 0;
380 9 : sBuffer.cur_index = 0;
381 :
382 1045 : for (int i = 0; i < GetXSize(); i++)
383 : {
384 : int bSuccess;
385 1037 : const int nRowNumber = USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
386 1037 : if (nRowNumber != 1)
387 1 : CPLDebug("USGSDEM", "i = %d, nRowNumber = %d", i, nRowNumber);
388 1037 : if (bSuccess)
389 : {
390 : const int nColNumber =
391 1037 : USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
392 1037 : if (nColNumber != i + 1)
393 : {
394 3 : CPLDebug("USGSDEM", "i = %d, nColNumber = %d", i, nColNumber);
395 : }
396 : }
397 : const int nCPoints =
398 1037 : (bSuccess) ? USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess) : 0;
399 : #ifdef DEBUG_VERBOSE
400 : CPLDebug("USGSDEM", "i = %d, nCPoints = %d", i, nCPoints);
401 : #endif
402 :
403 1037 : if (bSuccess)
404 : {
405 : const int nNumberOfCols =
406 1037 : USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
407 1037 : if (nNumberOfCols != 1)
408 : {
409 0 : CPLDebug("USGSDEM", "i = %d, nNumberOfCols = %d", i,
410 : nNumberOfCols);
411 : }
412 : }
413 :
414 : // x-start
415 1037 : if (bSuccess)
416 1037 : /* dxStart = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24,
417 : &bSuccess);
418 :
419 : double dyStart =
420 1037 : (bSuccess) ? USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess)
421 1037 : : 0;
422 : const double dfElevOffset =
423 1037 : (bSuccess) ? USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess)
424 1037 : : 0;
425 :
426 : // min z value
427 1037 : if (bSuccess)
428 1037 : /* djunk = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess);
429 :
430 : // max z value
431 1037 : if (bSuccess)
432 1036 : /* djunk = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess);
433 1037 : if (!bSuccess)
434 : {
435 1 : CPLFree(sBuffer.buffer);
436 1 : return CE_Failure;
437 : }
438 :
439 1036 : if (poGDS->m_oSRS.IsGeographic())
440 4 : dyStart = dyStart / 3600.0;
441 :
442 1036 : double dygap = (dfYMin - dyStart) / poGDS->m_gt[5] + 0.5;
443 1036 : if (dygap <= INT_MIN || dygap >= INT_MAX || !std::isfinite(dygap))
444 : {
445 0 : CPLFree(sBuffer.buffer);
446 0 : return CE_Failure;
447 : }
448 1036 : int lygap = static_cast<int>(dygap);
449 1036 : if (nCPoints <= 0)
450 0 : continue;
451 1036 : if (lygap > INT_MAX - nCPoints)
452 0 : lygap = INT_MAX - nCPoints;
453 1036 : if (lygap < 0 && GetYSize() > INT_MAX + lygap)
454 : {
455 0 : CPLFree(sBuffer.buffer);
456 0 : return CE_Failure;
457 : }
458 :
459 11339 : for (int j = lygap; j < (nCPoints + lygap); j++)
460 : {
461 10303 : const int iY = GetYSize() - j - 1;
462 :
463 10303 : const int nElev = USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
464 : #ifdef DEBUG_VERBOSE
465 : CPLDebug("USGSDEM", " j - lygap = %d, nElev = %d", j - lygap,
466 : nElev);
467 : #endif
468 :
469 10303 : if (!bSuccess)
470 : {
471 0 : CPLFree(sBuffer.buffer);
472 0 : return CE_Failure;
473 : }
474 :
475 10303 : if (iY < 0 || iY >= GetYSize())
476 : {
477 : /* bad = TRUE; */
478 : }
479 10303 : else if (nElev == USGSDEM_NODATA)
480 : /* leave in output buffer as nodata */;
481 : else
482 : {
483 6341 : const float fComputedElev =
484 6341 : static_cast<float>(nElev * poGDS->fVRes + dfElevOffset);
485 :
486 6341 : if (GetRasterDataType() == GDT_Int16)
487 : {
488 12560 : GUInt16 nVal = (fComputedElev < -32768) ? -32768
489 : : (fComputedElev > 32767)
490 : ? 32767
491 6280 : : static_cast<GInt16>(fComputedElev);
492 6280 : reinterpret_cast<GInt16 *>(pImage)[i + iY * GetXSize()] =
493 6280 : nVal;
494 : }
495 : else
496 : {
497 61 : reinterpret_cast<float *>(pImage)[i + iY * GetXSize()] =
498 : fComputedElev;
499 : }
500 : }
501 : }
502 :
503 1036 : if (poGDS->nDataStartOffset == 1024)
504 : {
505 : // Seek to the next 1024 byte boundary.
506 : // Some files have 'junk' profile values after the valid/declared
507 : // ones
508 10 : vsi_l_offset nCurPos = USGSDEMGetCurrentFilePos(&sBuffer);
509 10 : vsi_l_offset nNewPos = (nCurPos + 1023) / 1024 * 1024;
510 10 : if (nNewPos > nCurPos)
511 : {
512 10 : USGSDEMSetCurrentFilePos(&sBuffer, nNewPos);
513 : }
514 : }
515 : }
516 8 : CPLFree(sBuffer.buffer);
517 :
518 8 : return CE_None;
519 : }
520 :
521 : /************************************************************************/
522 : /* GetNoDataValue() */
523 : /************************************************************************/
524 :
525 0 : double USGSDEMRasterBand::GetNoDataValue(int *pbSuccess)
526 :
527 : {
528 0 : if (pbSuccess != nullptr)
529 0 : *pbSuccess = TRUE;
530 :
531 0 : return USGSDEM_NODATA;
532 : }
533 :
534 : /************************************************************************/
535 : /* GetUnitType() */
536 : /************************************************************************/
537 0 : const char *USGSDEMRasterBand::GetUnitType()
538 : {
539 0 : USGSDEMDataset *poGDS = cpl::down_cast<USGSDEMDataset *>(poDS);
540 :
541 0 : return poGDS->pszUnits;
542 : }
543 :
544 : /************************************************************************/
545 : /* ==================================================================== */
546 : /* USGSDEMDataset */
547 : /* ==================================================================== */
548 : /************************************************************************/
549 :
550 : /************************************************************************/
551 : /* USGSDEMDataset() */
552 : /************************************************************************/
553 :
554 18 : USGSDEMDataset::USGSDEMDataset()
555 : : nDataStartOffset(0), eNaturalDataFormat(GDT_Unknown), fVRes(0.0),
556 18 : pszUnits(nullptr), fp(nullptr)
557 : {
558 18 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
559 18 : }
560 :
561 : /************************************************************************/
562 : /* ~USGSDEMDataset() */
563 : /************************************************************************/
564 :
565 36 : USGSDEMDataset::~USGSDEMDataset()
566 :
567 : {
568 18 : FlushCache(true);
569 :
570 18 : if (fp != nullptr)
571 18 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
572 36 : }
573 :
574 : /************************************************************************/
575 : /* LoadFromFile() */
576 : /* */
577 : /* If the data from DEM is in meters, then values are stored as */
578 : /* shorts. If DEM data is in feet, then height data will be */
579 : /* stored in float, to preserve the precision of the original */
580 : /* data. returns true if the file was successfully opened and */
581 : /* read. */
582 : /************************************************************************/
583 :
584 18 : int USGSDEMDataset::LoadFromFile(VSILFILE *InDem)
585 : {
586 : // check for version of DEM format
587 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 864, 0));
588 :
589 : // Read DEM into matrix
590 18 : const int nRow = ReadInt(InDem);
591 18 : const int nColumn = ReadInt(InDem);
592 : const bool bNewFormat =
593 18 : VSIFTellL(InDem) >= 1024 || nRow != 1 || nColumn != 1;
594 18 : if (bNewFormat)
595 : {
596 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 1024, 0)); // New Format
597 18 : int i = ReadInt(InDem);
598 18 : int j = ReadInt(InDem);
599 18 : if (i != 1 || (j != 1 && j != 0)) // File OK?
600 : {
601 4 : CPL_IGNORE_RET_VAL(
602 4 : VSIFSeekL(InDem, 893, 0)); // Undocumented Format (39109h1.dem)
603 4 : i = ReadInt(InDem);
604 4 : j = ReadInt(InDem);
605 4 : if (i != 1 || j != 1) // File OK?
606 : {
607 2 : CPL_IGNORE_RET_VAL(VSIFSeekL(
608 : InDem, 918, 0)); // Latest iteration of the A record, such
609 : // as in fema06-140cm_2995441b.dem
610 2 : i = ReadInt(InDem);
611 2 : j = ReadInt(InDem);
612 2 : if (i != 1 || j != 1) // File OK?
613 : {
614 0 : CPLError(CE_Failure, CPLE_AppDefined,
615 : "Does not appear to be a USGS DEM file.");
616 0 : return FALSE;
617 : }
618 : else
619 2 : nDataStartOffset = 918;
620 : }
621 : else
622 2 : nDataStartOffset = 893;
623 : }
624 : else
625 : {
626 14 : nDataStartOffset = 1024;
627 :
628 : // Some files use 1025 byte records ending with a newline character.
629 : // See https://github.com/OSGeo/gdal/issues/5007
630 14 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 1024, 0));
631 : char c;
632 28 : if (VSIFReadL(&c, 1, 1, InDem) == 1 && c == '\n' &&
633 2 : VSIFSeekL(InDem, 1024 + 1024 + 1, 0) == 0 &&
634 28 : VSIFReadL(&c, 1, 1, InDem) == 1 && c == '\n')
635 : {
636 2 : nDataStartOffset = 1025;
637 : }
638 : }
639 : }
640 : else
641 0 : nDataStartOffset = 864;
642 :
643 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 156, 0));
644 18 : const int nCoordSystem = ReadInt(InDem);
645 18 : const int iUTMZone = ReadInt(InDem);
646 :
647 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 528, 0));
648 18 : const int nGUnit = ReadInt(InDem);
649 18 : const int nVUnit = ReadInt(InDem);
650 :
651 : // Vertical Units in meters
652 18 : if (nVUnit == 1)
653 0 : pszUnits = "ft";
654 : else
655 18 : pszUnits = "m";
656 :
657 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 816, 0));
658 18 : const double dxdelta = DConvert(InDem, 12);
659 18 : const double dydelta = DConvert(InDem, 12);
660 18 : if (dydelta == 0)
661 0 : return FALSE;
662 18 : fVRes = DConvert(InDem, 12);
663 :
664 : /* -------------------------------------------------------------------- */
665 : /* Should we treat this as floating point, or GInt16. */
666 : /* -------------------------------------------------------------------- */
667 18 : if (nVUnit == 1 || fVRes < 1.0)
668 4 : eNaturalDataFormat = GDT_Float32;
669 : else
670 14 : eNaturalDataFormat = GDT_Int16;
671 :
672 : /* -------------------------------------------------------------------- */
673 : /* Read four corner coordinates. */
674 : /* -------------------------------------------------------------------- */
675 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 546, 0));
676 : DPoint2 corners[4]; // SW, NW, NE, SE
677 90 : for (int i = 0; i < 4; i++)
678 : {
679 72 : corners[i].x = DConvert(InDem, 24);
680 72 : corners[i].y = DConvert(InDem, 24);
681 : }
682 :
683 : // find absolute extents of raw vales
684 : DPoint2 extent_min, extent_max;
685 18 : extent_min.x = std::min(corners[0].x, corners[1].x);
686 18 : extent_max.x = std::max(corners[2].x, corners[3].x);
687 18 : extent_min.y = std::min(corners[0].y, corners[3].y);
688 18 : extent_max.y = std::max(corners[1].y, corners[2].y);
689 :
690 18 : /* dElevMin = */ DConvert(InDem, 48);
691 18 : /* dElevMax = */ DConvert(InDem, 48);
692 :
693 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 858, 0));
694 18 : const int nProfiles = ReadInt(InDem);
695 :
696 : /* -------------------------------------------------------------------- */
697 : /* Collect the spatial reference system. */
698 : /* -------------------------------------------------------------------- */
699 36 : OGRSpatialReference sr;
700 18 : sr.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
701 18 : bool bNAD83 = true;
702 :
703 : // OLD format header ends at byte 864
704 18 : if (bNewFormat)
705 : {
706 : // year of data compilation
707 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 876, 0));
708 : char szDateBuffer[5];
709 18 : CPL_IGNORE_RET_VAL(VSIFReadL(szDateBuffer, 4, 1, InDem));
710 : /* szDateBuffer[4] = 0; */
711 :
712 : // Horizontal datum
713 : // 1=North American Datum 1927 (NAD 27)
714 : // 2=World Geodetic System 1972 (WGS 72)
715 : // 3=WGS 84
716 : // 4=NAD 83
717 : // 5=Old Hawaii Datum
718 : // 6=Puerto Rico Datum
719 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 890, 0));
720 :
721 : char szHorzDatum[3];
722 18 : CPL_IGNORE_RET_VAL(VSIFReadL(szHorzDatum, 1, 2, InDem));
723 18 : szHorzDatum[2] = '\0';
724 18 : const int datum = atoi(szHorzDatum);
725 18 : switch (datum)
726 : {
727 4 : case 1:
728 4 : sr.SetWellKnownGeogCS("NAD27");
729 4 : bNAD83 = false;
730 4 : break;
731 :
732 2 : case 2:
733 2 : sr.SetWellKnownGeogCS("WGS72");
734 2 : break;
735 :
736 0 : case 3:
737 0 : sr.SetWellKnownGeogCS("WGS84");
738 0 : break;
739 :
740 2 : case 4:
741 2 : sr.SetWellKnownGeogCS("NAD83");
742 2 : break;
743 :
744 0 : case -9:
745 0 : break;
746 :
747 10 : default:
748 10 : sr.SetWellKnownGeogCS("NAD27");
749 10 : break;
750 : }
751 : }
752 : else
753 : {
754 0 : sr.SetWellKnownGeogCS("NAD27");
755 0 : bNAD83 = false;
756 : }
757 :
758 18 : if (nCoordSystem == 1) // UTM
759 : {
760 12 : if (iUTMZone >= -60 && iUTMZone <= 60)
761 : {
762 12 : sr.SetUTM(abs(iUTMZone), iUTMZone >= 0);
763 12 : if (nGUnit == 1)
764 : {
765 0 : sr.SetLinearUnitsAndUpdateParameters(
766 : SRS_UL_US_FOOT, CPLAtof(SRS_UL_US_FOOT_CONV));
767 : char szUTMName[128];
768 0 : snprintf(szUTMName, sizeof(szUTMName),
769 : "UTM Zone %d, Northern Hemisphere, us-ft", iUTMZone);
770 0 : sr.SetNode("PROJCS", szUTMName);
771 : }
772 : }
773 : }
774 6 : else if (nCoordSystem == 2) // state plane
775 : {
776 0 : if (nGUnit == 1)
777 0 : sr.SetStatePlane(iUTMZone, bNAD83, "Foot",
778 : CPLAtof(SRS_UL_US_FOOT_CONV));
779 : else
780 0 : sr.SetStatePlane(iUTMZone, bNAD83);
781 : }
782 :
783 18 : m_oSRS = std::move(sr);
784 :
785 : /* -------------------------------------------------------------------- */
786 : /* For UTM we use the extents (really the UTM coordinates of */
787 : /* the lat/long corners of the quad) to determine the size in */
788 : /* pixels and lines, but we have to make the anchors be modulus */
789 : /* the pixel size which what really gets used. */
790 : /* -------------------------------------------------------------------- */
791 18 : if (nCoordSystem == 1 // UTM
792 6 : || nCoordSystem == 2 // State Plane
793 6 : || nCoordSystem == -9999) // unknown
794 : {
795 : // expand extents modulus the pixel size.
796 12 : extent_min.y = floor(extent_min.y / dydelta) * dydelta;
797 12 : extent_max.y = ceil(extent_max.y / dydelta) * dydelta;
798 :
799 : // Forcibly compute X extents based on first profile and pixelsize.
800 12 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, nDataStartOffset, 0));
801 12 : /* njunk = */ ReadInt(InDem);
802 12 : /* njunk = */ ReadInt(InDem);
803 12 : /* njunk = */ ReadInt(InDem);
804 12 : /* njunk = */ ReadInt(InDem);
805 12 : const double dxStart = DConvert(InDem, 24);
806 :
807 12 : nRasterYSize =
808 12 : static_cast<int>((extent_max.y - extent_min.y) / dydelta + 1.5);
809 12 : nRasterXSize = nProfiles;
810 :
811 12 : m_gt[0] = dxStart - dxdelta / 2.0;
812 12 : m_gt[1] = dxdelta;
813 12 : m_gt[2] = 0.0;
814 12 : m_gt[3] = extent_max.y + dydelta / 2.0;
815 12 : m_gt[4] = 0.0;
816 12 : m_gt[5] = -dydelta;
817 : }
818 : /* -------------------------------------------------------------------- */
819 : /* Geographic -- use corners directly. */
820 : /* -------------------------------------------------------------------- */
821 : else
822 : {
823 6 : nRasterYSize =
824 6 : static_cast<int>((extent_max.y - extent_min.y) / dydelta + 1.5);
825 6 : nRasterXSize = nProfiles;
826 :
827 : // Translate extents from arc-seconds to decimal degrees.
828 6 : m_gt[0] = (extent_min.x - dxdelta / 2.0) / 3600.0;
829 6 : m_gt[1] = dxdelta / 3600.0;
830 6 : m_gt[2] = 0.0;
831 6 : m_gt[3] = (extent_max.y + dydelta / 2.0) / 3600.0;
832 6 : m_gt[4] = 0.0;
833 6 : m_gt[5] = (-dydelta) / 3600.0;
834 : }
835 :
836 : // IReadBlock() not ready for more than INT_MAX pixels, and that
837 : // would behave badly
838 36 : if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize) ||
839 18 : nRasterXSize > INT_MAX / nRasterYSize)
840 : {
841 0 : return FALSE;
842 : }
843 :
844 18 : return TRUE;
845 : }
846 :
847 : /************************************************************************/
848 : /* GetGeoTransform() */
849 : /************************************************************************/
850 :
851 6 : CPLErr USGSDEMDataset::GetGeoTransform(GDALGeoTransform >) const
852 :
853 : {
854 6 : gt = m_gt;
855 6 : return CE_None;
856 : }
857 :
858 : /************************************************************************/
859 : /* GetSpatialRef() */
860 : /************************************************************************/
861 :
862 6 : const OGRSpatialReference *USGSDEMDataset::GetSpatialRef() const
863 :
864 : {
865 6 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
866 : }
867 :
868 : /************************************************************************/
869 : /* Identify() */
870 : /************************************************************************/
871 :
872 58349 : int USGSDEMDataset::Identify(GDALOpenInfo *poOpenInfo)
873 :
874 : {
875 58349 : if (poOpenInfo->nHeaderBytes < 200)
876 55149 : return FALSE;
877 :
878 3200 : if (!STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 0") &&
879 3192 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 1") &&
880 3168 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 2") &&
881 3168 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 3") &&
882 3168 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " -9999"))
883 3168 : return FALSE;
884 :
885 32 : if (!STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 150, " 1") &&
886 4 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 150, " 4"))
887 0 : return FALSE;
888 :
889 32 : return TRUE;
890 : }
891 :
892 : /************************************************************************/
893 : /* Open() */
894 : /************************************************************************/
895 :
896 18 : GDALDataset *USGSDEMDataset::Open(GDALOpenInfo *poOpenInfo)
897 :
898 : {
899 18 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
900 0 : return nullptr;
901 :
902 : /* -------------------------------------------------------------------- */
903 : /* Create a corresponding GDALDataset. */
904 : /* -------------------------------------------------------------------- */
905 18 : USGSDEMDataset *poDS = new USGSDEMDataset();
906 :
907 18 : poDS->fp = poOpenInfo->fpL;
908 18 : poOpenInfo->fpL = nullptr;
909 :
910 : /* -------------------------------------------------------------------- */
911 : /* Read the file. */
912 : /* -------------------------------------------------------------------- */
913 18 : if (!poDS->LoadFromFile(poDS->fp))
914 : {
915 0 : delete poDS;
916 0 : return nullptr;
917 : }
918 :
919 : /* -------------------------------------------------------------------- */
920 : /* Confirm the requested access is supported. */
921 : /* -------------------------------------------------------------------- */
922 18 : if (poOpenInfo->eAccess == GA_Update)
923 : {
924 0 : delete poDS;
925 0 : ReportUpdateNotSupportedByDriver("USGSDEM");
926 0 : return nullptr;
927 : }
928 :
929 : /* -------------------------------------------------------------------- */
930 : /* Create band information objects. */
931 : /* -------------------------------------------------------------------- */
932 18 : poDS->SetBand(1, new USGSDEMRasterBand(poDS));
933 :
934 18 : poDS->SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
935 :
936 : /* -------------------------------------------------------------------- */
937 : /* Initialize any PAM information. */
938 : /* -------------------------------------------------------------------- */
939 18 : poDS->SetDescription(poOpenInfo->pszFilename);
940 18 : poDS->TryLoadXML();
941 :
942 : /* -------------------------------------------------------------------- */
943 : /* Open overviews. */
944 : /* -------------------------------------------------------------------- */
945 18 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
946 :
947 18 : return poDS;
948 : }
949 :
950 : /************************************************************************/
951 : /* GDALRegister_USGSDEM() */
952 : /************************************************************************/
953 :
954 2033 : void GDALRegister_USGSDEM()
955 :
956 : {
957 2033 : if (GDALGetDriverByName("USGSDEM") != nullptr)
958 283 : return;
959 :
960 1750 : GDALDriver *poDriver = new GDALDriver();
961 :
962 1750 : poDriver->SetDescription("USGSDEM");
963 1750 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
964 1750 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "dem");
965 1750 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
966 1750 : "USGS Optional ASCII DEM (and CDED)");
967 1750 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
968 1750 : "drivers/raster/usgsdem.html");
969 :
970 1750 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
971 :
972 1750 : poDriver->pfnOpen = USGSDEMDataset::Open;
973 1750 : poDriver->pfnIdentify = USGSDEMDataset::Identify;
974 :
975 1750 : GetGDALDriverManager()->RegisterDriver(poDriver);
976 : }
|
