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