1 : /******************************************************************************
2 : * $Id: usgsdemdataset.cpp 25675 2013-02-23 15:33:45Z rouault $
3 : *
4 : * Project: USGS DEM Driver
5 : * Purpose: All reader for USGS DEM Reader
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : * Portions of this module derived from the VTP USGS DEM driver by Ben
9 : * Discoe, see http://www.vterrain.org
10 : *
11 : ******************************************************************************
12 : * Copyright (c) 2001, Frank Warmerdam <warmerdam@pobox.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_pam.h"
34 : #include "ogr_spatialref.h"
35 :
36 : CPL_CVSID("$Id: usgsdemdataset.cpp 25675 2013-02-23 15:33:45Z rouault $");
37 :
38 : CPL_C_START
39 : void GDALRegister_USGSDEM(void);
40 : CPL_C_END
41 :
42 : typedef struct {
43 : double x;
44 : double y;
45 : } DPoint2;
46 :
47 : #define USGSDEM_NODATA -32767
48 :
49 : GDALDataset *USGSDEMCreateCopy( const char *, GDALDataset *, int, char **,
50 : GDALProgressFunc pfnProgress,
51 : void * pProgressData );
52 :
53 :
54 : /************************************************************************/
55 : /* ReadInt() */
56 : /************************************************************************/
57 :
58 240 : static int ReadInt( VSILFILE *fp )
59 : {
60 240 : int nVal = 0;
61 : char c;
62 240 : int nRead = 0;
63 240 : vsi_l_offset nOffset = VSIFTellL(fp);
64 :
65 1604 : while (TRUE)
66 : {
67 1844 : if (VSIFReadL(&c, 1, 1, fp) != 1)
68 : {
69 0 : return 0;
70 : }
71 : else
72 1844 : nRead ++;
73 1844 : if (!isspace((int)c))
74 : break;
75 : }
76 :
77 240 : int nSign = 1;
78 240 : if (c == '-')
79 0 : nSign = -1;
80 240 : else if (c == '+')
81 0 : nSign = 1;
82 479 : else if (c >= '0' && c <= '9')
83 239 : nVal = c - '0';
84 : else
85 : {
86 1 : VSIFSeekL(fp, nOffset + nRead, SEEK_SET);
87 1 : return 0;
88 : }
89 :
90 45 : while (TRUE)
91 : {
92 284 : if (VSIFReadL(&c, 1, 1, fp) != 1)
93 0 : return nSign * nVal;
94 284 : nRead ++;
95 284 : if (c >= '0' && c <= '9')
96 45 : nVal = nVal * 10 + (c - '0');
97 : else
98 : {
99 239 : VSIFSeekL(fp, nOffset + (nRead - 1), SEEK_SET);
100 239 : return nSign * nVal;
101 : }
102 : }
103 : }
104 :
105 : typedef struct
106 : {
107 : VSILFILE *fp;
108 : int max_size;
109 : char* buffer;
110 : int buffer_size;
111 : int cur_index;
112 : } Buffer;
113 :
114 : /************************************************************************/
115 : /* USGSDEMRefillBuffer() */
116 : /************************************************************************/
117 :
118 16 : static void USGSDEMRefillBuffer( Buffer* psBuffer )
119 : {
120 : memmove(psBuffer->buffer, psBuffer->buffer + psBuffer->cur_index,
121 16 : psBuffer->buffer_size - psBuffer->cur_index);
122 :
123 16 : psBuffer->buffer_size -= psBuffer->cur_index;
124 : psBuffer->buffer_size += VSIFReadL(psBuffer->buffer + psBuffer->buffer_size,
125 : 1, psBuffer->max_size - psBuffer->buffer_size,
126 16 : psBuffer->fp);
127 16 : psBuffer->cur_index = 0;
128 16 : }
129 :
130 : /************************************************************************/
131 : /* USGSDEMReadIntFromBuffer() */
132 : /************************************************************************/
133 :
134 40262 : static int USGSDEMReadIntFromBuffer( Buffer* psBuffer, int* pbSuccess = NULL )
135 : {
136 40262 : int nVal = 0;
137 : char c;
138 :
139 156438 : while (TRUE)
140 : {
141 196700 : if (psBuffer->cur_index >= psBuffer->buffer_size)
142 : {
143 16 : USGSDEMRefillBuffer(psBuffer);
144 16 : if (psBuffer->cur_index >= psBuffer->buffer_size)
145 : {
146 0 : if( pbSuccess ) *pbSuccess = FALSE;
147 0 : return 0;
148 : }
149 : }
150 :
151 196700 : c = psBuffer->buffer[psBuffer->cur_index];
152 196700 : psBuffer->cur_index ++;
153 196700 : if (!isspace((int)c))
154 : break;
155 : }
156 :
157 40262 : int nSign = 1;
158 40262 : if (c == '-')
159 6192 : nSign = -1;
160 34070 : else if (c == '+')
161 0 : nSign = 1;
162 68140 : else if (c >= '0' && c <= '9')
163 34070 : nVal = c - '0';
164 : else
165 : {
166 0 : if( pbSuccess ) *pbSuccess = FALSE;
167 0 : return 0;
168 : }
169 :
170 84005 : while (TRUE)
171 : {
172 124267 : if (psBuffer->cur_index >= psBuffer->buffer_size)
173 : {
174 0 : USGSDEMRefillBuffer(psBuffer);
175 0 : if (psBuffer->cur_index >= psBuffer->buffer_size)
176 : {
177 0 : if( pbSuccess ) *pbSuccess = TRUE;
178 0 : return nSign * nVal;
179 : }
180 : }
181 :
182 124267 : c = psBuffer->buffer[psBuffer->cur_index];
183 124267 : if (c >= '0' && c <= '9')
184 : {
185 84005 : psBuffer->cur_index ++;
186 84005 : nVal = nVal * 10 + (c - '0');
187 : }
188 : else
189 : {
190 40262 : if( pbSuccess ) *pbSuccess = TRUE;
191 40262 : return nSign * nVal;
192 : }
193 : }
194 : }
195 :
196 : /************************************************************************/
197 : /* USGSDEMReadDoubleFromBuffer() */
198 : /************************************************************************/
199 :
200 1280 : static double USGSDEMReadDoubleFromBuffer( Buffer* psBuffer, int nCharCount )
201 :
202 : {
203 : int i;
204 :
205 1280 : if (psBuffer->cur_index + nCharCount > psBuffer->buffer_size)
206 : {
207 0 : USGSDEMRefillBuffer(psBuffer);
208 0 : if (psBuffer->cur_index + nCharCount > psBuffer->buffer_size)
209 : {
210 0 : return 0;
211 : }
212 : }
213 :
214 1280 : char* szPtr = psBuffer->buffer + psBuffer->cur_index;
215 1280 : char backupC = szPtr[nCharCount];
216 1280 : szPtr[nCharCount] = 0;
217 32000 : for( i = 0; i < nCharCount; i++ )
218 : {
219 30720 : if( szPtr[i] == 'D' )
220 1266 : szPtr[i] = 'E';
221 : }
222 :
223 1280 : double dfVal = CPLAtof(szPtr);
224 1280 : szPtr[nCharCount] = backupC;
225 1280 : psBuffer->cur_index += nCharCount;
226 :
227 1280 : return dfVal;
228 : }
229 :
230 : /************************************************************************/
231 : /* DConvert() */
232 : /************************************************************************/
233 :
234 318 : static double DConvert( VSILFILE *fp, int nCharCount )
235 :
236 : {
237 : char szBuffer[100];
238 : int i;
239 :
240 318 : VSIFReadL( szBuffer, nCharCount, 1, fp );
241 318 : szBuffer[nCharCount] = '\0';
242 :
243 8238 : for( i = 0; i < nCharCount; i++ )
244 : {
245 7920 : if( szBuffer[i] == 'D' )
246 309 : szBuffer[i] = 'E';
247 : }
248 :
249 318 : return CPLAtof(szBuffer);
250 : }
251 :
252 : /************************************************************************/
253 : /* ==================================================================== */
254 : /* USGSDEMDataset */
255 : /* ==================================================================== */
256 : /************************************************************************/
257 :
258 : class USGSDEMRasterBand;
259 :
260 : class USGSDEMDataset : public GDALPamDataset
261 : {
262 : friend class USGSDEMRasterBand;
263 :
264 : int nDataStartOffset;
265 : GDALDataType eNaturalDataFormat;
266 :
267 : double adfGeoTransform[6];
268 : char *pszProjection;
269 :
270 : double fVRes;
271 :
272 : const char *pszUnits;
273 :
274 : int LoadFromFile( VSILFILE * );
275 :
276 : VSILFILE *fp;
277 :
278 : public:
279 : USGSDEMDataset();
280 : ~USGSDEMDataset();
281 :
282 : static int Identify( GDALOpenInfo * );
283 : static GDALDataset *Open( GDALOpenInfo * );
284 : CPLErr GetGeoTransform( double * padfTransform );
285 : const char *GetProjectionRef();
286 : };
287 :
288 : /************************************************************************/
289 : /* ==================================================================== */
290 : /* USGSDEMRasterBand */
291 : /* ==================================================================== */
292 : /************************************************************************/
293 :
294 : class USGSDEMRasterBand : public GDALPamRasterBand
295 24 : {
296 : friend class USGSDEMDataset;
297 :
298 : public:
299 :
300 : USGSDEMRasterBand( USGSDEMDataset * );
301 :
302 : virtual const char *GetUnitType();
303 : virtual double GetNoDataValue( int *pbSuccess = NULL );
304 : virtual CPLErr IReadBlock( int, int, void * );
305 : };
306 :
307 :
308 : /************************************************************************/
309 : /* USGSDEMRasterBand() */
310 : /************************************************************************/
311 :
312 24 : USGSDEMRasterBand::USGSDEMRasterBand( USGSDEMDataset *poDS )
313 :
314 : {
315 24 : this->poDS = poDS;
316 24 : this->nBand = 1;
317 :
318 24 : eDataType = poDS->eNaturalDataFormat;
319 :
320 24 : nBlockXSize = poDS->GetRasterXSize();
321 24 : nBlockYSize = poDS->GetRasterYSize();
322 :
323 24 : }
324 :
325 : /************************************************************************/
326 : /* IReadBlock() */
327 : /************************************************************************/
328 :
329 10 : CPLErr USGSDEMRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
330 : void * pImage )
331 :
332 : {
333 : double dfYMin;
334 10 : int bad = FALSE;
335 10 : USGSDEMDataset *poGDS = (USGSDEMDataset *) poDS;
336 :
337 : /* -------------------------------------------------------------------- */
338 : /* Initialize image buffer to nodata value. */
339 : /* -------------------------------------------------------------------- */
340 40678 : for( int k = GetXSize() * GetYSize() - 1; k >= 0; k-- )
341 : {
342 40668 : if( GetRasterDataType() == GDT_Int16 )
343 37846 : ((GInt16 *) pImage)[k] = USGSDEM_NODATA;
344 : else
345 2822 : ((float *) pImage)[k] = USGSDEM_NODATA;
346 : }
347 :
348 : /* -------------------------------------------------------------------- */
349 : /* Seek to data. */
350 : /* -------------------------------------------------------------------- */
351 10 : VSIFSeekL(poGDS->fp, poGDS->nDataStartOffset, 0);
352 :
353 10 : dfYMin = poGDS->adfGeoTransform[3]
354 10 : + (GetYSize()-0.5) * poGDS->adfGeoTransform[5];
355 :
356 : /* -------------------------------------------------------------------- */
357 : /* Read all the profiles into the image buffer. */
358 : /* -------------------------------------------------------------------- */
359 :
360 : Buffer sBuffer;
361 10 : sBuffer.max_size = 32768;
362 10 : sBuffer.buffer = (char*) CPLMalloc(sBuffer.max_size + 1);
363 10 : sBuffer.fp = poGDS->fp;
364 10 : sBuffer.buffer_size = 0;
365 10 : sBuffer.cur_index = 0;
366 :
367 266 : for( int i = 0; i < GetXSize(); i++)
368 : {
369 : int njunk, nCPoints, lygap;
370 : double djunk, dxStart, dyStart, dfElevOffset;
371 :
372 256 : njunk = USGSDEMReadIntFromBuffer(&sBuffer);
373 256 : njunk = USGSDEMReadIntFromBuffer(&sBuffer);
374 256 : nCPoints = USGSDEMReadIntFromBuffer(&sBuffer);
375 256 : njunk = USGSDEMReadIntFromBuffer(&sBuffer);
376 :
377 256 : dxStart = USGSDEMReadDoubleFromBuffer(&sBuffer, 24);
378 256 : dyStart = USGSDEMReadDoubleFromBuffer(&sBuffer, 24);
379 256 : dfElevOffset = USGSDEMReadDoubleFromBuffer(&sBuffer, 24);
380 256 : djunk = USGSDEMReadDoubleFromBuffer(&sBuffer, 24);
381 256 : djunk = USGSDEMReadDoubleFromBuffer(&sBuffer, 24);
382 :
383 256 : if( EQUALN(poGDS->pszProjection,"GEOGCS",6) )
384 246 : dyStart = dyStart / 3600.0;
385 :
386 256 : lygap = (int)((dfYMin - dyStart)/poGDS->adfGeoTransform[5]+ 0.5);
387 :
388 39494 : for (int j=lygap; j < (nCPoints+(int)lygap); j++)
389 : {
390 39238 : int iY = GetYSize() - j - 1;
391 : int nElev;
392 : int bSuccess;
393 :
394 39238 : nElev = USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
395 39238 : if( !bSuccess )
396 : {
397 0 : CPLFree(sBuffer.buffer);
398 0 : return CE_Failure;
399 : }
400 :
401 39238 : if (iY < 0 || iY >= GetYSize() )
402 0 : bad = TRUE;
403 39238 : else if( nElev == USGSDEM_NODATA )
404 : /* leave in output buffer as nodata */;
405 : else
406 : {
407 : float fComputedElev =
408 33846 : (float)(nElev * poGDS->fVRes + dfElevOffset);
409 :
410 33846 : if( GetRasterDataType() == GDT_Int16 )
411 : {
412 33785 : ((GInt16 *) pImage)[i + iY*GetXSize()] =
413 33785 : (GInt16) fComputedElev;
414 : }
415 : else
416 : {
417 61 : ((float *) pImage)[i + iY*GetXSize()] = fComputedElev;
418 : }
419 : }
420 : }
421 : }
422 10 : CPLFree(sBuffer.buffer);
423 :
424 10 : return CE_None;
425 : }
426 :
427 : /************************************************************************/
428 : /* GetNoDataValue() */
429 : /************************************************************************/
430 :
431 16 : double USGSDEMRasterBand::GetNoDataValue( int *pbSuccess )
432 :
433 : {
434 16 : if( pbSuccess != NULL )
435 14 : *pbSuccess = TRUE;
436 :
437 16 : return USGSDEM_NODATA;
438 : }
439 :
440 : /************************************************************************/
441 : /* GetUnitType() */
442 : /************************************************************************/
443 9 : const char *USGSDEMRasterBand::GetUnitType()
444 : {
445 9 : USGSDEMDataset *poGDS = (USGSDEMDataset *) poDS;
446 :
447 9 : return poGDS->pszUnits;
448 : }
449 :
450 : /************************************************************************/
451 : /* ==================================================================== */
452 : /* USGSDEMDataset */
453 : /* ==================================================================== */
454 : /************************************************************************/
455 :
456 : /************************************************************************/
457 : /* USGSDEMDataset() */
458 : /************************************************************************/
459 :
460 24 : USGSDEMDataset::USGSDEMDataset()
461 :
462 : {
463 24 : fp = NULL;
464 24 : pszProjection = NULL;
465 24 : }
466 :
467 : /************************************************************************/
468 : /* ~USGSDEMDataset() */
469 : /************************************************************************/
470 :
471 24 : USGSDEMDataset::~USGSDEMDataset()
472 :
473 : {
474 24 : FlushCache();
475 :
476 24 : CPLFree( pszProjection );
477 24 : if( fp != NULL )
478 24 : VSIFCloseL( fp );
479 24 : }
480 :
481 : /************************************************************************/
482 : /* LoadFromFile() */
483 : /* */
484 : /* If the data from DEM is in meters, then values are stored as */
485 : /* shorts. If DEM data is in feet, then height data will be */
486 : /* stored in float, to preserve the precision of the original */
487 : /* data. returns true if the file was successfully opened and */
488 : /* read. */
489 : /************************************************************************/
490 :
491 24 : int USGSDEMDataset::LoadFromFile(VSILFILE *InDem)
492 : {
493 : int i, j;
494 : int nRow, nColumn;
495 : int nVUnit, nGUnit;
496 : double dxdelta, dydelta;
497 : double dElevMax, dElevMin;
498 : int bNewFormat;
499 : int nCoordSystem;
500 : int nProfiles;
501 : char szDateBuffer[5];
502 : DPoint2 corners[4]; // SW, NW, NE, SE
503 : DPoint2 extent_min, extent_max;
504 : int iUTMZone;
505 :
506 : // check for version of DEM format
507 24 : VSIFSeekL(InDem, 864, 0);
508 :
509 : // Read DEM into matrix
510 24 : nRow = ReadInt(InDem);
511 24 : nColumn = ReadInt(InDem);
512 24 : bNewFormat = ((nRow!=1)||(nColumn!=1));
513 24 : if (bNewFormat)
514 : {
515 23 : VSIFSeekL(InDem, 1024, 0); // New Format
516 23 : i = ReadInt(InDem);
517 23 : j = ReadInt(InDem);
518 24 : if ((i!=1)||(j!=1 && j != 0)) // File OK?
519 : {
520 1 : VSIFSeekL(InDem, 893, 0); // Undocumented Format (39109h1.dem)
521 1 : i = ReadInt(InDem);
522 1 : j = ReadInt(InDem);
523 1 : if ((i!=1)||(j!=1)) // File OK?
524 : {
525 : CPLError( CE_Failure, CPLE_AppDefined,
526 0 : "Does not appear to be a USGS DEM file." );
527 0 : return FALSE;
528 : }
529 : else
530 1 : nDataStartOffset = 893;
531 : }
532 : else
533 22 : nDataStartOffset = 1024;
534 : }
535 : else
536 1 : nDataStartOffset = 864;
537 :
538 24 : VSIFSeekL(InDem, 156, 0);
539 24 : nCoordSystem = ReadInt(InDem);
540 24 : iUTMZone = ReadInt(InDem);
541 :
542 24 : VSIFSeekL(InDem, 528, 0);
543 24 : nGUnit = ReadInt(InDem);
544 24 : nVUnit = ReadInt(InDem);
545 :
546 : // Vertical Units in meters
547 24 : if (nVUnit==1)
548 0 : pszUnits = "ft";
549 : else
550 24 : pszUnits = "m";
551 :
552 24 : VSIFSeekL(InDem, 816, 0);
553 24 : dxdelta = DConvert(InDem, 12);
554 24 : dydelta = DConvert(InDem, 12);
555 24 : fVRes = DConvert(InDem, 12);
556 :
557 : /* -------------------------------------------------------------------- */
558 : /* Should we treat this as floating point, or GInt16. */
559 : /* -------------------------------------------------------------------- */
560 25 : if (nVUnit==1 || fVRes < 1.0)
561 1 : eNaturalDataFormat = GDT_Float32;
562 : else
563 23 : eNaturalDataFormat = GDT_Int16;
564 :
565 : /* -------------------------------------------------------------------- */
566 : /* Read four corner coordinates. */
567 : /* -------------------------------------------------------------------- */
568 24 : VSIFSeekL(InDem, 546, 0);
569 120 : for (i = 0; i < 4; i++)
570 : {
571 96 : corners[i].x = DConvert(InDem, 24);
572 96 : corners[i].y = DConvert(InDem, 24);
573 : }
574 :
575 : // find absolute extents of raw vales
576 24 : extent_min.x = MIN(corners[0].x, corners[1].x);
577 24 : extent_max.x = MAX(corners[2].x, corners[3].x);
578 24 : extent_min.y = MIN(corners[0].y, corners[3].y);
579 24 : extent_max.y = MAX(corners[1].y, corners[2].y);
580 :
581 24 : dElevMin = DConvert(InDem, 48);
582 24 : dElevMax = DConvert(InDem, 48);
583 :
584 24 : VSIFSeekL(InDem, 858, 0);
585 24 : nProfiles = ReadInt(InDem);
586 :
587 : /* -------------------------------------------------------------------- */
588 : /* Collect the spatial reference system. */
589 : /* -------------------------------------------------------------------- */
590 24 : OGRSpatialReference sr;
591 24 : int bNAD83 =TRUE;
592 :
593 : // OLD format header ends at byte 864
594 24 : if (bNewFormat)
595 : {
596 : char szHorzDatum[3];
597 :
598 : // year of data compilation
599 23 : VSIFSeekL(InDem, 876, 0);
600 23 : VSIFReadL(szDateBuffer, 4, 1, InDem);
601 23 : szDateBuffer[4] = 0;
602 :
603 : // Horizontal datum
604 : // 1=North American Datum 1927 (NAD 27)
605 : // 2=World Geodetic System 1972 (WGS 72)
606 : // 3=WGS 84
607 : // 4=NAD 83
608 : // 5=Old Hawaii Datum
609 : // 6=Puerto Rico Datum
610 : int datum;
611 23 : VSIFSeekL(InDem, 890, 0);
612 23 : VSIFReadL( szHorzDatum, 1, 2, InDem );
613 23 : szHorzDatum[2] = '\0';
614 23 : datum = atoi(szHorzDatum);
615 23 : switch (datum)
616 : {
617 : case 1:
618 1 : sr.SetWellKnownGeogCS( "NAD27" );
619 1 : bNAD83 = FALSE;
620 1 : break;
621 :
622 : case 2:
623 5 : sr.SetWellKnownGeogCS( "WGS72" );
624 5 : break;
625 :
626 : case 3:
627 14 : sr.SetWellKnownGeogCS( "WGS84" );
628 14 : break;
629 :
630 : case 4:
631 1 : sr.SetWellKnownGeogCS( "NAD83" );
632 1 : break;
633 :
634 : case -9:
635 0 : break;
636 :
637 : default:
638 2 : sr.SetWellKnownGeogCS( "NAD27" );
639 : break;
640 : }
641 : }
642 : else
643 : {
644 1 : sr.SetWellKnownGeogCS( "NAD27" );
645 1 : bNAD83 = FALSE;
646 : }
647 :
648 24 : if (nCoordSystem == 1) // UTM
649 6 : sr.SetUTM( iUTMZone, TRUE );
650 :
651 18 : else if (nCoordSystem == 2) // state plane
652 : {
653 0 : if( nGUnit == 1 )
654 : sr.SetStatePlane( iUTMZone, bNAD83,
655 0 : "Foot", CPLAtof(SRS_UL_US_FOOT_CONV) );
656 : else
657 0 : sr.SetStatePlane( iUTMZone, bNAD83 );
658 : }
659 :
660 24 : sr.exportToWkt( &pszProjection );
661 :
662 : /* -------------------------------------------------------------------- */
663 : /* For UTM we use the extents (really the UTM coordinates of */
664 : /* the lat/long corners of the quad) to determine the size in */
665 : /* pixels and lines, but we have to make the anchors be modulus */
666 : /* the pixel size which what really gets used. */
667 : /* -------------------------------------------------------------------- */
668 30 : if (nCoordSystem == 1 // UTM
669 : || nCoordSystem == 2 // State Plane
670 : || nCoordSystem == -9999 ) // unknown
671 : {
672 : int njunk;
673 : double dxStart;
674 :
675 : // expand extents modulus the pixel size.
676 6 : extent_min.y = floor(extent_min.y/dydelta) * dydelta;
677 6 : extent_max.y = ceil(extent_max.y/dydelta) * dydelta;
678 :
679 : // Forceably compute X extents based on first profile and pixelsize.
680 6 : VSIFSeekL(InDem, nDataStartOffset, 0);
681 6 : njunk = ReadInt(InDem);
682 6 : njunk = ReadInt(InDem);
683 6 : njunk = ReadInt(InDem);
684 6 : njunk = ReadInt(InDem);
685 6 : dxStart = DConvert(InDem, 24);
686 :
687 6 : nRasterYSize = (int) ((extent_max.y - extent_min.y)/dydelta + 1.5);
688 6 : nRasterXSize = nProfiles;
689 :
690 6 : adfGeoTransform[0] = dxStart - dxdelta/2.0;
691 6 : adfGeoTransform[1] = dxdelta;
692 6 : adfGeoTransform[2] = 0.0;
693 6 : adfGeoTransform[3] = extent_max.y + dydelta/2.0;
694 6 : adfGeoTransform[4] = 0.0;
695 6 : adfGeoTransform[5] = -dydelta;
696 : }
697 : /* -------------------------------------------------------------------- */
698 : /* Geographic -- use corners directly. */
699 : /* -------------------------------------------------------------------- */
700 : else
701 : {
702 18 : nRasterYSize = (int) ((extent_max.y - extent_min.y)/dydelta + 1.5);
703 18 : nRasterXSize = nProfiles;
704 :
705 : // Translate extents from arc-seconds to decimal degrees.
706 18 : adfGeoTransform[0] = (extent_min.x - dxdelta/2.0) / 3600.0;
707 18 : adfGeoTransform[1] = dxdelta / 3600.0;
708 18 : adfGeoTransform[2] = 0.0;
709 18 : adfGeoTransform[3] = (extent_max.y + dydelta/2.0) / 3600.0;
710 18 : adfGeoTransform[4] = 0.0;
711 18 : adfGeoTransform[5] = (-dydelta) / 3600.0;
712 : }
713 :
714 24 : if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
715 : {
716 0 : return FALSE;
717 : }
718 :
719 24 : return TRUE;
720 : }
721 :
722 : /************************************************************************/
723 : /* GetGeoTransform() */
724 : /************************************************************************/
725 :
726 31 : CPLErr USGSDEMDataset::GetGeoTransform( double * padfTransform )
727 :
728 : {
729 31 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
730 31 : return CE_None;
731 : }
732 :
733 : /************************************************************************/
734 : /* GetProjectionRef() */
735 : /************************************************************************/
736 :
737 47 : const char *USGSDEMDataset::GetProjectionRef()
738 :
739 : {
740 47 : return pszProjection;
741 : }
742 :
743 :
744 : /************************************************************************/
745 : /* Identify() */
746 : /************************************************************************/
747 :
748 11903 : int USGSDEMDataset::Identify( GDALOpenInfo * poOpenInfo )
749 :
750 : {
751 11903 : if( poOpenInfo->nHeaderBytes < 200 )
752 11381 : return FALSE;
753 :
754 522 : if( !EQUALN((const char *) poOpenInfo->pabyHeader+156, " 0",6)
755 : && !EQUALN((const char *) poOpenInfo->pabyHeader+156, " 1",6)
756 : && !EQUALN((const char *) poOpenInfo->pabyHeader+156, " 2",6)
757 : && !EQUALN((const char *) poOpenInfo->pabyHeader+156, " 3",6)
758 : && !EQUALN((const char *) poOpenInfo->pabyHeader+156, " -9999",6) )
759 498 : return FALSE;
760 :
761 24 : if( !EQUALN((const char *) poOpenInfo->pabyHeader+150, " 1",6)
762 : && !EQUALN((const char *) poOpenInfo->pabyHeader+150, " 4",6))
763 0 : return FALSE;
764 :
765 24 : return TRUE;
766 : }
767 :
768 : /************************************************************************/
769 : /* Open() */
770 : /************************************************************************/
771 :
772 1774 : GDALDataset *USGSDEMDataset::Open( GDALOpenInfo * poOpenInfo )
773 :
774 : {
775 1774 : if( !Identify( poOpenInfo ) )
776 1750 : return NULL;
777 :
778 24 : VSILFILE* fp = VSIFOpenL(poOpenInfo->pszFilename, "rb");
779 24 : if (fp == NULL)
780 0 : return NULL;
781 :
782 : /* -------------------------------------------------------------------- */
783 : /* Create a corresponding GDALDataset. */
784 : /* -------------------------------------------------------------------- */
785 : USGSDEMDataset *poDS;
786 :
787 24 : poDS = new USGSDEMDataset();
788 :
789 24 : poDS->fp = fp;
790 :
791 : /* -------------------------------------------------------------------- */
792 : /* Read the file. */
793 : /* -------------------------------------------------------------------- */
794 24 : if( !poDS->LoadFromFile( poDS->fp ) )
795 : {
796 0 : delete poDS;
797 0 : return NULL;
798 : }
799 :
800 : /* -------------------------------------------------------------------- */
801 : /* Confirm the requested access is supported. */
802 : /* -------------------------------------------------------------------- */
803 24 : if( poOpenInfo->eAccess == GA_Update )
804 : {
805 0 : delete poDS;
806 : CPLError( CE_Failure, CPLE_NotSupported,
807 : "The USGSDEM driver does not support update access to existing"
808 0 : " datasets.\n" );
809 0 : return NULL;
810 : }
811 :
812 : /* -------------------------------------------------------------------- */
813 : /* Create band information objects. */
814 : /* -------------------------------------------------------------------- */
815 24 : poDS->SetBand( 1, new USGSDEMRasterBand( poDS ));
816 :
817 24 : poDS->SetMetadataItem( GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT );
818 :
819 : /* -------------------------------------------------------------------- */
820 : /* Initialize any PAM information. */
821 : /* -------------------------------------------------------------------- */
822 24 : poDS->SetDescription( poOpenInfo->pszFilename );
823 24 : poDS->TryLoadXML();
824 :
825 : /* -------------------------------------------------------------------- */
826 : /* Open overviews. */
827 : /* -------------------------------------------------------------------- */
828 24 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
829 :
830 24 : return( poDS );
831 : }
832 :
833 : /************************************************************************/
834 : /* GDALRegister_USGSDEM() */
835 : /************************************************************************/
836 :
837 610 : void GDALRegister_USGSDEM()
838 :
839 : {
840 : GDALDriver *poDriver;
841 :
842 610 : if( GDALGetDriverByName( "USGSDEM" ) == NULL )
843 : {
844 588 : poDriver = new GDALDriver();
845 :
846 588 : poDriver->SetDescription( "USGSDEM" );
847 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION,
848 588 : "dem" );
849 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
850 588 : "USGS Optional ASCII DEM (and CDED)" );
851 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
852 588 : "frmt_usgsdem.html" );
853 588 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, "Int16" );
854 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
855 : "<CreationOptionList>"
856 : " <Option name='PRODUCT' type='string-select' description='Specific Product Type'>"
857 : " <Value>DEFAULT</Value>"
858 : " <Value>CDED50K</Value>"
859 : " </Option>"
860 : " <Option name='TOPLEFT' type='string' description='Top left product corner (ie. 117d15w,52d30n'/>"
861 : " <Option name='RESAMPLE' type='string-select' description='Resampling kernel to use if resampled.'>"
862 : " <Value>Nearest</Value>"
863 : " <Value>Bilinear</Value>"
864 : " <Value>Cubic</Value>"
865 : " <Value>CubicSpline</Value>"
866 : " </Option>"
867 : " <Option name='TEMPLATE' type='string' description='File to default metadata from.'/>"
868 : " <Option name='DEMLevelCode' type='int' description='DEM Level (1, 2 or 3 if set)'/>"
869 : " <Option name='DataSpecVersion' type='int' description='Data and Specification version/revision (eg. 1020)'/>"
870 : " <Option name='PRODUCER' type='string' description='Producer Agency (up to 60 characters)'/>"
871 : " <Option name='OriginCode' type='string' description='Origin code (up to 4 characters, YT for Yukon)'/>"
872 : " <Option name='ProcessCode' type='string' description='Processing Code (8=ANUDEM, 9=FME, A=TopoGrid)'/>"
873 : " <Option name='ZRESOLUTION' type='float' description='Scaling factor for elevation values'/>"
874 : " <Option name='NTS' type='string' description='NTS Mapsheet name, used to derive TOPLEFT.'/>"
875 : " <Option name='INTERNALNAME' type='string' description='Dataset name written into file header.'/>"
876 588 : "</CreationOptionList>" );
877 :
878 588 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
879 :
880 588 : poDriver->pfnOpen = USGSDEMDataset::Open;
881 588 : poDriver->pfnCreateCopy = USGSDEMCreateCopy;
882 588 : poDriver->pfnIdentify = USGSDEMDataset::Identify;
883 :
884 588 : GetGDALDriverManager()->RegisterDriver( poDriver );
885 : }
886 610 : }
|
