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 : double adfGeoTransform[6];
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(double *padfTransform) 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 = reinterpret_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->adfGeoTransform[3] +
366 9 : (GetYSize() - 0.5) * poGDS->adfGeoTransform[5];
367 :
368 : /* -------------------------------------------------------------------- */
369 : /* Read all the profiles into the image buffer. */
370 : /* -------------------------------------------------------------------- */
371 :
372 : Buffer sBuffer;
373 9 : sBuffer.max_size = 32768;
374 9 : sBuffer.buffer = reinterpret_cast<char *>(CPLMalloc(sBuffer.max_size + 1));
375 9 : sBuffer.fp = poGDS->fp;
376 9 : sBuffer.buffer_size = 0;
377 9 : sBuffer.cur_index = 0;
378 :
379 1045 : for (int i = 0; i < GetXSize(); i++)
380 : {
381 : int bSuccess;
382 1037 : const int nRowNumber = USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
383 1037 : if (nRowNumber != 1)
384 1 : CPLDebug("USGSDEM", "i = %d, nRowNumber = %d", i, nRowNumber);
385 1037 : if (bSuccess)
386 : {
387 : const int nColNumber =
388 1037 : USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
389 1037 : if (nColNumber != i + 1)
390 : {
391 3 : CPLDebug("USGSDEM", "i = %d, nColNumber = %d", i, nColNumber);
392 : }
393 : }
394 : const int nCPoints =
395 1037 : (bSuccess) ? USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess) : 0;
396 : #ifdef DEBUG_VERBOSE
397 : CPLDebug("USGSDEM", "i = %d, nCPoints = %d", i, nCPoints);
398 : #endif
399 :
400 1037 : if (bSuccess)
401 : {
402 : const int nNumberOfCols =
403 1037 : USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
404 1037 : if (nNumberOfCols != 1)
405 : {
406 0 : CPLDebug("USGSDEM", "i = %d, nNumberOfCols = %d", i,
407 : nNumberOfCols);
408 : }
409 : }
410 :
411 : // x-start
412 1037 : if (bSuccess)
413 1037 : /* dxStart = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24,
414 : &bSuccess);
415 :
416 : double dyStart =
417 1037 : (bSuccess) ? USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess)
418 1037 : : 0;
419 : const double dfElevOffset =
420 1037 : (bSuccess) ? USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess)
421 1037 : : 0;
422 :
423 : // min z value
424 1037 : if (bSuccess)
425 1037 : /* djunk = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess);
426 :
427 : // max z value
428 1037 : if (bSuccess)
429 1036 : /* djunk = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess);
430 1037 : if (!bSuccess)
431 : {
432 1 : CPLFree(sBuffer.buffer);
433 1 : return CE_Failure;
434 : }
435 :
436 1036 : if (poGDS->m_oSRS.IsGeographic())
437 4 : dyStart = dyStart / 3600.0;
438 :
439 1036 : double dygap = (dfYMin - dyStart) / poGDS->adfGeoTransform[5] + 0.5;
440 1036 : if (dygap <= INT_MIN || dygap >= INT_MAX || !std::isfinite(dygap))
441 : {
442 0 : CPLFree(sBuffer.buffer);
443 0 : return CE_Failure;
444 : }
445 1036 : int lygap = static_cast<int>(dygap);
446 1036 : if (nCPoints <= 0)
447 0 : continue;
448 1036 : if (lygap > INT_MAX - nCPoints)
449 0 : lygap = INT_MAX - nCPoints;
450 1036 : if (lygap < 0 && GetYSize() > INT_MAX + lygap)
451 : {
452 0 : CPLFree(sBuffer.buffer);
453 0 : return CE_Failure;
454 : }
455 :
456 11339 : for (int j = lygap; j < (nCPoints + lygap); j++)
457 : {
458 10303 : const int iY = GetYSize() - j - 1;
459 :
460 10303 : const int nElev = USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
461 : #ifdef DEBUG_VERBOSE
462 : CPLDebug("USGSDEM", " j - lygap = %d, nElev = %d", j - lygap,
463 : nElev);
464 : #endif
465 :
466 10303 : if (!bSuccess)
467 : {
468 0 : CPLFree(sBuffer.buffer);
469 0 : return CE_Failure;
470 : }
471 :
472 10303 : if (iY < 0 || iY >= GetYSize())
473 : {
474 : /* bad = TRUE; */
475 : }
476 10303 : else if (nElev == USGSDEM_NODATA)
477 : /* leave in output buffer as nodata */;
478 : else
479 : {
480 6341 : const float fComputedElev =
481 6341 : static_cast<float>(nElev * poGDS->fVRes + dfElevOffset);
482 :
483 6341 : if (GetRasterDataType() == GDT_Int16)
484 : {
485 12560 : GUInt16 nVal = (fComputedElev < -32768) ? -32768
486 : : (fComputedElev > 32767)
487 : ? 32767
488 6280 : : static_cast<GInt16>(fComputedElev);
489 6280 : reinterpret_cast<GInt16 *>(pImage)[i + iY * GetXSize()] =
490 6280 : nVal;
491 : }
492 : else
493 : {
494 61 : reinterpret_cast<float *>(pImage)[i + iY * GetXSize()] =
495 : fComputedElev;
496 : }
497 : }
498 : }
499 :
500 1036 : if (poGDS->nDataStartOffset == 1024)
501 : {
502 : // Seek to the next 1024 byte boundary.
503 : // Some files have 'junk' profile values after the valid/declared
504 : // ones
505 10 : vsi_l_offset nCurPos = USGSDEMGetCurrentFilePos(&sBuffer);
506 10 : vsi_l_offset nNewPos = (nCurPos + 1023) / 1024 * 1024;
507 10 : if (nNewPos > nCurPos)
508 : {
509 10 : USGSDEMSetCurrentFilePos(&sBuffer, nNewPos);
510 : }
511 : }
512 : }
513 8 : CPLFree(sBuffer.buffer);
514 :
515 8 : return CE_None;
516 : }
517 :
518 : /************************************************************************/
519 : /* GetNoDataValue() */
520 : /************************************************************************/
521 :
522 0 : double USGSDEMRasterBand::GetNoDataValue(int *pbSuccess)
523 :
524 : {
525 0 : if (pbSuccess != nullptr)
526 0 : *pbSuccess = TRUE;
527 :
528 0 : return USGSDEM_NODATA;
529 : }
530 :
531 : /************************************************************************/
532 : /* GetUnitType() */
533 : /************************************************************************/
534 0 : const char *USGSDEMRasterBand::GetUnitType()
535 : {
536 0 : USGSDEMDataset *poGDS = reinterpret_cast<USGSDEMDataset *>(poDS);
537 :
538 0 : return poGDS->pszUnits;
539 : }
540 :
541 : /************************************************************************/
542 : /* ==================================================================== */
543 : /* USGSDEMDataset */
544 : /* ==================================================================== */
545 : /************************************************************************/
546 :
547 : /************************************************************************/
548 : /* USGSDEMDataset() */
549 : /************************************************************************/
550 :
551 18 : USGSDEMDataset::USGSDEMDataset()
552 : : nDataStartOffset(0), eNaturalDataFormat(GDT_Unknown), fVRes(0.0),
553 18 : pszUnits(nullptr), fp(nullptr)
554 : {
555 18 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
556 18 : memset(adfGeoTransform, 0, sizeof(adfGeoTransform));
557 18 : }
558 :
559 : /************************************************************************/
560 : /* ~USGSDEMDataset() */
561 : /************************************************************************/
562 :
563 36 : USGSDEMDataset::~USGSDEMDataset()
564 :
565 : {
566 18 : FlushCache(true);
567 :
568 18 : if (fp != nullptr)
569 18 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
570 36 : }
571 :
572 : /************************************************************************/
573 : /* LoadFromFile() */
574 : /* */
575 : /* If the data from DEM is in meters, then values are stored as */
576 : /* shorts. If DEM data is in feet, then height data will be */
577 : /* stored in float, to preserve the precision of the original */
578 : /* data. returns true if the file was successfully opened and */
579 : /* read. */
580 : /************************************************************************/
581 :
582 18 : int USGSDEMDataset::LoadFromFile(VSILFILE *InDem)
583 : {
584 : // check for version of DEM format
585 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 864, 0));
586 :
587 : // Read DEM into matrix
588 18 : const int nRow = ReadInt(InDem);
589 18 : const int nColumn = ReadInt(InDem);
590 : const bool bNewFormat =
591 18 : VSIFTellL(InDem) >= 1024 || nRow != 1 || nColumn != 1;
592 18 : if (bNewFormat)
593 : {
594 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 1024, 0)); // New Format
595 18 : int i = ReadInt(InDem);
596 18 : int j = ReadInt(InDem);
597 18 : if (i != 1 || (j != 1 && j != 0)) // File OK?
598 : {
599 4 : CPL_IGNORE_RET_VAL(
600 4 : VSIFSeekL(InDem, 893, 0)); // Undocumented Format (39109h1.dem)
601 4 : i = ReadInt(InDem);
602 4 : j = ReadInt(InDem);
603 4 : if (i != 1 || j != 1) // File OK?
604 : {
605 2 : CPL_IGNORE_RET_VAL(VSIFSeekL(
606 : InDem, 918, 0)); // Latest iteration of the A record, such
607 : // as in fema06-140cm_2995441b.dem
608 2 : i = ReadInt(InDem);
609 2 : j = ReadInt(InDem);
610 2 : if (i != 1 || j != 1) // File OK?
611 : {
612 0 : CPLError(CE_Failure, CPLE_AppDefined,
613 : "Does not appear to be a USGS DEM file.");
614 0 : return FALSE;
615 : }
616 : else
617 2 : nDataStartOffset = 918;
618 : }
619 : else
620 2 : nDataStartOffset = 893;
621 : }
622 : else
623 : {
624 14 : nDataStartOffset = 1024;
625 :
626 : // Some files use 1025 byte records ending with a newline character.
627 : // See https://github.com/OSGeo/gdal/issues/5007
628 14 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 1024, 0));
629 : char c;
630 28 : if (VSIFReadL(&c, 1, 1, InDem) == 1 && c == '\n' &&
631 2 : VSIFSeekL(InDem, 1024 + 1024 + 1, 0) == 0 &&
632 28 : VSIFReadL(&c, 1, 1, InDem) == 1 && c == '\n')
633 : {
634 2 : nDataStartOffset = 1025;
635 : }
636 : }
637 : }
638 : else
639 0 : nDataStartOffset = 864;
640 :
641 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 156, 0));
642 18 : const int nCoordSystem = ReadInt(InDem);
643 18 : const int iUTMZone = ReadInt(InDem);
644 :
645 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 528, 0));
646 18 : const int nGUnit = ReadInt(InDem);
647 18 : const int nVUnit = ReadInt(InDem);
648 :
649 : // Vertical Units in meters
650 18 : if (nVUnit == 1)
651 0 : pszUnits = "ft";
652 : else
653 18 : pszUnits = "m";
654 :
655 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 816, 0));
656 18 : const double dxdelta = DConvert(InDem, 12);
657 18 : const double dydelta = DConvert(InDem, 12);
658 18 : if (dydelta == 0)
659 0 : return FALSE;
660 18 : fVRes = DConvert(InDem, 12);
661 :
662 : /* -------------------------------------------------------------------- */
663 : /* Should we treat this as floating point, or GInt16. */
664 : /* -------------------------------------------------------------------- */
665 18 : if (nVUnit == 1 || fVRes < 1.0)
666 4 : eNaturalDataFormat = GDT_Float32;
667 : else
668 14 : eNaturalDataFormat = GDT_Int16;
669 :
670 : /* -------------------------------------------------------------------- */
671 : /* Read four corner coordinates. */
672 : /* -------------------------------------------------------------------- */
673 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 546, 0));
674 : DPoint2 corners[4]; // SW, NW, NE, SE
675 90 : for (int i = 0; i < 4; i++)
676 : {
677 72 : corners[i].x = DConvert(InDem, 24);
678 72 : corners[i].y = DConvert(InDem, 24);
679 : }
680 :
681 : // find absolute extents of raw vales
682 : DPoint2 extent_min, extent_max;
683 18 : extent_min.x = std::min(corners[0].x, corners[1].x);
684 18 : extent_max.x = std::max(corners[2].x, corners[3].x);
685 18 : extent_min.y = std::min(corners[0].y, corners[3].y);
686 18 : extent_max.y = std::max(corners[1].y, corners[2].y);
687 :
688 18 : /* dElevMin = */ DConvert(InDem, 48);
689 18 : /* dElevMax = */ DConvert(InDem, 48);
690 :
691 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 858, 0));
692 18 : const int nProfiles = ReadInt(InDem);
693 :
694 : /* -------------------------------------------------------------------- */
695 : /* Collect the spatial reference system. */
696 : /* -------------------------------------------------------------------- */
697 36 : OGRSpatialReference sr;
698 18 : sr.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
699 18 : bool bNAD83 = true;
700 :
701 : // OLD format header ends at byte 864
702 18 : if (bNewFormat)
703 : {
704 : // year of data compilation
705 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 876, 0));
706 : char szDateBuffer[5];
707 18 : CPL_IGNORE_RET_VAL(VSIFReadL(szDateBuffer, 4, 1, InDem));
708 : /* szDateBuffer[4] = 0; */
709 :
710 : // Horizontal datum
711 : // 1=North American Datum 1927 (NAD 27)
712 : // 2=World Geodetic System 1972 (WGS 72)
713 : // 3=WGS 84
714 : // 4=NAD 83
715 : // 5=Old Hawaii Datum
716 : // 6=Puerto Rico Datum
717 18 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 890, 0));
718 :
719 : char szHorzDatum[3];
720 18 : CPL_IGNORE_RET_VAL(VSIFReadL(szHorzDatum, 1, 2, InDem));
721 18 : szHorzDatum[2] = '\0';
722 18 : const int datum = atoi(szHorzDatum);
723 18 : switch (datum)
724 : {
725 4 : case 1:
726 4 : sr.SetWellKnownGeogCS("NAD27");
727 4 : bNAD83 = false;
728 4 : break;
729 :
730 2 : case 2:
731 2 : sr.SetWellKnownGeogCS("WGS72");
732 2 : break;
733 :
734 0 : case 3:
735 0 : sr.SetWellKnownGeogCS("WGS84");
736 0 : break;
737 :
738 2 : case 4:
739 2 : sr.SetWellKnownGeogCS("NAD83");
740 2 : break;
741 :
742 0 : case -9:
743 0 : break;
744 :
745 10 : default:
746 10 : sr.SetWellKnownGeogCS("NAD27");
747 10 : break;
748 : }
749 : }
750 : else
751 : {
752 0 : sr.SetWellKnownGeogCS("NAD27");
753 0 : bNAD83 = false;
754 : }
755 :
756 18 : if (nCoordSystem == 1) // UTM
757 : {
758 12 : if (iUTMZone >= -60 && iUTMZone <= 60)
759 : {
760 12 : sr.SetUTM(abs(iUTMZone), iUTMZone >= 0);
761 12 : if (nGUnit == 1)
762 : {
763 0 : sr.SetLinearUnitsAndUpdateParameters(
764 : SRS_UL_US_FOOT, CPLAtof(SRS_UL_US_FOOT_CONV));
765 : char szUTMName[128];
766 0 : snprintf(szUTMName, sizeof(szUTMName),
767 : "UTM Zone %d, Northern Hemisphere, us-ft", iUTMZone);
768 0 : sr.SetNode("PROJCS", szUTMName);
769 : }
770 : }
771 : }
772 6 : else if (nCoordSystem == 2) // state plane
773 : {
774 0 : if (nGUnit == 1)
775 0 : sr.SetStatePlane(iUTMZone, bNAD83, "Foot",
776 : CPLAtof(SRS_UL_US_FOOT_CONV));
777 : else
778 0 : sr.SetStatePlane(iUTMZone, bNAD83);
779 : }
780 :
781 18 : m_oSRS = std::move(sr);
782 :
783 : /* -------------------------------------------------------------------- */
784 : /* For UTM we use the extents (really the UTM coordinates of */
785 : /* the lat/long corners of the quad) to determine the size in */
786 : /* pixels and lines, but we have to make the anchors be modulus */
787 : /* the pixel size which what really gets used. */
788 : /* -------------------------------------------------------------------- */
789 18 : if (nCoordSystem == 1 // UTM
790 6 : || nCoordSystem == 2 // State Plane
791 6 : || nCoordSystem == -9999) // unknown
792 : {
793 : // expand extents modulus the pixel size.
794 12 : extent_min.y = floor(extent_min.y / dydelta) * dydelta;
795 12 : extent_max.y = ceil(extent_max.y / dydelta) * dydelta;
796 :
797 : // Forcibly compute X extents based on first profile and pixelsize.
798 12 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, nDataStartOffset, 0));
799 12 : /* njunk = */ ReadInt(InDem);
800 12 : /* njunk = */ ReadInt(InDem);
801 12 : /* njunk = */ ReadInt(InDem);
802 12 : /* njunk = */ ReadInt(InDem);
803 12 : const double dxStart = DConvert(InDem, 24);
804 :
805 12 : nRasterYSize =
806 12 : static_cast<int>((extent_max.y - extent_min.y) / dydelta + 1.5);
807 12 : nRasterXSize = nProfiles;
808 :
809 12 : adfGeoTransform[0] = dxStart - dxdelta / 2.0;
810 12 : adfGeoTransform[1] = dxdelta;
811 12 : adfGeoTransform[2] = 0.0;
812 12 : adfGeoTransform[3] = extent_max.y + dydelta / 2.0;
813 12 : adfGeoTransform[4] = 0.0;
814 12 : adfGeoTransform[5] = -dydelta;
815 : }
816 : /* -------------------------------------------------------------------- */
817 : /* Geographic -- use corners directly. */
818 : /* -------------------------------------------------------------------- */
819 : else
820 : {
821 6 : nRasterYSize =
822 6 : static_cast<int>((extent_max.y - extent_min.y) / dydelta + 1.5);
823 6 : nRasterXSize = nProfiles;
824 :
825 : // Translate extents from arc-seconds to decimal degrees.
826 6 : adfGeoTransform[0] = (extent_min.x - dxdelta / 2.0) / 3600.0;
827 6 : adfGeoTransform[1] = dxdelta / 3600.0;
828 6 : adfGeoTransform[2] = 0.0;
829 6 : adfGeoTransform[3] = (extent_max.y + dydelta / 2.0) / 3600.0;
830 6 : adfGeoTransform[4] = 0.0;
831 6 : adfGeoTransform[5] = (-dydelta) / 3600.0;
832 : }
833 :
834 : // IReadBlock() not ready for more than INT_MAX pixels, and that
835 : // would behave badly
836 36 : if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize) ||
837 18 : nRasterXSize > INT_MAX / nRasterYSize)
838 : {
839 0 : return FALSE;
840 : }
841 :
842 18 : return TRUE;
843 : }
844 :
845 : /************************************************************************/
846 : /* GetGeoTransform() */
847 : /************************************************************************/
848 :
849 6 : CPLErr USGSDEMDataset::GetGeoTransform(double *padfTransform)
850 :
851 : {
852 6 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
853 6 : return CE_None;
854 : }
855 :
856 : /************************************************************************/
857 : /* GetSpatialRef() */
858 : /************************************************************************/
859 :
860 6 : const OGRSpatialReference *USGSDEMDataset::GetSpatialRef() const
861 :
862 : {
863 6 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
864 : }
865 :
866 : /************************************************************************/
867 : /* Identify() */
868 : /************************************************************************/
869 :
870 49463 : int USGSDEMDataset::Identify(GDALOpenInfo *poOpenInfo)
871 :
872 : {
873 49463 : if (poOpenInfo->nHeaderBytes < 200)
874 46542 : return FALSE;
875 :
876 2921 : if (!STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 0") &&
877 2909 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 1") &&
878 2885 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 2") &&
879 2885 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 3") &&
880 2885 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " -9999"))
881 2885 : return FALSE;
882 :
883 36 : if (!STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 150, " 1") &&
884 4 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 150, " 4"))
885 0 : return FALSE;
886 :
887 36 : return TRUE;
888 : }
889 :
890 : /************************************************************************/
891 : /* Open() */
892 : /************************************************************************/
893 :
894 18 : GDALDataset *USGSDEMDataset::Open(GDALOpenInfo *poOpenInfo)
895 :
896 : {
897 18 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
898 0 : return nullptr;
899 :
900 : /* -------------------------------------------------------------------- */
901 : /* Create a corresponding GDALDataset. */
902 : /* -------------------------------------------------------------------- */
903 18 : USGSDEMDataset *poDS = new USGSDEMDataset();
904 :
905 18 : poDS->fp = poOpenInfo->fpL;
906 18 : poOpenInfo->fpL = nullptr;
907 :
908 : /* -------------------------------------------------------------------- */
909 : /* Read the file. */
910 : /* -------------------------------------------------------------------- */
911 18 : if (!poDS->LoadFromFile(poDS->fp))
912 : {
913 0 : delete poDS;
914 0 : return nullptr;
915 : }
916 :
917 : /* -------------------------------------------------------------------- */
918 : /* Confirm the requested access is supported. */
919 : /* -------------------------------------------------------------------- */
920 18 : if (poOpenInfo->eAccess == GA_Update)
921 : {
922 0 : delete poDS;
923 0 : ReportUpdateNotSupportedByDriver("USGSDEM");
924 0 : return nullptr;
925 : }
926 :
927 : /* -------------------------------------------------------------------- */
928 : /* Create band information objects. */
929 : /* -------------------------------------------------------------------- */
930 18 : poDS->SetBand(1, new USGSDEMRasterBand(poDS));
931 :
932 18 : poDS->SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
933 :
934 : /* -------------------------------------------------------------------- */
935 : /* Initialize any PAM information. */
936 : /* -------------------------------------------------------------------- */
937 18 : poDS->SetDescription(poOpenInfo->pszFilename);
938 18 : poDS->TryLoadXML();
939 :
940 : /* -------------------------------------------------------------------- */
941 : /* Open overviews. */
942 : /* -------------------------------------------------------------------- */
943 18 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
944 :
945 18 : return poDS;
946 : }
947 :
948 : /************************************************************************/
949 : /* GDALRegister_USGSDEM() */
950 : /************************************************************************/
951 :
952 1667 : void GDALRegister_USGSDEM()
953 :
954 : {
955 1667 : if (GDALGetDriverByName("USGSDEM") != nullptr)
956 282 : return;
957 :
958 1385 : GDALDriver *poDriver = new GDALDriver();
959 :
960 1385 : poDriver->SetDescription("USGSDEM");
961 1385 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
962 1385 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "dem");
963 1385 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
964 1385 : "USGS Optional ASCII DEM (and CDED)");
965 1385 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
966 1385 : "drivers/raster/usgsdem.html");
967 :
968 1385 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
969 :
970 1385 : poDriver->pfnOpen = USGSDEMDataset::Open;
971 1385 : poDriver->pfnIdentify = USGSDEMDataset::Identify;
972 :
973 1385 : GetGDALDriverManager()->RegisterDriver(poDriver);
974 : }
|
