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 248 : static int ReadInt(VSILFILE *fp)
45 : {
46 : char c;
47 248 : int nRead = 0;
48 : char szBuffer[12];
49 248 : bool bInProlog = true;
50 :
51 : while (true)
52 : {
53 2480 : if (VSIFReadL(&c, 1, 1, fp) != 1)
54 : {
55 4 : return 0;
56 : }
57 2476 : if (bInProlog)
58 : {
59 2152 : if (!isspace(static_cast<unsigned char>(c)))
60 : {
61 244 : bInProlog = false;
62 : }
63 : }
64 2476 : if (!bInProlog)
65 : {
66 568 : if (c != '-' && c != '+' && !(c >= '0' && c <= '9'))
67 : {
68 244 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, VSIFTellL(fp) - 1, SEEK_SET));
69 244 : break;
70 : }
71 324 : if (nRead < 11)
72 324 : szBuffer[nRead] = c;
73 324 : nRead++;
74 : }
75 : }
76 244 : szBuffer[std::min(nRead, 11)] = 0;
77 244 : 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 75 : static void USGSDEMRefillBuffer(Buffer *psBuffer)
94 : {
95 75 : memmove(psBuffer->buffer, psBuffer->buffer + psBuffer->cur_index,
96 75 : psBuffer->buffer_size - psBuffer->cur_index);
97 :
98 75 : psBuffer->buffer_size -= psBuffer->cur_index;
99 75 : psBuffer->buffer_size += static_cast<int>(
100 150 : VSIFReadL(psBuffer->buffer + psBuffer->buffer_size, 1,
101 75 : psBuffer->max_size - psBuffer->buffer_size, psBuffer->fp));
102 75 : psBuffer->cur_index = 0;
103 75 : }
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 1863830 : static int USGSDEMReadIntFromBuffer(Buffer *psBuffer, int *pbSuccess = nullptr)
140 : {
141 : char c;
142 :
143 : while (true)
144 : {
145 1863830 : if (psBuffer->cur_index >= psBuffer->buffer_size)
146 : {
147 67 : USGSDEMRefillBuffer(psBuffer);
148 67 : if (psBuffer->cur_index >= psBuffer->buffer_size)
149 : {
150 0 : if (pbSuccess)
151 0 : *pbSuccess = FALSE;
152 0 : return 0;
153 : }
154 : }
155 :
156 1863830 : c = psBuffer->buffer[psBuffer->cur_index];
157 1863830 : psBuffer->cur_index++;
158 1863830 : if (!isspace(static_cast<unsigned char>(c)))
159 20595 : break;
160 : }
161 :
162 20595 : GIntBig nVal = 0;
163 20595 : int nSign = 1;
164 20595 : if (c == '-')
165 4944 : nSign = -1;
166 15651 : else if (c == '+')
167 0 : nSign = 1;
168 15651 : else if (c >= '0' && c <= '9')
169 15651 : 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 59465 : 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 59465 : c = psBuffer->buffer[psBuffer->cur_index];
191 59465 : if (c >= '0' && c <= '9')
192 : {
193 38870 : psBuffer->cur_index++;
194 38870 : if (nVal * nSign < INT_MAX && nVal * nSign > INT_MIN)
195 : {
196 38870 : nVal = nVal * 10 + (c - '0');
197 38870 : if (nVal * nSign > INT_MAX)
198 : {
199 0 : nVal = INT_MAX;
200 0 : nSign = 1;
201 : }
202 38870 : else if (nVal * nSign < INT_MIN)
203 : {
204 0 : nVal = INT_MIN;
205 0 : nSign = 1;
206 : }
207 : }
208 : }
209 : else
210 : {
211 20595 : if (pbSuccess)
212 20595 : *pbSuccess = TRUE;
213 20595 : return static_cast<int>(nSign * nVal);
214 : }
215 : }
216 : }
217 :
218 : /************************************************************************/
219 : /* USGSDEMReadDoubleFromBuffer() */
220 : /************************************************************************/
221 :
222 10304 : static double USGSDEMReadDoubleFromBuffer(Buffer *psBuffer, int nCharCount,
223 : int *pbSuccess = nullptr)
224 :
225 : {
226 10304 : if (psBuffer->cur_index + nCharCount > psBuffer->buffer_size)
227 : {
228 8 : USGSDEMRefillBuffer(psBuffer);
229 8 : 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 10303 : char *szPtr = psBuffer->buffer + psBuffer->cur_index;
238 10303 : char backupC = szPtr[nCharCount];
239 10303 : szPtr[nCharCount] = 0;
240 257575 : for (int i = 0; i < nCharCount; i++)
241 : {
242 247272 : if (szPtr[i] == 'D')
243 10284 : szPtr[i] = 'E';
244 : }
245 :
246 10303 : double dfVal = CPLAtof(szPtr);
247 10303 : szPtr[nCharCount] = backupC;
248 10303 : psBuffer->cur_index += nCharCount;
249 :
250 10303 : if (pbSuccess)
251 10303 : *pbSuccess = TRUE;
252 10303 : return dfVal;
253 : }
254 :
255 : /************************************************************************/
256 : /* DConvert() */
257 : /************************************************************************/
258 :
259 274 : static double DConvert(VSILFILE *fp, int nCharCount)
260 :
261 : {
262 : char szBuffer[100];
263 :
264 274 : CPL_IGNORE_RET_VAL(VSIFReadL(szBuffer, nCharCount, 1, fp));
265 274 : szBuffer[nCharCount] = '\0';
266 :
267 7090 : for (int i = 0; i < nCharCount; i++)
268 : {
269 6816 : if (szBuffer[i] == 'D')
270 208 : szBuffer[i] = 'E';
271 : }
272 :
273 548 : 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 : vsi_l_offset 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 20 : USGSDEMRasterBand::USGSDEMRasterBand(USGSDEMDataset *poDSIn)
335 :
336 : {
337 20 : this->poDS = poDSIn;
338 20 : this->nBand = 1;
339 :
340 20 : eDataType = poDSIn->eNaturalDataFormat;
341 :
342 20 : nBlockXSize = poDSIn->GetRasterXSize();
343 20 : nBlockYSize = poDSIn->GetRasterYSize();
344 20 : }
345 :
346 : /************************************************************************/
347 : /* IReadBlock() */
348 : /************************************************************************/
349 :
350 10 : CPLErr USGSDEMRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
351 : CPL_UNUSED int nBlockYOff, void *pImage)
352 :
353 : {
354 : /* int bad = FALSE; */
355 10 : USGSDEMDataset *poGDS = cpl::down_cast<USGSDEMDataset *>(poDS);
356 :
357 : /* -------------------------------------------------------------------- */
358 : /* Initialize image buffer to nodata value. */
359 : /* -------------------------------------------------------------------- */
360 10 : GDALCopyWords(&USGSDEM_NODATA, GDT_Int32, 0, pImage, GetRasterDataType(),
361 : GDALGetDataTypeSizeBytes(GetRasterDataType()),
362 10 : GetXSize() * GetYSize());
363 :
364 : /* -------------------------------------------------------------------- */
365 : /* Seek to data. */
366 : /* -------------------------------------------------------------------- */
367 10 : CPL_IGNORE_RET_VAL(VSIFSeekL(poGDS->fp, poGDS->nDataStartOffset, 0));
368 :
369 10 : double dfYMin = poGDS->m_gt.yorig + (GetYSize() - 0.5) * poGDS->m_gt.yscale;
370 :
371 : /* -------------------------------------------------------------------- */
372 : /* Read all the profiles into the image buffer. */
373 : /* -------------------------------------------------------------------- */
374 :
375 : Buffer sBuffer;
376 10 : sBuffer.max_size = 32768;
377 10 : sBuffer.buffer = static_cast<char *>(CPLMalloc(sBuffer.max_size + 1));
378 10 : sBuffer.fp = poGDS->fp;
379 10 : sBuffer.buffer_size = 0;
380 10 : sBuffer.cur_index = 0;
381 :
382 2070 : for (int i = 0; i < GetXSize(); i++)
383 : {
384 : int bSuccess;
385 2061 : const int nRowNumber = USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
386 2061 : if (nRowNumber != 1)
387 1 : CPLDebug("USGSDEM", "i = %d, nRowNumber = %d", i, nRowNumber);
388 2061 : if (bSuccess)
389 : {
390 : const int nColNumber =
391 2061 : USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
392 2061 : if (nColNumber != i + 1)
393 : {
394 3 : CPLDebug("USGSDEM", "i = %d, nColNumber = %d", i, nColNumber);
395 : }
396 : }
397 : const int nCPoints =
398 2061 : (bSuccess) ? USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess) : 0;
399 : #ifdef DEBUG_VERBOSE
400 : CPLDebug("USGSDEM", "i = %d, nCPoints = %d", i, nCPoints);
401 : #endif
402 :
403 2061 : if (bSuccess)
404 : {
405 : const int nNumberOfCols =
406 2061 : USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
407 2061 : if (nNumberOfCols != 1)
408 : {
409 0 : CPLDebug("USGSDEM", "i = %d, nNumberOfCols = %d", i,
410 : nNumberOfCols);
411 : }
412 : }
413 :
414 : // x-start
415 2061 : if (bSuccess)
416 2061 : /* dxStart = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24,
417 : &bSuccess);
418 :
419 : double dyStart =
420 2061 : (bSuccess) ? USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess)
421 2061 : : 0;
422 : const double dfElevOffset =
423 2061 : (bSuccess) ? USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess)
424 2061 : : 0;
425 :
426 : // min z value
427 2061 : if (bSuccess)
428 2061 : /* djunk = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess);
429 :
430 : // max z value
431 2061 : if (bSuccess)
432 2060 : /* djunk = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess);
433 2061 : if (!bSuccess)
434 : {
435 1 : CPLFree(sBuffer.buffer);
436 1 : return CE_Failure;
437 : }
438 :
439 2060 : if (poGDS->m_oSRS.IsGeographic())
440 4 : dyStart = dyStart / 3600.0;
441 :
442 2060 : double dygap = (dfYMin - dyStart) / poGDS->m_gt.yscale + 0.5;
443 2060 : if (dygap <= INT_MIN || dygap >= INT_MAX || !std::isfinite(dygap))
444 : {
445 0 : CPLFree(sBuffer.buffer);
446 0 : return CE_Failure;
447 : }
448 2060 : int lygap = static_cast<int>(dygap);
449 2060 : if (nCPoints <= 0)
450 0 : continue;
451 2060 : if (lygap > INT_MAX - nCPoints)
452 0 : lygap = INT_MAX - nCPoints;
453 2060 : if (lygap < 0 && GetYSize() > INT_MAX + lygap)
454 : {
455 0 : CPLFree(sBuffer.buffer);
456 0 : return CE_Failure;
457 : }
458 :
459 14411 : for (int j = lygap; j < (nCPoints + lygap); j++)
460 : {
461 12351 : const int iY = GetYSize() - j - 1;
462 :
463 12351 : 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 12351 : if (!bSuccess)
470 : {
471 0 : CPLFree(sBuffer.buffer);
472 0 : return CE_Failure;
473 : }
474 :
475 12351 : if (iY < 0 || iY >= GetYSize())
476 : {
477 : /* bad = TRUE; */
478 : }
479 12351 : else if (nElev == USGSDEM_NODATA)
480 : /* leave in output buffer as nodata */;
481 : else
482 : {
483 8389 : const float fComputedElev =
484 8389 : static_cast<float>(nElev * poGDS->fVRes + dfElevOffset);
485 :
486 8389 : if (GetRasterDataType() == GDT_Int16)
487 : {
488 16656 : GUInt16 nVal = (fComputedElev < -32768) ? -32768
489 : : (fComputedElev > 32767)
490 : ? 32767
491 8328 : : static_cast<GInt16>(fComputedElev);
492 8328 : reinterpret_cast<GInt16 *>(pImage)[i + iY * GetXSize()] =
493 8328 : nVal;
494 : }
495 : else
496 : {
497 61 : reinterpret_cast<float *>(pImage)[i + iY * GetXSize()] =
498 : fComputedElev;
499 : }
500 : }
501 : }
502 :
503 2060 : 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 9 : CPLFree(sBuffer.buffer);
517 :
518 9 : 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 20 : USGSDEMDataset::USGSDEMDataset()
555 : : nDataStartOffset(0), eNaturalDataFormat(GDT_Unknown), fVRes(0.0),
556 20 : pszUnits(nullptr), fp(nullptr)
557 : {
558 20 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
559 20 : }
560 :
561 : /************************************************************************/
562 : /* ~USGSDEMDataset() */
563 : /************************************************************************/
564 :
565 40 : USGSDEMDataset::~USGSDEMDataset()
566 :
567 : {
568 20 : FlushCache(true);
569 :
570 20 : if (fp != nullptr)
571 20 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
572 40 : }
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 20 : int USGSDEMDataset::LoadFromFile(VSILFILE *InDem)
585 : {
586 : // check for version of DEM format
587 20 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 864, 0));
588 :
589 : // Read DEM into matrix
590 20 : const int nRow = ReadInt(InDem);
591 20 : const int nColumn = ReadInt(InDem);
592 : const bool bNewFormat =
593 20 : VSIFTellL(InDem) >= 1024 || nRow != 1 || nColumn != 1;
594 20 : if (bNewFormat)
595 : {
596 20 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 1024, 0)); // New Format
597 20 : int i = ReadInt(InDem);
598 20 : int j = ReadInt(InDem);
599 20 : 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 2 : 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 16 : nDataStartOffset = 1024;
627 : }
628 : }
629 : else
630 0 : nDataStartOffset = 864;
631 :
632 20 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, nDataStartOffset, 0));
633 : char c;
634 20 : if (VSIFReadL(&c, 1, 1, InDem) == 1)
635 : {
636 2 : if (c == '\n' && VSIFSeekL(InDem, 2 * nDataStartOffset + 1, 0) == 0 &&
637 22 : VSIFReadL(&c, 1, 1, InDem) == 1 && c == '\n')
638 : {
639 : // Some files use 1025 byte records ending with a newline
640 : // character.
641 : // See https://github.com/OSGeo/gdal/issues/5007
642 2 : nDataStartOffset++;
643 : }
644 2 : else if (c == '\r' && VSIFReadL(&c, 1, 1, InDem) == 1 && c == '\n' &&
645 2 : VSIFSeekL(InDem, 2 * nDataStartOffset + 2, 0) == 0 &&
646 2 : VSIFReadL(&c, 1, 1, InDem) == 1 && c == '\r' &&
647 20 : VSIFReadL(&c, 1, 1, InDem) == 1 && c == '\n')
648 : {
649 2 : nDataStartOffset += 2;
650 : }
651 : }
652 :
653 20 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 156, 0));
654 20 : const int nCoordSystem = ReadInt(InDem);
655 20 : const int iUTMZone = ReadInt(InDem);
656 :
657 20 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 528, 0));
658 20 : const int nGUnit = ReadInt(InDem);
659 20 : const int nVUnit = ReadInt(InDem);
660 :
661 : // Vertical Units in meters
662 20 : if (nVUnit == 1)
663 0 : pszUnits = "ft";
664 : else
665 20 : pszUnits = "m";
666 :
667 20 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 816, 0));
668 20 : const double dxdelta = DConvert(InDem, 12);
669 20 : const double dydelta = DConvert(InDem, 12);
670 20 : if (dydelta == 0)
671 0 : return FALSE;
672 20 : fVRes = DConvert(InDem, 12);
673 :
674 : /* -------------------------------------------------------------------- */
675 : /* Should we treat this as floating point, or GInt16. */
676 : /* -------------------------------------------------------------------- */
677 20 : if (nVUnit == 1 || fVRes < 1.0)
678 4 : eNaturalDataFormat = GDT_Float32;
679 : else
680 16 : eNaturalDataFormat = GDT_Int16;
681 :
682 : /* -------------------------------------------------------------------- */
683 : /* Read four corner coordinates. */
684 : /* -------------------------------------------------------------------- */
685 20 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 546, 0));
686 : DPoint2 corners[4]; // SW, NW, NE, SE
687 100 : for (int i = 0; i < 4; i++)
688 : {
689 80 : corners[i].x = DConvert(InDem, 24);
690 80 : corners[i].y = DConvert(InDem, 24);
691 : }
692 :
693 : // find absolute extents of raw vales
694 : DPoint2 extent_min, extent_max;
695 20 : extent_min.x = std::min(corners[0].x, corners[1].x);
696 20 : extent_max.x = std::max(corners[2].x, corners[3].x);
697 20 : extent_min.y = std::min(corners[0].y, corners[3].y);
698 20 : extent_max.y = std::max(corners[1].y, corners[2].y);
699 :
700 20 : /* dElevMin = */ DConvert(InDem, 48);
701 20 : /* dElevMax = */ DConvert(InDem, 48);
702 :
703 20 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 858, 0));
704 20 : const int nProfiles = ReadInt(InDem);
705 :
706 : /* -------------------------------------------------------------------- */
707 : /* Collect the spatial reference system. */
708 : /* -------------------------------------------------------------------- */
709 40 : OGRSpatialReference sr;
710 20 : sr.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
711 20 : bool bNAD83 = true;
712 :
713 : // OLD format header ends at byte 864
714 20 : if (bNewFormat)
715 : {
716 : // year of data compilation
717 20 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 876, 0));
718 : char szDateBuffer[5];
719 20 : CPL_IGNORE_RET_VAL(VSIFReadL(szDateBuffer, 4, 1, InDem));
720 : /* szDateBuffer[4] = 0; */
721 :
722 : // Horizontal datum
723 : // 1=North American Datum 1927 (NAD 27)
724 : // 2=World Geodetic System 1972 (WGS 72)
725 : // 3=WGS 84
726 : // 4=NAD 83
727 : // 5=Old Hawaii Datum
728 : // 6=Puerto Rico Datum
729 20 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 890, 0));
730 :
731 : char szHorzDatum[3];
732 20 : CPL_IGNORE_RET_VAL(VSIFReadL(szHorzDatum, 1, 2, InDem));
733 20 : szHorzDatum[2] = '\0';
734 20 : const int datum = atoi(szHorzDatum);
735 20 : switch (datum)
736 : {
737 6 : case 1:
738 6 : sr.SetWellKnownGeogCS("NAD27");
739 6 : bNAD83 = false;
740 6 : break;
741 :
742 2 : case 2:
743 2 : sr.SetWellKnownGeogCS("WGS72");
744 2 : break;
745 :
746 0 : case 3:
747 0 : sr.SetWellKnownGeogCS("WGS84");
748 0 : break;
749 :
750 2 : case 4:
751 2 : sr.SetWellKnownGeogCS("NAD83");
752 2 : break;
753 :
754 0 : case -9:
755 0 : break;
756 :
757 10 : default:
758 10 : sr.SetWellKnownGeogCS("NAD27");
759 10 : break;
760 : }
761 : }
762 : else
763 : {
764 0 : sr.SetWellKnownGeogCS("NAD27");
765 0 : bNAD83 = false;
766 : }
767 :
768 20 : if (nCoordSystem == 1) // UTM
769 : {
770 14 : if (iUTMZone >= -60 && iUTMZone <= 60)
771 : {
772 14 : sr.SetUTM(abs(iUTMZone), iUTMZone >= 0);
773 14 : if (nGUnit == 1)
774 : {
775 0 : sr.SetLinearUnitsAndUpdateParameters(
776 : SRS_UL_US_FOOT, CPLAtof(SRS_UL_US_FOOT_CONV));
777 : char szUTMName[128];
778 0 : snprintf(szUTMName, sizeof(szUTMName),
779 : "UTM Zone %d, Northern Hemisphere, us-ft", iUTMZone);
780 0 : sr.SetNode("PROJCS", szUTMName);
781 : }
782 : }
783 : }
784 6 : else if (nCoordSystem == 2) // state plane
785 : {
786 0 : if (nGUnit == 1)
787 0 : sr.SetStatePlane(iUTMZone, bNAD83, "Foot",
788 : CPLAtof(SRS_UL_US_FOOT_CONV));
789 : else
790 0 : sr.SetStatePlane(iUTMZone, bNAD83);
791 : }
792 :
793 20 : m_oSRS = std::move(sr);
794 :
795 : /* -------------------------------------------------------------------- */
796 : /* For UTM we use the extents (really the UTM coordinates of */
797 : /* the lat/long corners of the quad) to determine the size in */
798 : /* pixels and lines, but we have to make the anchors be modulus */
799 : /* the pixel size which what really gets used. */
800 : /* -------------------------------------------------------------------- */
801 20 : if (nCoordSystem == 1 // UTM
802 6 : || nCoordSystem == 2 // State Plane
803 6 : || nCoordSystem == -9999) // unknown
804 : {
805 : // expand extents modulus the pixel size.
806 14 : extent_min.y = floor(extent_min.y / dydelta) * dydelta;
807 14 : extent_max.y = ceil(extent_max.y / dydelta) * dydelta;
808 :
809 : // Forcibly compute X extents based on first profile and pixelsize.
810 14 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, nDataStartOffset, 0));
811 14 : /* njunk = */ ReadInt(InDem);
812 14 : /* njunk = */ ReadInt(InDem);
813 14 : /* njunk = */ ReadInt(InDem);
814 14 : /* njunk = */ ReadInt(InDem);
815 14 : const double dxStart = DConvert(InDem, 24);
816 :
817 14 : const double dfRasterYSize =
818 14 : (extent_max.y - extent_min.y) / dydelta + 1.5;
819 28 : if (dfRasterYSize <= INT_MIN || dfRasterYSize >= INT_MAX ||
820 14 : !std::isfinite(dfRasterYSize))
821 0 : return FALSE;
822 14 : nRasterYSize = static_cast<int>(dfRasterYSize);
823 14 : nRasterXSize = nProfiles;
824 :
825 14 : m_gt.xorig = dxStart - dxdelta / 2.0;
826 14 : m_gt.xscale = dxdelta;
827 14 : m_gt.xrot = 0.0;
828 14 : m_gt.yorig = extent_max.y + dydelta / 2.0;
829 14 : m_gt.yrot = 0.0;
830 14 : m_gt.yscale = -dydelta;
831 : }
832 : /* -------------------------------------------------------------------- */
833 : /* Geographic -- use corners directly. */
834 : /* -------------------------------------------------------------------- */
835 : else
836 : {
837 6 : const double dfRasterYSize =
838 6 : (extent_max.y - extent_min.y) / dydelta + 1.5;
839 12 : if (dfRasterYSize <= INT_MIN || dfRasterYSize >= INT_MAX ||
840 6 : !std::isfinite(dfRasterYSize))
841 0 : return FALSE;
842 6 : nRasterYSize = static_cast<int>(dfRasterYSize);
843 6 : nRasterXSize = nProfiles;
844 :
845 : // Translate extents from arc-seconds to decimal degrees.
846 6 : m_gt.xorig = (extent_min.x - dxdelta / 2.0) / 3600.0;
847 6 : m_gt.xscale = dxdelta / 3600.0;
848 6 : m_gt.xrot = 0.0;
849 6 : m_gt.yorig = (extent_max.y + dydelta / 2.0) / 3600.0;
850 6 : m_gt.yrot = 0.0;
851 6 : m_gt.yscale = (-dydelta) / 3600.0;
852 : }
853 :
854 : // IReadBlock() not ready for more than INT_MAX pixels, and that
855 : // would behave badly
856 40 : if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize) ||
857 20 : nRasterXSize > INT_MAX / nRasterYSize)
858 : {
859 0 : return FALSE;
860 : }
861 :
862 20 : return TRUE;
863 : }
864 :
865 : /************************************************************************/
866 : /* GetGeoTransform() */
867 : /************************************************************************/
868 :
869 6 : CPLErr USGSDEMDataset::GetGeoTransform(GDALGeoTransform >) const
870 :
871 : {
872 6 : gt = m_gt;
873 6 : return CE_None;
874 : }
875 :
876 : /************************************************************************/
877 : /* GetSpatialRef() */
878 : /************************************************************************/
879 :
880 6 : const OGRSpatialReference *USGSDEMDataset::GetSpatialRef() const
881 :
882 : {
883 6 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
884 : }
885 :
886 : /************************************************************************/
887 : /* Identify() */
888 : /************************************************************************/
889 :
890 63131 : int USGSDEMDataset::Identify(GDALOpenInfo *poOpenInfo)
891 :
892 : {
893 63131 : if (poOpenInfo->nHeaderBytes < 200)
894 58229 : return FALSE;
895 :
896 4902 : if (!STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 0") &&
897 4890 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 1") &&
898 4862 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 2") &&
899 4862 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 3") &&
900 4862 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " -9999"))
901 4862 : return FALSE;
902 :
903 40 : if (!STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 150, " 1") &&
904 4 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 150, " 4"))
905 0 : return FALSE;
906 :
907 40 : return TRUE;
908 : }
909 :
910 : /************************************************************************/
911 : /* Open() */
912 : /************************************************************************/
913 :
914 20 : GDALDataset *USGSDEMDataset::Open(GDALOpenInfo *poOpenInfo)
915 :
916 : {
917 20 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
918 0 : return nullptr;
919 :
920 : /* -------------------------------------------------------------------- */
921 : /* Create a corresponding GDALDataset. */
922 : /* -------------------------------------------------------------------- */
923 20 : USGSDEMDataset *poDS = new USGSDEMDataset();
924 :
925 20 : poDS->fp = poOpenInfo->fpL;
926 20 : poOpenInfo->fpL = nullptr;
927 :
928 : /* -------------------------------------------------------------------- */
929 : /* Read the file. */
930 : /* -------------------------------------------------------------------- */
931 20 : if (!poDS->LoadFromFile(poDS->fp))
932 : {
933 0 : delete poDS;
934 0 : return nullptr;
935 : }
936 :
937 : /* -------------------------------------------------------------------- */
938 : /* Confirm the requested access is supported. */
939 : /* -------------------------------------------------------------------- */
940 20 : if (poOpenInfo->eAccess == GA_Update)
941 : {
942 0 : delete poDS;
943 0 : ReportUpdateNotSupportedByDriver("USGSDEM");
944 0 : return nullptr;
945 : }
946 :
947 : /* -------------------------------------------------------------------- */
948 : /* Create band information objects. */
949 : /* -------------------------------------------------------------------- */
950 20 : poDS->SetBand(1, new USGSDEMRasterBand(poDS));
951 :
952 20 : poDS->SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
953 :
954 : /* -------------------------------------------------------------------- */
955 : /* Initialize any PAM information. */
956 : /* -------------------------------------------------------------------- */
957 20 : poDS->SetDescription(poOpenInfo->pszFilename);
958 20 : poDS->TryLoadXML();
959 :
960 : /* -------------------------------------------------------------------- */
961 : /* Open overviews. */
962 : /* -------------------------------------------------------------------- */
963 20 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
964 :
965 20 : return poDS;
966 : }
967 :
968 : /************************************************************************/
969 : /* GDALRegister_USGSDEM() */
970 : /************************************************************************/
971 :
972 2062 : void GDALRegister_USGSDEM()
973 :
974 : {
975 2062 : if (GDALGetDriverByName("USGSDEM") != nullptr)
976 263 : return;
977 :
978 1799 : GDALDriver *poDriver = new GDALDriver();
979 :
980 1799 : poDriver->SetDescription("USGSDEM");
981 1799 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
982 1799 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "dem");
983 1799 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
984 1799 : "USGS Optional ASCII DEM (and CDED)");
985 1799 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
986 1799 : "drivers/raster/usgsdem.html");
987 :
988 1799 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
989 :
990 1799 : poDriver->pfnOpen = USGSDEMDataset::Open;
991 1799 : poDriver->pfnIdentify = USGSDEMDataset::Identify;
992 :
993 1799 : GetGDALDriverManager()->RegisterDriver(poDriver);
994 : }
|
