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