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