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