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 409 : static int ReadInt(VSILFILE *fp)
41 : {
42 : char c;
43 409 : int nRead = 0;
44 : char szBuffer[12];
45 409 : bool bInProlog = true;
46 :
47 : while (true)
48 : {
49 3833 : if (VSIFReadL(&c, 1, 1, fp) != 1)
50 : {
51 4 : return 0;
52 : }
53 3829 : if (bInProlog)
54 : {
55 3321 : if (!isspace(static_cast<unsigned char>(c)))
56 : {
57 405 : bInProlog = false;
58 : }
59 : }
60 3829 : if (!bInProlog)
61 : {
62 913 : if (c != '-' && c != '+' && !(c >= '0' && c <= '9'))
63 : {
64 405 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, VSIFTellL(fp) - 1, SEEK_SET));
65 405 : break;
66 : }
67 508 : if (nRead < 11)
68 508 : szBuffer[nRead] = c;
69 508 : nRead++;
70 : }
71 : }
72 405 : szBuffer[std::min(nRead, 11)] = 0;
73 405 : 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 53 : static void USGSDEMRefillBuffer(Buffer *psBuffer)
90 : {
91 53 : memmove(psBuffer->buffer, psBuffer->buffer + psBuffer->cur_index,
92 53 : psBuffer->buffer_size - psBuffer->cur_index);
93 :
94 53 : psBuffer->buffer_size -= psBuffer->cur_index;
95 53 : psBuffer->buffer_size += static_cast<int>(
96 106 : VSIFReadL(psBuffer->buffer + psBuffer->buffer_size, 1,
97 53 : psBuffer->max_size - psBuffer->buffer_size, psBuffer->fp));
98 53 : psBuffer->cur_index = 0;
99 53 : }
100 :
101 : /************************************************************************/
102 : /* USGSDEMGetCurrentFilePos() */
103 : /************************************************************************/
104 :
105 258 : static vsi_l_offset USGSDEMGetCurrentFilePos(const Buffer *psBuffer)
106 : {
107 258 : return VSIFTellL(psBuffer->fp) - psBuffer->buffer_size +
108 258 : psBuffer->cur_index;
109 : }
110 :
111 : /************************************************************************/
112 : /* USGSDEMSetCurrentFilePos() */
113 : /************************************************************************/
114 :
115 258 : static void USGSDEMSetCurrentFilePos(Buffer *psBuffer, vsi_l_offset nNewPos)
116 : {
117 258 : vsi_l_offset nCurPosFP = VSIFTellL(psBuffer->fp);
118 258 : if (nNewPos >= nCurPosFP - psBuffer->buffer_size && nNewPos < nCurPosFP)
119 : {
120 242 : psBuffer->cur_index =
121 242 : static_cast<int>(nNewPos - (nCurPosFP - psBuffer->buffer_size));
122 : }
123 : else
124 : {
125 16 : CPL_IGNORE_RET_VAL(VSIFSeekL(psBuffer->fp, nNewPos, SEEK_SET));
126 16 : psBuffer->buffer_size = 0;
127 16 : psBuffer->cur_index = 0;
128 : }
129 258 : }
130 :
131 : /************************************************************************/
132 : /* USGSDEMReadIntFromBuffer() */
133 : /************************************************************************/
134 :
135 1080060 : static int USGSDEMReadIntFromBuffer(Buffer *psBuffer, int *pbSuccess = nullptr)
136 : {
137 : char c;
138 :
139 : while (true)
140 : {
141 1080060 : if (psBuffer->cur_index >= psBuffer->buffer_size)
142 : {
143 49 : USGSDEMRefillBuffer(psBuffer);
144 49 : if (psBuffer->cur_index >= psBuffer->buffer_size)
145 : {
146 0 : if (pbSuccess)
147 0 : *pbSuccess = FALSE;
148 0 : return 0;
149 : }
150 : }
151 :
152 1080060 : c = psBuffer->buffer[psBuffer->cur_index];
153 1080060 : psBuffer->cur_index++;
154 1080060 : if (!isspace(static_cast<unsigned char>(c)))
155 46830 : break;
156 : }
157 :
158 46830 : GIntBig nVal = 0;
159 46830 : int nSign = 1;
160 46830 : if (c == '-')
161 6374 : nSign = -1;
162 40456 : else if (c == '+')
163 0 : nSign = 1;
164 40456 : else if (c >= '0' && c <= '9')
165 40456 : 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 136899 : 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 136899 : c = psBuffer->buffer[psBuffer->cur_index];
187 136899 : if (c >= '0' && c <= '9')
188 : {
189 90069 : psBuffer->cur_index++;
190 90069 : if (nVal * nSign < INT_MAX && nVal * nSign > INT_MIN)
191 : {
192 90069 : nVal = nVal * 10 + (c - '0');
193 90069 : if (nVal * nSign > INT_MAX)
194 : {
195 0 : nVal = INT_MAX;
196 0 : nSign = 1;
197 : }
198 90069 : else if (nVal * nSign < INT_MIN)
199 : {
200 0 : nVal = INT_MIN;
201 0 : nSign = 1;
202 : }
203 : }
204 : }
205 : else
206 : {
207 46830 : if (pbSuccess)
208 46830 : *pbSuccess = TRUE;
209 46830 : return static_cast<int>(nSign * nVal);
210 : }
211 : }
212 : }
213 :
214 : /************************************************************************/
215 : /* USGSDEMReadDoubleFromBuffer() */
216 : /************************************************************************/
217 :
218 6424 : static double USGSDEMReadDoubleFromBuffer(Buffer *psBuffer, int nCharCount,
219 : int *pbSuccess = nullptr)
220 :
221 : {
222 6424 : 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 6423 : char *szPtr = psBuffer->buffer + psBuffer->cur_index;
234 6423 : char backupC = szPtr[nCharCount];
235 6423 : szPtr[nCharCount] = 0;
236 160575 : for (int i = 0; i < nCharCount; i++)
237 : {
238 154152 : if (szPtr[i] == 'D')
239 6404 : szPtr[i] = 'E';
240 : }
241 :
242 6423 : double dfVal = CPLAtof(szPtr);
243 6423 : szPtr[nCharCount] = backupC;
244 6423 : psBuffer->cur_index += nCharCount;
245 :
246 6423 : if (pbSuccess)
247 6423 : *pbSuccess = TRUE;
248 6423 : return dfVal;
249 : }
250 :
251 : /************************************************************************/
252 : /* DConvert() */
253 : /************************************************************************/
254 :
255 497 : static double DConvert(VSILFILE *fp, int nCharCount)
256 :
257 : {
258 : char szBuffer[100];
259 :
260 497 : CPL_IGNORE_RET_VAL(VSIFReadL(szBuffer, nCharCount, 1, fp));
261 497 : szBuffer[nCharCount] = '\0';
262 :
263 12869 : for (int i = 0; i < nCharCount; i++)
264 : {
265 12372 : if (szBuffer[i] == 'D')
266 449 : szBuffer[i] = 'E';
267 : }
268 :
269 994 : 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 37 : USGSDEMRasterBand::USGSDEMRasterBand(USGSDEMDataset *poDSIn)
331 :
332 : {
333 37 : this->poDS = poDSIn;
334 37 : this->nBand = 1;
335 :
336 37 : eDataType = poDSIn->eNaturalDataFormat;
337 :
338 37 : nBlockXSize = poDSIn->GetRasterXSize();
339 37 : nBlockYSize = poDSIn->GetRasterYSize();
340 37 : }
341 :
342 : /************************************************************************/
343 : /* IReadBlock() */
344 : /************************************************************************/
345 :
346 14 : CPLErr USGSDEMRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
347 : CPL_UNUSED int nBlockYOff, void *pImage)
348 :
349 : {
350 : /* int bad = FALSE; */
351 14 : USGSDEMDataset *poGDS = reinterpret_cast<USGSDEMDataset *>(poDS);
352 :
353 : /* -------------------------------------------------------------------- */
354 : /* Initialize image buffer to nodata value. */
355 : /* -------------------------------------------------------------------- */
356 14 : GDALCopyWords(&USGSDEM_NODATA, GDT_Int32, 0, pImage, GetRasterDataType(),
357 : GDALGetDataTypeSizeBytes(GetRasterDataType()),
358 14 : GetXSize() * GetYSize());
359 :
360 : /* -------------------------------------------------------------------- */
361 : /* Seek to data. */
362 : /* -------------------------------------------------------------------- */
363 14 : CPL_IGNORE_RET_VAL(VSIFSeekL(poGDS->fp, poGDS->nDataStartOffset, 0));
364 :
365 14 : double dfYMin = poGDS->adfGeoTransform[3] +
366 14 : (GetYSize() - 0.5) * poGDS->adfGeoTransform[5];
367 :
368 : /* -------------------------------------------------------------------- */
369 : /* Read all the profiles into the image buffer. */
370 : /* -------------------------------------------------------------------- */
371 :
372 : Buffer sBuffer;
373 14 : sBuffer.max_size = 32768;
374 14 : sBuffer.buffer = reinterpret_cast<char *>(CPLMalloc(sBuffer.max_size + 1));
375 14 : sBuffer.fp = poGDS->fp;
376 14 : sBuffer.buffer_size = 0;
377 14 : sBuffer.cur_index = 0;
378 :
379 1298 : for (int i = 0; i < GetXSize(); i++)
380 : {
381 : int bSuccess;
382 1285 : const int nRowNumber = USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
383 1285 : if (nRowNumber != 1)
384 1 : CPLDebug("USGSDEM", "i = %d, nRowNumber = %d", i, nRowNumber);
385 1285 : if (bSuccess)
386 : {
387 : const int nColNumber =
388 1285 : USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
389 1285 : if (nColNumber != i + 1)
390 : {
391 5 : CPLDebug("USGSDEM", "i = %d, nColNumber = %d", i, nColNumber);
392 : }
393 : }
394 : const int nCPoints =
395 1285 : (bSuccess) ? USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess) : 0;
396 : #ifdef DEBUG_VERBOSE
397 : CPLDebug("USGSDEM", "i = %d, nCPoints = %d", i, nCPoints);
398 : #endif
399 :
400 1285 : if (bSuccess)
401 : {
402 : const int nNumberOfCols =
403 1285 : USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
404 1285 : if (nNumberOfCols != 1)
405 : {
406 0 : CPLDebug("USGSDEM", "i = %d, nNumberOfCols = %d", i,
407 : nNumberOfCols);
408 : }
409 : }
410 :
411 : // x-start
412 1285 : if (bSuccess)
413 1285 : /* dxStart = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24,
414 : &bSuccess);
415 :
416 : double dyStart =
417 1285 : (bSuccess) ? USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess)
418 1285 : : 0;
419 : const double dfElevOffset =
420 1285 : (bSuccess) ? USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess)
421 1285 : : 0;
422 :
423 : // min z value
424 1285 : if (bSuccess)
425 1285 : /* djunk = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess);
426 :
427 : // max z value
428 1285 : if (bSuccess)
429 1284 : /* djunk = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess);
430 1285 : if (!bSuccess)
431 : {
432 1 : CPLFree(sBuffer.buffer);
433 1 : return CE_Failure;
434 : }
435 :
436 1284 : if (poGDS->m_oSRS.IsGeographic())
437 246 : dyStart = dyStart / 3600.0;
438 :
439 1284 : double dygap = (dfYMin - dyStart) / poGDS->adfGeoTransform[5] + 0.5;
440 1284 : if (dygap <= INT_MIN || dygap >= INT_MAX || !std::isfinite(dygap))
441 : {
442 0 : CPLFree(sBuffer.buffer);
443 0 : return CE_Failure;
444 : }
445 1284 : int lygap = static_cast<int>(dygap);
446 1284 : if (nCPoints <= 0)
447 0 : continue;
448 1284 : if (lygap > INT_MAX - nCPoints)
449 0 : lygap = INT_MAX - nCPoints;
450 1284 : if (lygap < 0 && GetYSize() > INT_MAX + lygap)
451 : {
452 0 : CPLFree(sBuffer.buffer);
453 0 : return CE_Failure;
454 : }
455 :
456 42974 : for (int j = lygap; j < (nCPoints + lygap); j++)
457 : {
458 41690 : const int iY = GetYSize() - j - 1;
459 :
460 41690 : 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 41690 : if (!bSuccess)
467 : {
468 0 : CPLFree(sBuffer.buffer);
469 0 : return CE_Failure;
470 : }
471 :
472 41690 : if (iY < 0 || iY >= GetYSize())
473 : {
474 : /* bad = TRUE; */
475 : }
476 41690 : else if (nElev == USGSDEM_NODATA)
477 : /* leave in output buffer as nodata */;
478 : else
479 : {
480 36298 : const float fComputedElev =
481 36298 : static_cast<float>(nElev * poGDS->fVRes + dfElevOffset);
482 :
483 36298 : if (GetRasterDataType() == GDT_Int16)
484 : {
485 72474 : GUInt16 nVal = (fComputedElev < -32768) ? -32768
486 : : (fComputedElev > 32767)
487 : ? 32767
488 36237 : : static_cast<GInt16>(fComputedElev);
489 36237 : reinterpret_cast<GInt16 *>(pImage)[i + iY * GetXSize()] =
490 36237 : nVal;
491 : }
492 : else
493 : {
494 61 : reinterpret_cast<float *>(pImage)[i + iY * GetXSize()] =
495 : fComputedElev;
496 : }
497 : }
498 : }
499 :
500 1284 : 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 258 : vsi_l_offset nCurPos = USGSDEMGetCurrentFilePos(&sBuffer);
506 258 : vsi_l_offset nNewPos = (nCurPos + 1023) / 1024 * 1024;
507 258 : if (nNewPos > nCurPos)
508 : {
509 258 : USGSDEMSetCurrentFilePos(&sBuffer, nNewPos);
510 : }
511 : }
512 : }
513 13 : CPLFree(sBuffer.buffer);
514 :
515 13 : return CE_None;
516 : }
517 :
518 : /************************************************************************/
519 : /* GetNoDataValue() */
520 : /************************************************************************/
521 :
522 16 : double USGSDEMRasterBand::GetNoDataValue(int *pbSuccess)
523 :
524 : {
525 16 : if (pbSuccess != nullptr)
526 14 : *pbSuccess = TRUE;
527 :
528 16 : return USGSDEM_NODATA;
529 : }
530 :
531 : /************************************************************************/
532 : /* GetUnitType() */
533 : /************************************************************************/
534 9 : const char *USGSDEMRasterBand::GetUnitType()
535 : {
536 9 : USGSDEMDataset *poGDS = reinterpret_cast<USGSDEMDataset *>(poDS);
537 :
538 9 : return poGDS->pszUnits;
539 : }
540 :
541 : /************************************************************************/
542 : /* ==================================================================== */
543 : /* USGSDEMDataset */
544 : /* ==================================================================== */
545 : /************************************************************************/
546 :
547 : /************************************************************************/
548 : /* USGSDEMDataset() */
549 : /************************************************************************/
550 :
551 37 : USGSDEMDataset::USGSDEMDataset()
552 : : nDataStartOffset(0), eNaturalDataFormat(GDT_Unknown), fVRes(0.0),
553 37 : pszUnits(nullptr), fp(nullptr)
554 : {
555 37 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
556 37 : memset(adfGeoTransform, 0, sizeof(adfGeoTransform));
557 37 : }
558 :
559 : /************************************************************************/
560 : /* ~USGSDEMDataset() */
561 : /************************************************************************/
562 :
563 74 : USGSDEMDataset::~USGSDEMDataset()
564 :
565 : {
566 37 : FlushCache(true);
567 :
568 37 : if (fp != nullptr)
569 37 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
570 74 : }
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 37 : int USGSDEMDataset::LoadFromFile(VSILFILE *InDem)
583 : {
584 : // check for version of DEM format
585 37 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 864, 0));
586 :
587 : // Read DEM into matrix
588 37 : const int nRow = ReadInt(InDem);
589 37 : const int nColumn = ReadInt(InDem);
590 : const bool bNewFormat =
591 37 : VSIFTellL(InDem) >= 1024 || nRow != 1 || nColumn != 1;
592 37 : if (bNewFormat)
593 : {
594 37 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 1024, 0)); // New Format
595 37 : int i = ReadInt(InDem);
596 37 : int j = ReadInt(InDem);
597 37 : 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 33 : 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 33 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 1024, 0));
629 : char c;
630 66 : if (VSIFReadL(&c, 1, 1, InDem) == 1 && c == '\n' &&
631 2 : VSIFSeekL(InDem, 1024 + 1024 + 1, 0) == 0 &&
632 66 : VSIFReadL(&c, 1, 1, InDem) == 1 && c == '\n')
633 : {
634 2 : nDataStartOffset = 1025;
635 : }
636 : }
637 : }
638 : else
639 0 : nDataStartOffset = 864;
640 :
641 37 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 156, 0));
642 37 : const int nCoordSystem = ReadInt(InDem);
643 37 : const int iUTMZone = ReadInt(InDem);
644 :
645 37 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 528, 0));
646 37 : const int nGUnit = ReadInt(InDem);
647 37 : const int nVUnit = ReadInt(InDem);
648 :
649 : // Vertical Units in meters
650 37 : if (nVUnit == 1)
651 0 : pszUnits = "ft";
652 : else
653 37 : pszUnits = "m";
654 :
655 37 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 816, 0));
656 37 : const double dxdelta = DConvert(InDem, 12);
657 37 : const double dydelta = DConvert(InDem, 12);
658 37 : if (dydelta == 0)
659 0 : return FALSE;
660 37 : fVRes = DConvert(InDem, 12);
661 :
662 : /* -------------------------------------------------------------------- */
663 : /* Should we treat this as floating point, or GInt16. */
664 : /* -------------------------------------------------------------------- */
665 37 : if (nVUnit == 1 || fVRes < 1.0)
666 4 : eNaturalDataFormat = GDT_Float32;
667 : else
668 33 : eNaturalDataFormat = GDT_Int16;
669 :
670 : /* -------------------------------------------------------------------- */
671 : /* Read four corner coordinates. */
672 : /* -------------------------------------------------------------------- */
673 37 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 546, 0));
674 : DPoint2 corners[4]; // SW, NW, NE, SE
675 185 : for (int i = 0; i < 4; i++)
676 : {
677 148 : corners[i].x = DConvert(InDem, 24);
678 148 : corners[i].y = DConvert(InDem, 24);
679 : }
680 :
681 : // find absolute extents of raw vales
682 : DPoint2 extent_min, extent_max;
683 37 : extent_min.x = std::min(corners[0].x, corners[1].x);
684 37 : extent_max.x = std::max(corners[2].x, corners[3].x);
685 37 : extent_min.y = std::min(corners[0].y, corners[3].y);
686 37 : extent_max.y = std::max(corners[1].y, corners[2].y);
687 :
688 37 : /* dElevMin = */ DConvert(InDem, 48);
689 37 : /* dElevMax = */ DConvert(InDem, 48);
690 :
691 37 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 858, 0));
692 37 : const int nProfiles = ReadInt(InDem);
693 :
694 : /* -------------------------------------------------------------------- */
695 : /* Collect the spatial reference system. */
696 : /* -------------------------------------------------------------------- */
697 74 : OGRSpatialReference sr;
698 37 : sr.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
699 37 : bool bNAD83 = true;
700 :
701 : // OLD format header ends at byte 864
702 37 : if (bNewFormat)
703 : {
704 : // year of data compilation
705 37 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 876, 0));
706 : char szDateBuffer[5];
707 37 : 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 37 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 890, 0));
718 :
719 : char szHorzDatum[3];
720 37 : CPL_IGNORE_RET_VAL(VSIFReadL(szHorzDatum, 1, 2, InDem));
721 37 : szHorzDatum[2] = '\0';
722 37 : const int datum = atoi(szHorzDatum);
723 37 : switch (datum)
724 : {
725 4 : case 1:
726 4 : sr.SetWellKnownGeogCS("NAD27");
727 4 : bNAD83 = false;
728 4 : break;
729 :
730 6 : case 2:
731 6 : sr.SetWellKnownGeogCS("WGS72");
732 6 : break;
733 :
734 14 : case 3:
735 14 : sr.SetWellKnownGeogCS("WGS84");
736 14 : break;
737 :
738 3 : case 4:
739 3 : sr.SetWellKnownGeogCS("NAD83");
740 3 : 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 37 : if (nCoordSystem == 1) // UTM
757 : {
758 16 : if (iUTMZone >= -60 && iUTMZone <= 60)
759 : {
760 16 : sr.SetUTM(abs(iUTMZone), iUTMZone >= 0);
761 16 : 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 21 : 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 37 : 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 37 : if (nCoordSystem == 1 // UTM
790 21 : || nCoordSystem == 2 // State Plane
791 21 : || nCoordSystem == -9999) // unknown
792 : {
793 : // expand extents modulus the pixel size.
794 16 : extent_min.y = floor(extent_min.y / dydelta) * dydelta;
795 16 : extent_max.y = ceil(extent_max.y / dydelta) * dydelta;
796 :
797 : // Forcibly compute X extents based on first profile and pixelsize.
798 16 : CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, nDataStartOffset, 0));
799 16 : /* njunk = */ ReadInt(InDem);
800 16 : /* njunk = */ ReadInt(InDem);
801 16 : /* njunk = */ ReadInt(InDem);
802 16 : /* njunk = */ ReadInt(InDem);
803 16 : const double dxStart = DConvert(InDem, 24);
804 :
805 16 : nRasterYSize =
806 16 : static_cast<int>((extent_max.y - extent_min.y) / dydelta + 1.5);
807 16 : nRasterXSize = nProfiles;
808 :
809 16 : adfGeoTransform[0] = dxStart - dxdelta / 2.0;
810 16 : adfGeoTransform[1] = dxdelta;
811 16 : adfGeoTransform[2] = 0.0;
812 16 : adfGeoTransform[3] = extent_max.y + dydelta / 2.0;
813 16 : adfGeoTransform[4] = 0.0;
814 16 : adfGeoTransform[5] = -dydelta;
815 : }
816 : /* -------------------------------------------------------------------- */
817 : /* Geographic -- use corners directly. */
818 : /* -------------------------------------------------------------------- */
819 : else
820 : {
821 21 : nRasterYSize =
822 21 : static_cast<int>((extent_max.y - extent_min.y) / dydelta + 1.5);
823 21 : nRasterXSize = nProfiles;
824 :
825 : // Translate extents from arc-seconds to decimal degrees.
826 21 : adfGeoTransform[0] = (extent_min.x - dxdelta / 2.0) / 3600.0;
827 21 : adfGeoTransform[1] = dxdelta / 3600.0;
828 21 : adfGeoTransform[2] = 0.0;
829 21 : adfGeoTransform[3] = (extent_max.y + dydelta / 2.0) / 3600.0;
830 21 : adfGeoTransform[4] = 0.0;
831 21 : adfGeoTransform[5] = (-dydelta) / 3600.0;
832 : }
833 :
834 : // IReadBlock() not ready for more than INT_MAX pixels, and that
835 : // would behave badly
836 74 : if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize) ||
837 37 : nRasterXSize > INT_MAX / nRasterYSize)
838 : {
839 0 : return FALSE;
840 : }
841 :
842 37 : return TRUE;
843 : }
844 :
845 : /************************************************************************/
846 : /* GetGeoTransform() */
847 : /************************************************************************/
848 :
849 33 : CPLErr USGSDEMDataset::GetGeoTransform(double *padfTransform)
850 :
851 : {
852 33 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
853 33 : return CE_None;
854 : }
855 :
856 : /************************************************************************/
857 : /* GetSpatialRef() */
858 : /************************************************************************/
859 :
860 32 : const OGRSpatialReference *USGSDEMDataset::GetSpatialRef() const
861 :
862 : {
863 32 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
864 : }
865 :
866 : /************************************************************************/
867 : /* Identify() */
868 : /************************************************************************/
869 :
870 51321 : int USGSDEMDataset::Identify(GDALOpenInfo *poOpenInfo)
871 :
872 : {
873 51321 : if (poOpenInfo->nHeaderBytes < 200)
874 48389 : return FALSE;
875 :
876 2932 : if (!STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 0") &&
877 2890 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 1") &&
878 2858 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 2") &&
879 2858 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " 3") &&
880 2858 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " -9999"))
881 2858 : return FALSE;
882 :
883 74 : if (!STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 150, " 1") &&
884 6 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 150, " 4"))
885 0 : return FALSE;
886 :
887 74 : return TRUE;
888 : }
889 :
890 : /************************************************************************/
891 : /* Open() */
892 : /************************************************************************/
893 :
894 37 : GDALDataset *USGSDEMDataset::Open(GDALOpenInfo *poOpenInfo)
895 :
896 : {
897 37 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
898 0 : return nullptr;
899 :
900 : /* -------------------------------------------------------------------- */
901 : /* Create a corresponding GDALDataset. */
902 : /* -------------------------------------------------------------------- */
903 37 : USGSDEMDataset *poDS = new USGSDEMDataset();
904 :
905 37 : poDS->fp = poOpenInfo->fpL;
906 37 : poOpenInfo->fpL = nullptr;
907 :
908 : /* -------------------------------------------------------------------- */
909 : /* Read the file. */
910 : /* -------------------------------------------------------------------- */
911 37 : 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 37 : if (poOpenInfo->eAccess == GA_Update)
921 : {
922 0 : delete poDS;
923 0 : CPLError(CE_Failure, CPLE_NotSupported,
924 : "The USGSDEM driver does not support update access to existing"
925 : " datasets.\n");
926 0 : return nullptr;
927 : }
928 :
929 : /* -------------------------------------------------------------------- */
930 : /* Create band information objects. */
931 : /* -------------------------------------------------------------------- */
932 37 : poDS->SetBand(1, new USGSDEMRasterBand(poDS));
933 :
934 37 : poDS->SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
935 :
936 : /* -------------------------------------------------------------------- */
937 : /* Initialize any PAM information. */
938 : /* -------------------------------------------------------------------- */
939 37 : poDS->SetDescription(poOpenInfo->pszFilename);
940 37 : poDS->TryLoadXML();
941 :
942 : /* -------------------------------------------------------------------- */
943 : /* Open overviews. */
944 : /* -------------------------------------------------------------------- */
945 37 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
946 :
947 37 : return poDS;
948 : }
949 :
950 : /************************************************************************/
951 : /* GDALRegister_USGSDEM() */
952 : /************************************************************************/
953 :
954 1595 : void GDALRegister_USGSDEM()
955 :
956 : {
957 1595 : if (GDALGetDriverByName("USGSDEM") != nullptr)
958 302 : return;
959 :
960 1293 : GDALDriver *poDriver = new GDALDriver();
961 :
962 1293 : poDriver->SetDescription("USGSDEM");
963 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
964 1293 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "dem");
965 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
966 1293 : "USGS Optional ASCII DEM (and CDED)");
967 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
968 1293 : "drivers/raster/usgsdem.html");
969 1293 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Int16");
970 1293 : poDriver->SetMetadataItem(
971 : GDAL_DMD_CREATIONOPTIONLIST,
972 : "<CreationOptionList>"
973 : " <Option name='PRODUCT' type='string-select' description='Specific "
974 : "Product Type'>"
975 : " <Value>DEFAULT</Value>"
976 : " <Value>CDED50K</Value>"
977 : " </Option>"
978 : " <Option name='TOPLEFT' type='string' description='Top left product "
979 : "corner (i.e. 117d15w,52d30n'/>"
980 : " <Option name='RESAMPLE' type='string-select' "
981 : "description='Resampling kernel to use if resampled.'>"
982 : " <Value>Nearest</Value>"
983 : " <Value>Bilinear</Value>"
984 : " <Value>Cubic</Value>"
985 : " <Value>CubicSpline</Value>"
986 : " </Option>"
987 : " <Option name='TEMPLATE' type='string' description='File to default "
988 : "metadata from.'/>"
989 : " <Option name='DEMLevelCode' type='int' description='DEM Level (1, "
990 : "2 or 3 if set)'/>"
991 : " <Option name='DataSpecVersion' type='int' description='Data and "
992 : "Specification version/revision (eg. 1020)'/>"
993 : " <Option name='PRODUCER' type='string' description='Producer Agency "
994 : "(up to 60 characters)'/>"
995 : " <Option name='OriginCode' type='string' description='Origin code "
996 : "(up to 4 characters, YT for Yukon)'/>"
997 : " <Option name='ProcessCode' type='string' description='Processing "
998 : "Code (8=ANUDEM, 9=FME, A=TopoGrid)'/>"
999 : " <Option name='ZRESOLUTION' type='float' description='Scaling "
1000 : "factor for elevation values'/>"
1001 : " <Option name='NTS' type='string' description='NTS Mapsheet name, "
1002 : "used to derive TOPLEFT.'/>"
1003 : " <Option name='INTERNALNAME' type='string' description='Dataset "
1004 : "name written into file header.'/>"
1005 1293 : "</CreationOptionList>");
1006 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1007 :
1008 1293 : poDriver->pfnOpen = USGSDEMDataset::Open;
1009 1293 : poDriver->pfnCreateCopy = USGSDEMCreateCopy;
1010 1293 : poDriver->pfnIdentify = USGSDEMDataset::Identify;
1011 :
1012 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1013 : }
|
