1 : /******************************************************************************
2 : * $Id: usgsdem_create.cpp 21680 2011-02-11 21:12:07Z warmerdam $
3 : *
4 : * Project: USGS DEM Driver
5 : * Purpose: CreateCopy() implementation.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : * This writing code based on the format specification:
9 : * Canadian Digital Elevation Data Product Specification - Edition 2.0
10 : *
11 : ******************************************************************************
12 : * Copyright (c) 2004, 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 : #include "cpl_string.h"
36 : #include "gdalwarper.h"
37 : #include "cpl_csv.h"
38 :
39 : CPL_CVSID("$Id: usgsdem_create.cpp 21680 2011-02-11 21:12:07Z warmerdam $");
40 :
41 : typedef struct
42 : {
43 : GDALDataset *poSrcDS;
44 : char *pszFilename;
45 : int nXSize, nYSize;
46 :
47 : char *pszDstSRS;
48 :
49 : double dfLLX, dfLLY; // These are adjusted in to center of
50 : double dfULX, dfULY; // corner pixels, and in decimal degrees.
51 : double dfURX, dfURY;
52 : double dfLRX, dfLRY;
53 :
54 : int utmzone;
55 : char horizdatum[2];
56 :
57 : double dfHorizStepSize;
58 : double dfVertStepSize;
59 : double dfElevStepSize;
60 :
61 : char **papszOptions;
62 : int bStrict;
63 :
64 : VSILFILE *fp;
65 :
66 : GInt16 *panData;
67 :
68 : } USGSDEMWriteInfo;
69 :
70 : #define DEM_NODATA -32767
71 :
72 : /************************************************************************/
73 : /* USGSDEMWriteCleanup() */
74 : /************************************************************************/
75 :
76 36 : static void USGSDEMWriteCleanup( USGSDEMWriteInfo *psWInfo )
77 :
78 : {
79 36 : CSLDestroy( psWInfo->papszOptions );
80 36 : CPLFree( psWInfo->pszDstSRS );
81 36 : CPLFree( psWInfo->pszFilename );
82 36 : if( psWInfo->fp != NULL )
83 32 : VSIFCloseL( psWInfo->fp );
84 36 : if( psWInfo->panData != NULL )
85 36 : VSIFree( psWInfo->panData );
86 36 : }
87 :
88 : /************************************************************************/
89 : /* USGSDEMDectoPackedDMS() */
90 : /************************************************************************/
91 60 : const char *USGSDEMDecToPackedDMS( double dfDec )
92 : {
93 : double dfSeconds;
94 : int nDegrees, nMinutes, nSign;
95 : static char szPackBuf[100];
96 :
97 60 : nSign = ( dfDec < 0.0 )? -1 : 1;
98 :
99 60 : dfDec = ABS( dfDec );
100 : /* If the difference between the value and the nearest degree
101 : is less than 1e-5 second, then we force to round to the
102 : nearest degree, to avoid result strings like '40 59 60.0000' instead of '41'.
103 : This is of general interest, but was mainly done to workaround a strange
104 : Valgrind bug when running usgsdem_6 where the value of psDInfo->dfULCornerY
105 : computed in DTEDOpen() differ between Valgrind and non-Valgrind executions.
106 : */
107 60 : if (fabs(dfDec - (int) floor( dfDec + .5)) < 1e-5 / 3600)
108 14 : dfDec = nDegrees = (int) floor( dfDec + .5);
109 : else
110 46 : nDegrees = (int) floor( dfDec );
111 60 : nMinutes = (int) floor( ( dfDec - nDegrees ) * 60.0 );
112 60 : dfSeconds = (dfDec - nDegrees) * 3600.0 - nMinutes * 60.0;
113 :
114 : sprintf( szPackBuf, "%4d%2d%7.4f",
115 60 : nSign * nDegrees, nMinutes, dfSeconds );
116 60 : return szPackBuf;
117 : }
118 :
119 : /************************************************************************/
120 : /* TextFill() */
121 : /************************************************************************/
122 :
123 550 : static void TextFill( char *pszTarget, unsigned int nMaxChars,
124 : const char *pszSrc )
125 :
126 : {
127 550 : if( strlen(pszSrc) < nMaxChars )
128 : {
129 456 : memcpy( pszTarget, pszSrc, strlen(pszSrc) );
130 456 : memset( pszTarget + strlen(pszSrc), ' ', nMaxChars - strlen(pszSrc));
131 : }
132 : else
133 : {
134 94 : memcpy( pszTarget, pszSrc, nMaxChars );
135 : }
136 550 : }
137 :
138 : /************************************************************************/
139 : /* TextFillR() */
140 : /* */
141 : /* Right justified. */
142 : /************************************************************************/
143 :
144 3008332 : static void TextFillR( char *pszTarget, unsigned int nMaxChars,
145 : const char *pszSrc )
146 :
147 : {
148 3008332 : if( strlen(pszSrc) < nMaxChars )
149 : {
150 2993046 : memset( pszTarget, ' ', nMaxChars - strlen(pszSrc) );
151 : memcpy( pszTarget + nMaxChars - strlen(pszSrc), pszSrc,
152 2993046 : strlen(pszSrc) );
153 : }
154 : else
155 15286 : memcpy( pszTarget, pszSrc, nMaxChars );
156 3008332 : }
157 :
158 : /************************************************************************/
159 : /* USGSDEMPrintDouble() */
160 : /* */
161 : /* The MSVC C runtime library uses 3 digits */
162 : /* for the exponent. This causes various problems, so we try */
163 : /* to correct it here. */
164 : /************************************************************************/
165 :
166 : #if defined(_MSC_VER) || defined(__MSVCRT__)
167 : # define MSVC_HACK
168 : #endif
169 :
170 13728 : static void USGSDEMPrintDouble( char *pszBuffer, double dfValue )
171 :
172 : {
173 : #define DOUBLE_BUFFER_SIZE 64
174 :
175 : char szTemp[DOUBLE_BUFFER_SIZE];
176 : int i;
177 : #ifdef MSVC_HACK
178 : const char *pszFormat = "%25.15e";
179 : #else
180 13728 : const char *pszFormat = "%24.15e";
181 : #endif
182 :
183 13728 : if ( !pszBuffer )
184 0 : return;
185 :
186 : #if defined(HAVE_SNPRINTF)
187 13728 : snprintf( szTemp, DOUBLE_BUFFER_SIZE, pszFormat, dfValue );
188 : #else
189 : sprintf( szTemp, pszFormat, dfValue );
190 : #endif
191 13728 : szTemp[DOUBLE_BUFFER_SIZE - 1] = '\0';
192 :
193 343200 : for( i = 0; szTemp[i] != '\0'; i++ )
194 : {
195 329472 : if( szTemp[i] == 'E' || szTemp[i] == 'e' )
196 13728 : szTemp[i] = 'D';
197 : #ifdef MSVC_HACK
198 : if( (szTemp[i] == '+' || szTemp[i] == '-')
199 : && szTemp[i+1] == '0' && isdigit(szTemp[i+2])
200 : && isdigit(szTemp[i+3]) && szTemp[i+4] == '\0' )
201 : {
202 : memmove( szTemp+i+1, szTemp+i+2, 2 );
203 : szTemp[i+3] = '\0';
204 : break;
205 : }
206 : #endif
207 : }
208 :
209 13728 : TextFillR( pszBuffer, 24, szTemp );
210 : }
211 :
212 : /************************************************************************/
213 : /* USGSDEMPrintSingle() */
214 : /* */
215 : /* The MSVC C runtime library uses 3 digits */
216 : /* for the exponent. This causes various problems, so we try */
217 : /* to correct it here. */
218 : /************************************************************************/
219 :
220 96 : static void USGSDEMPrintSingle( char *pszBuffer, double dfValue )
221 :
222 : {
223 : #define DOUBLE_BUFFER_SIZE 64
224 :
225 : char szTemp[DOUBLE_BUFFER_SIZE];
226 : int i;
227 : #ifdef MSVC_HACK
228 : const char *pszFormat = "%13.6e";
229 : #else
230 96 : const char *pszFormat = "%12.6e";
231 : #endif
232 :
233 96 : if ( !pszBuffer )
234 0 : return;
235 :
236 : #if defined(HAVE_SNPRINTF)
237 96 : snprintf( szTemp, DOUBLE_BUFFER_SIZE, pszFormat, dfValue );
238 : #else
239 : sprintf( szTemp, pszFormat, dfValue );
240 : #endif
241 96 : szTemp[DOUBLE_BUFFER_SIZE - 1] = '\0';
242 :
243 1248 : for( i = 0; szTemp[i] != '\0'; i++ )
244 : {
245 1152 : if( szTemp[i] == 'E' || szTemp[i] == 'e' )
246 96 : szTemp[i] = 'D';
247 : #ifdef MSVC_HACK
248 : if( (szTemp[i] == '+' || szTemp[i] == '-')
249 : && szTemp[i+1] == '0' && isdigit(szTemp[i+2])
250 : && isdigit(szTemp[i+3]) && szTemp[i+4] == '\0' )
251 : {
252 : memmove( szTemp+i+1, szTemp+i+2, 2 );
253 : szTemp[i+3] = '\0';
254 : break;
255 : }
256 : #endif
257 : }
258 :
259 96 : TextFillR( pszBuffer, 12, szTemp );
260 : }
261 :
262 : /************************************************************************/
263 : /* USGSDEMWriteARecord() */
264 : /************************************************************************/
265 :
266 32 : static int USGSDEMWriteARecord( USGSDEMWriteInfo *psWInfo )
267 :
268 : {
269 : char achARec[1024];
270 : int i;
271 : const char *pszOption;
272 :
273 : /* -------------------------------------------------------------------- */
274 : /* Init to blanks. */
275 : /* -------------------------------------------------------------------- */
276 32 : memset( achARec, ' ', sizeof(achARec) );
277 :
278 : /* -------------------------------------------------------------------- */
279 : /* Load template file, if one is indicated. */
280 : /* -------------------------------------------------------------------- */
281 : const char *pszTemplate =
282 32 : CSLFetchNameValue( psWInfo->papszOptions, "TEMPLATE" );
283 32 : if( pszTemplate != NULL )
284 : {
285 : VSILFILE *fpTemplate;
286 :
287 2 : fpTemplate = VSIFOpenL( pszTemplate, "rb" );
288 2 : if( fpTemplate == NULL )
289 : {
290 : CPLError( CE_Failure, CPLE_OpenFailed,
291 : "Unable to open template file '%s'.\n%s",
292 0 : pszTemplate, VSIStrerror( errno ) );
293 0 : return FALSE;
294 : }
295 :
296 2 : if( VSIFReadL( achARec, 1, 1024, fpTemplate ) != 1024 )
297 : {
298 : CPLError( CE_Failure, CPLE_FileIO,
299 : "Unable to read 1024 byte A Record from template file '%s'.\n%s",
300 0 : pszTemplate, VSIStrerror( errno ) );
301 0 : return FALSE;
302 : }
303 2 : VSIFCloseL( fpTemplate );
304 : }
305 :
306 : /* -------------------------------------------------------------------- */
307 : /* Filename (right justify) */
308 : /* -------------------------------------------------------------------- */
309 32 : TextFillR( achARec + 0, 40, CPLGetFilename( psWInfo->pszFilename ) );
310 :
311 : /* -------------------------------------------------------------------- */
312 : /* Producer */
313 : /* -------------------------------------------------------------------- */
314 32 : pszOption = CSLFetchNameValue( psWInfo->papszOptions, "PRODUCER" );
315 :
316 32 : if( pszOption != NULL )
317 2 : TextFillR( achARec + 40, 60, pszOption );
318 :
319 30 : else if( pszTemplate == NULL )
320 28 : TextFill( achARec + 40, 60, "" );
321 :
322 : /* -------------------------------------------------------------------- */
323 : /* Filler */
324 : /* -------------------------------------------------------------------- */
325 32 : TextFill( achARec + 100, 9, "" );
326 :
327 : /* -------------------------------------------------------------------- */
328 : /* SW Geographic Corner - SDDDMMSS.SSSS - longitude then latitude */
329 : /* -------------------------------------------------------------------- */
330 32 : if ( ! psWInfo->utmzone )
331 : {
332 : TextFill( achARec + 109, 13,
333 30 : USGSDEMDecToPackedDMS( psWInfo->dfLLX ) ); // longitude
334 : TextFill( achARec + 122, 13,
335 30 : USGSDEMDecToPackedDMS( psWInfo->dfLLY ) ); // latitude
336 : }
337 : /* this may not be best according to the spec. But for now,
338 : * we won't try to convert the UTM coordinates to lat/lon
339 : */
340 :
341 : /* -------------------------------------------------------------------- */
342 : /* Process code. */
343 : /* -------------------------------------------------------------------- */
344 32 : pszOption = CSLFetchNameValue( psWInfo->papszOptions, "ProcessCode" );
345 :
346 32 : if( pszOption != NULL )
347 2 : TextFill( achARec + 135, 1, pszOption );
348 :
349 30 : else if( pszTemplate == NULL )
350 28 : TextFill( achARec + 135, 1, " " );
351 :
352 : /* -------------------------------------------------------------------- */
353 : /* Filler */
354 : /* -------------------------------------------------------------------- */
355 32 : TextFill( achARec + 136, 1, "" );
356 :
357 : /* -------------------------------------------------------------------- */
358 : /* Sectional indicator */
359 : /* -------------------------------------------------------------------- */
360 32 : if( pszTemplate == NULL )
361 30 : TextFill( achARec + 137, 3, "" );
362 :
363 : /* -------------------------------------------------------------------- */
364 : /* Origin code */
365 : /* -------------------------------------------------------------------- */
366 32 : pszOption = CSLFetchNameValue( psWInfo->papszOptions, "OriginCode" );
367 :
368 32 : if( pszOption != NULL )
369 2 : TextFill( achARec + 140, 4, pszOption ); // Should be YT for Yukon.
370 :
371 30 : else if( pszTemplate == NULL )
372 28 : TextFill( achARec + 140, 4, "" );
373 :
374 : /* -------------------------------------------------------------------- */
375 : /* DEM level code (right justify) */
376 : /* -------------------------------------------------------------------- */
377 32 : pszOption = CSLFetchNameValue( psWInfo->papszOptions, "DEMLevelCode" );
378 :
379 32 : if( pszOption != NULL )
380 2 : TextFillR( achARec + 144, 6, pszOption ); // 1, 2 or 3.
381 :
382 30 : else if( pszTemplate == NULL )
383 28 : TextFillR( achARec + 144, 6, "1" ); // 1, 2 or 3.
384 : /* some DEM readers require a value, 1 seems to be a
385 : * default
386 : */
387 :
388 : /* -------------------------------------------------------------------- */
389 : /* Elevation Pattern */
390 : /* -------------------------------------------------------------------- */
391 32 : TextFillR( achARec + 150, 6, "1" ); // "1" for regular (random is 2)
392 :
393 : /* -------------------------------------------------------------------- */
394 : /* Horizontal Reference System. */
395 : /* */
396 : /* 0 = Geographic */
397 : /* 1 = UTM */
398 : /* 2 = Stateplane */
399 : /* -------------------------------------------------------------------- */
400 32 : if ( ! psWInfo->utmzone )
401 : {
402 30 : TextFillR( achARec + 156, 6, "0" );
403 : }
404 : else
405 : {
406 2 : TextFillR( achARec + 156, 6, "1" );
407 : }
408 :
409 : /* -------------------------------------------------------------------- */
410 : /* UTM / State Plane zone. */
411 : /* -------------------------------------------------------------------- */
412 32 : if ( ! psWInfo->utmzone )
413 : {
414 30 : TextFillR( achARec + 162, 6, "0");
415 : }
416 : else
417 : {
418 : TextFillR( achARec + 162, 6,
419 2 : CPLSPrintf( "%02d", psWInfo->utmzone) );
420 : }
421 :
422 : /* -------------------------------------------------------------------- */
423 : /* Map Projection Parameters (all 0.0). */
424 : /* -------------------------------------------------------------------- */
425 512 : for( i = 0; i < 15; i++ )
426 480 : TextFillR( achARec + 168 + i*24, 24, "0.0" );
427 :
428 : /* -------------------------------------------------------------------- */
429 : /* Horizontal Unit of Measure */
430 : /* 0 = radians */
431 : /* 1 = feet */
432 : /* 2 = meters */
433 : /* 3 = arc seconds */
434 : /* -------------------------------------------------------------------- */
435 32 : if ( ! psWInfo->utmzone )
436 : {
437 30 : TextFillR( achARec + 528, 6, "3" );
438 : }
439 : else
440 : {
441 2 : TextFillR( achARec + 528, 6, "2" );
442 : }
443 :
444 : /* -------------------------------------------------------------------- */
445 : /* Vertical unit of measure. */
446 : /* 1 = feet */
447 : /* 2 = meters */
448 : /* -------------------------------------------------------------------- */
449 32 : TextFillR( achARec + 534, 6, "2" );
450 :
451 : /* -------------------------------------------------------------------- */
452 : /* Number of sides in coverage polygon (always 4) */
453 : /* -------------------------------------------------------------------- */
454 32 : TextFillR( achARec + 540, 6, "4" );
455 :
456 : /* -------------------------------------------------------------------- */
457 : /* 4 corner coordinates: SW, NW, NE, SE */
458 : /* Corners are in 24.15 format in arc seconds. */
459 : /* -------------------------------------------------------------------- */
460 32 : if ( ! psWInfo->utmzone )
461 : {
462 : // SW - longitude
463 30 : USGSDEMPrintDouble( achARec + 546, psWInfo->dfLLX * 3600.0 );
464 : // SW - latitude
465 30 : USGSDEMPrintDouble( achARec + 570, psWInfo->dfLLY * 3600.0 );
466 :
467 : // NW - longitude
468 30 : USGSDEMPrintDouble( achARec + 594, psWInfo->dfULX * 3600.0 );
469 : // NW - latitude
470 30 : USGSDEMPrintDouble( achARec + 618, psWInfo->dfULY * 3600.0 );
471 :
472 : // NE - longitude
473 30 : USGSDEMPrintDouble( achARec + 642, psWInfo->dfURX * 3600.0 );
474 : // NE - latitude
475 30 : USGSDEMPrintDouble( achARec + 666, psWInfo->dfURY * 3600.0 );
476 :
477 : // SE - longitude
478 30 : USGSDEMPrintDouble( achARec + 690, psWInfo->dfLRX * 3600.0 );
479 : // SE - latitude
480 30 : USGSDEMPrintDouble( achARec + 714, psWInfo->dfLRY * 3600.0 );
481 : }
482 : else
483 : {
484 : // SW - easting
485 2 : USGSDEMPrintDouble( achARec + 546, psWInfo->dfLLX );
486 : // SW - northing
487 2 : USGSDEMPrintDouble( achARec + 570, psWInfo->dfLLY );
488 :
489 : // NW - easting
490 2 : USGSDEMPrintDouble( achARec + 594, psWInfo->dfULX );
491 : // NW - northing
492 2 : USGSDEMPrintDouble( achARec + 618, psWInfo->dfULY );
493 :
494 : // NE - easting
495 2 : USGSDEMPrintDouble( achARec + 642, psWInfo->dfURX );
496 : // NE - northing
497 2 : USGSDEMPrintDouble( achARec + 666, psWInfo->dfURY );
498 :
499 : // SE - easting
500 2 : USGSDEMPrintDouble( achARec + 690, psWInfo->dfLRX );
501 : // SE - northing
502 2 : USGSDEMPrintDouble( achARec + 714, psWInfo->dfLRY );
503 : }
504 :
505 : /* -------------------------------------------------------------------- */
506 : /* Minimum and Maximum elevations for this cell. */
507 : /* 24.15 format. */
508 : /* -------------------------------------------------------------------- */
509 32 : GInt16 nMin = DEM_NODATA, nMax = DEM_NODATA;
510 32 : int nVoid = 0;
511 :
512 2976760 : for( i = psWInfo->nXSize*psWInfo->nYSize-1; i >= 0; i-- )
513 : {
514 2976728 : if( psWInfo->panData[i] != DEM_NODATA )
515 : {
516 2975298 : if( nMin == DEM_NODATA )
517 : {
518 32 : nMin = nMax = psWInfo->panData[i];
519 : }
520 : else
521 : {
522 2975266 : nMin = MIN(nMin,psWInfo->panData[i]);
523 2975266 : nMax = MAX(nMax,psWInfo->panData[i]);
524 : }
525 : }
526 : else
527 1430 : nVoid++;
528 : }
529 :
530 : /* take into account z resolutions that are not 1.0 */
531 32 : nMin = (GInt16) floor(nMin * psWInfo->dfElevStepSize);
532 32 : nMax = (GInt16) ceil(nMax * psWInfo->dfElevStepSize);
533 :
534 32 : USGSDEMPrintDouble( achARec + 738, (double) nMin );
535 32 : USGSDEMPrintDouble( achARec + 762, (double) nMax );
536 :
537 : /* -------------------------------------------------------------------- */
538 : /* Counter Clockwise angle (in radians). Normally 0 */
539 : /* -------------------------------------------------------------------- */
540 32 : TextFillR( achARec + 786, 24, "0.0" );
541 :
542 : /* -------------------------------------------------------------------- */
543 : /* Accurancy code for elevations. 0 means there will be no C */
544 : /* record. */
545 : /* -------------------------------------------------------------------- */
546 32 : TextFillR( achARec + 810, 6, "0" );
547 :
548 : /* -------------------------------------------------------------------- */
549 : /* Spatial Resolution (x, y and z). 12.6 format. */
550 : /* -------------------------------------------------------------------- */
551 32 : if ( ! psWInfo->utmzone )
552 : {
553 : USGSDEMPrintSingle( achARec + 816,
554 30 : psWInfo->dfHorizStepSize*3600.0 );
555 : USGSDEMPrintSingle( achARec + 828,
556 30 : psWInfo->dfVertStepSize*3600.0 );
557 : }
558 : else
559 : {
560 : USGSDEMPrintSingle( achARec + 816,
561 2 : psWInfo->dfHorizStepSize );
562 : USGSDEMPrintSingle( achARec + 828,
563 2 : psWInfo->dfVertStepSize );
564 : }
565 :
566 32 : USGSDEMPrintSingle( achARec + 840, psWInfo->dfElevStepSize);
567 :
568 : /* -------------------------------------------------------------------- */
569 : /* Rows and Columns of profiles. */
570 : /* -------------------------------------------------------------------- */
571 32 : TextFillR( achARec + 852, 6, CPLSPrintf( "%d", 1 ) );
572 32 : TextFillR( achARec + 858, 6, CPLSPrintf( "%d", psWInfo->nXSize ) );
573 :
574 : /* -------------------------------------------------------------------- */
575 : /* Largest primary contour interval (blank). */
576 : /* -------------------------------------------------------------------- */
577 32 : TextFill( achARec + 864, 5, "" );
578 :
579 : /* -------------------------------------------------------------------- */
580 : /* Largest source contour internal unit (blank). */
581 : /* -------------------------------------------------------------------- */
582 32 : TextFill( achARec + 869, 1, "" );
583 :
584 : /* -------------------------------------------------------------------- */
585 : /* Smallest primary contour interval. */
586 : /* -------------------------------------------------------------------- */
587 32 : TextFill( achARec + 870, 5, "" );
588 :
589 : /* -------------------------------------------------------------------- */
590 : /* Smallest source contour interval unit. */
591 : /* -------------------------------------------------------------------- */
592 32 : TextFill( achARec + 875, 1, "" );
593 :
594 : /* -------------------------------------------------------------------- */
595 : /* Data source data - YYMM */
596 : /* -------------------------------------------------------------------- */
597 32 : if( pszTemplate == NULL )
598 30 : TextFill( achARec + 876, 4, "" );
599 :
600 : /* -------------------------------------------------------------------- */
601 : /* Data inspection/revision data (YYMM). */
602 : /* -------------------------------------------------------------------- */
603 32 : if( pszTemplate == NULL )
604 30 : TextFill( achARec + 880, 4, "" );
605 :
606 : /* -------------------------------------------------------------------- */
607 : /* Inspection revision flag (I or R) (blank) */
608 : /* -------------------------------------------------------------------- */
609 32 : if( pszTemplate == NULL )
610 30 : TextFill( achARec + 884, 1, "" );
611 :
612 : /* -------------------------------------------------------------------- */
613 : /* Data validation flag. */
614 : /* -------------------------------------------------------------------- */
615 32 : if( pszTemplate == NULL )
616 30 : TextFill( achARec + 885, 1, "" );
617 :
618 : /* -------------------------------------------------------------------- */
619 : /* Suspect and void area flag. */
620 : /* 0 = none */
621 : /* 1 = suspect areas */
622 : /* 2 = void areas */
623 : /* 3 = suspect and void areas */
624 : /* -------------------------------------------------------------------- */
625 32 : if( nVoid > 0 )
626 2 : TextFillR( achARec + 886, 2, "2" );
627 : else
628 30 : TextFillR( achARec + 886, 2, "0" );
629 :
630 : /* -------------------------------------------------------------------- */
631 : /* Vertical datum */
632 : /* 1 = MSL */
633 : /* 2 = NGVD29 */
634 : /* 3 = NAVD88 */
635 : /* -------------------------------------------------------------------- */
636 32 : if( pszTemplate == NULL )
637 30 : TextFillR( achARec + 888, 2, "1" );
638 :
639 : /* -------------------------------------------------------------------- */
640 : /* Horizonal Datum */
641 : /* 1 = NAD27 */
642 : /* 2 = WGS72 */
643 : /* 3 = WGS84 */
644 : /* 4 = NAD83 */
645 : /* -------------------------------------------------------------------- */
646 32 : if( strlen( psWInfo->horizdatum ) == 0) {
647 0 : if( pszTemplate == NULL )
648 0 : TextFillR( achARec + 890, 2, "4" );
649 : }
650 : else
651 : {
652 32 : if( pszTemplate == NULL )
653 30 : TextFillR( achARec + 890, 2, psWInfo->horizdatum );
654 : }
655 :
656 : /* -------------------------------------------------------------------- */
657 : /* Data edition/version, specification edition/version. */
658 : /* -------------------------------------------------------------------- */
659 32 : pszOption = CSLFetchNameValue( psWInfo->papszOptions, "DataSpecVersion" );
660 :
661 32 : if( pszOption != NULL )
662 2 : TextFill( achARec + 892, 4, pszOption );
663 :
664 30 : else if( pszTemplate == NULL )
665 28 : TextFill( achARec + 892, 4, "" );
666 :
667 : /* -------------------------------------------------------------------- */
668 : /* Percent void. */
669 : /* */
670 : /* Round to nearest integer percentage. */
671 : /* -------------------------------------------------------------------- */
672 : int nPercent;
673 :
674 : nPercent = (int)
675 32 : (((nVoid * 100.0) / (psWInfo->nXSize * psWInfo->nYSize)) + 0.5);
676 :
677 32 : TextFillR( achARec + 896, 4, CPLSPrintf( "%4d", nPercent ) );
678 :
679 : /* -------------------------------------------------------------------- */
680 : /* Edge matching flags. */
681 : /* -------------------------------------------------------------------- */
682 32 : if( pszTemplate == NULL )
683 30 : TextFill( achARec + 900, 8, "" );
684 :
685 : /* -------------------------------------------------------------------- */
686 : /* Vertical datum shift (F7.2). */
687 : /* -------------------------------------------------------------------- */
688 32 : TextFillR( achARec + 908, 7, "" );
689 :
690 : /* -------------------------------------------------------------------- */
691 : /* Write to file. */
692 : /* -------------------------------------------------------------------- */
693 32 : if( VSIFWriteL( achARec, 1, 1024, psWInfo->fp ) != 1024 )
694 : {
695 : CPLError( CE_Failure, CPLE_FileIO,
696 : "Error writing DEM/CDED A record.\n%s",
697 0 : VSIStrerror( errno ) );
698 0 : return FALSE;
699 : }
700 :
701 32 : return TRUE;
702 : }
703 :
704 : /************************************************************************/
705 : /* USGSDEMWriteProfile() */
706 : /* */
707 : /* Write B logical record. Split into 1024 byte chunks. */
708 : /************************************************************************/
709 :
710 3352 : static int USGSDEMWriteProfile( USGSDEMWriteInfo *psWInfo, int iProfile )
711 :
712 : {
713 : char achBuffer[1024];
714 :
715 3352 : memset( achBuffer, ' ', sizeof(achBuffer) );
716 :
717 : /* -------------------------------------------------------------------- */
718 : /* Row #. */
719 : /* -------------------------------------------------------------------- */
720 3352 : TextFillR( achBuffer + 0, 6, "1" );
721 :
722 : /* -------------------------------------------------------------------- */
723 : /* Column #. */
724 : /* -------------------------------------------------------------------- */
725 3352 : TextFillR( achBuffer + 6, 6, CPLSPrintf( "%d", iProfile + 1 ) );
726 :
727 : /* -------------------------------------------------------------------- */
728 : /* Number of data items. */
729 : /* -------------------------------------------------------------------- */
730 3352 : TextFillR( achBuffer + 12, 6, CPLSPrintf( "%d", psWInfo->nYSize ) );
731 3352 : TextFillR( achBuffer + 18, 6, "1" );
732 :
733 : /* -------------------------------------------------------------------- */
734 : /* Location of center of bottom most sample in profile. */
735 : /* Format D24.15. In arc-seconds if geographic, meters */
736 : /* if UTM. */
737 : /* -------------------------------------------------------------------- */
738 3352 : if( ! psWInfo->utmzone )
739 : {
740 : // longitude
741 : USGSDEMPrintDouble( achBuffer + 24,
742 : 3600 * (psWInfo->dfLLX
743 3348 : + iProfile * psWInfo->dfHorizStepSize) );
744 :
745 : // latitude
746 3348 : USGSDEMPrintDouble( achBuffer + 48, psWInfo->dfLLY * 3600.0 );
747 : }
748 : else
749 : {
750 : // easting
751 : USGSDEMPrintDouble( achBuffer + 24,
752 : (psWInfo->dfLLX
753 4 : + iProfile * psWInfo->dfHorizStepSize) );
754 :
755 : // northing
756 4 : USGSDEMPrintDouble( achBuffer + 48, psWInfo->dfLLY );
757 : }
758 :
759 :
760 : /* -------------------------------------------------------------------- */
761 : /* Local vertical datum offset. */
762 : /* -------------------------------------------------------------------- */
763 3352 : TextFillR( achBuffer + 72, 24, "0.000000D+00" );
764 :
765 : /* -------------------------------------------------------------------- */
766 : /* Min/Max elevation values for this profile. */
767 : /* -------------------------------------------------------------------- */
768 : int iY;
769 3352 : GInt16 nMin = DEM_NODATA, nMax = DEM_NODATA;
770 :
771 2980080 : for( iY = 0; iY < psWInfo->nYSize; iY++ )
772 : {
773 2976728 : int iData = (psWInfo->nYSize-iY-1) * psWInfo->nXSize + iProfile;
774 :
775 2976728 : if( psWInfo->panData[iData] != DEM_NODATA )
776 : {
777 2975298 : if( nMin == DEM_NODATA )
778 : {
779 3352 : nMin = nMax = psWInfo->panData[iData];
780 : }
781 : else
782 : {
783 2971946 : nMin = MIN(nMin,psWInfo->panData[iData]);
784 2971946 : nMax = MAX(nMax,psWInfo->panData[iData]);
785 : }
786 : }
787 : }
788 :
789 : /* take into account z resolutions that are not 1.0 */
790 3352 : nMin = (GInt16) floor(nMin * psWInfo->dfElevStepSize);
791 3352 : nMax = (GInt16) ceil(nMax * psWInfo->dfElevStepSize);
792 :
793 3352 : USGSDEMPrintDouble( achBuffer + 96, (double) nMin );
794 3352 : USGSDEMPrintDouble( achBuffer + 120, (double) nMax );
795 :
796 : /* -------------------------------------------------------------------- */
797 : /* Output all the actually elevation values, flushing blocks */
798 : /* when they fill up. */
799 : /* -------------------------------------------------------------------- */
800 3352 : int iOffset = 144;
801 :
802 2980080 : for( iY = 0; iY < psWInfo->nYSize; iY++ )
803 : {
804 2976728 : int iData = (psWInfo->nYSize-iY-1) * psWInfo->nXSize + iProfile;
805 : char szWord[10];
806 :
807 2976728 : if( iOffset + 6 > 1024 )
808 : {
809 16822 : if( VSIFWriteL( achBuffer, 1, 1024, psWInfo->fp ) != 1024 )
810 : {
811 : CPLError( CE_Failure, CPLE_FileIO,
812 : "Failure writing profile to disk.\n%s",
813 0 : VSIStrerror( errno ) );
814 0 : return FALSE;
815 : }
816 16822 : iOffset = 0;
817 16822 : memset( achBuffer, ' ', 1024 );
818 : }
819 :
820 2976728 : sprintf( szWord, "%d", psWInfo->panData[iData] );
821 2976728 : TextFillR( achBuffer + iOffset, 6, szWord );
822 :
823 2976728 : iOffset += 6;
824 : }
825 :
826 : /* -------------------------------------------------------------------- */
827 : /* Flush final partial block. */
828 : /* -------------------------------------------------------------------- */
829 3352 : if( VSIFWriteL( achBuffer, 1, 1024, psWInfo->fp ) != 1024 )
830 : {
831 : CPLError( CE_Failure, CPLE_FileIO,
832 : "Failure writing profile to disk.\n%s",
833 0 : VSIStrerror( errno ) );
834 0 : return FALSE;
835 : }
836 :
837 3352 : return TRUE;
838 : }
839 :
840 : /************************************************************************/
841 : /* USGSDEM_LookupNTSByLoc() */
842 : /************************************************************************/
843 :
844 : static int
845 4 : USGSDEM_LookupNTSByLoc( double dfULLong, double dfULLat,
846 : char *pszTile, char *pszName )
847 :
848 : {
849 : /* -------------------------------------------------------------------- */
850 : /* Access NTS 1:50k sheet CSV file. */
851 : /* -------------------------------------------------------------------- */
852 4 : const char *pszNTSFilename = CSVFilename( "NTS-50kindex.csv" );
853 : FILE *fpNTS;
854 :
855 4 : fpNTS = VSIFOpen( pszNTSFilename, "rb" );
856 4 : if( fpNTS == NULL )
857 : {
858 : CPLError( CE_Failure, CPLE_FileIO, "Unable to find NTS mapsheet lookup file: %s",
859 4 : pszNTSFilename );
860 4 : return FALSE;
861 : }
862 :
863 : /* -------------------------------------------------------------------- */
864 : /* Skip column titles line. */
865 : /* -------------------------------------------------------------------- */
866 0 : CSLDestroy( CSVReadParseLine( fpNTS ) );
867 :
868 : /* -------------------------------------------------------------------- */
869 : /* Find desired sheet. */
870 : /* -------------------------------------------------------------------- */
871 0 : int bGotHit = FALSE;
872 : char **papszTokens;
873 :
874 0 : while( !bGotHit
875 : && (papszTokens = CSVReadParseLine( fpNTS )) != NULL )
876 : {
877 0 : if( CSLCount( papszTokens ) != 4 )
878 0 : continue;
879 :
880 0 : if( ABS(dfULLong - atof(papszTokens[2])) < 0.01
881 0 : && ABS(dfULLat - atof(papszTokens[3])) < 0.01 )
882 : {
883 0 : bGotHit = TRUE;
884 0 : strncpy( pszTile, papszTokens[0], 7 );
885 0 : if( pszName != NULL )
886 0 : strncpy( pszName, papszTokens[1], 100 );
887 : }
888 :
889 0 : CSLDestroy( papszTokens );
890 : }
891 :
892 0 : VSIFClose( fpNTS );
893 :
894 0 : return bGotHit;
895 : }
896 :
897 : /************************************************************************/
898 : /* USGSDEM_LookupNTSByTile() */
899 : /************************************************************************/
900 :
901 : static int
902 0 : USGSDEM_LookupNTSByTile( const char *pszTile, char *pszName,
903 : double *pdfULLong, double *pdfULLat )
904 :
905 : {
906 : /* -------------------------------------------------------------------- */
907 : /* Access NTS 1:50k sheet CSV file. */
908 : /* -------------------------------------------------------------------- */
909 0 : const char *pszNTSFilename = CSVFilename( "NTS-50kindex.csv" );
910 : FILE *fpNTS;
911 :
912 0 : fpNTS = VSIFOpen( pszNTSFilename, "rb" );
913 0 : if( fpNTS == NULL )
914 : {
915 : CPLError( CE_Failure, CPLE_FileIO, "Unable to find NTS mapsheet lookup file: %s",
916 0 : pszNTSFilename );
917 0 : return FALSE;
918 : }
919 :
920 : /* -------------------------------------------------------------------- */
921 : /* Skip column titles line. */
922 : /* -------------------------------------------------------------------- */
923 0 : CSLDestroy( CSVReadParseLine( fpNTS ) );
924 :
925 : /* -------------------------------------------------------------------- */
926 : /* Find desired sheet. */
927 : /* -------------------------------------------------------------------- */
928 0 : int bGotHit = FALSE;
929 : char **papszTokens;
930 :
931 0 : while( !bGotHit
932 : && (papszTokens = CSVReadParseLine( fpNTS )) != NULL )
933 : {
934 0 : if( CSLCount( papszTokens ) != 4 )
935 0 : continue;
936 :
937 0 : if( EQUAL(pszTile,papszTokens[0]) )
938 : {
939 0 : bGotHit = TRUE;
940 0 : if( pszName != NULL )
941 0 : strncpy( pszName, papszTokens[1], 100 );
942 0 : *pdfULLong = atof(papszTokens[2]);
943 0 : *pdfULLat = atof(papszTokens[3]);
944 : }
945 :
946 0 : CSLDestroy( papszTokens );
947 : }
948 :
949 0 : VSIFClose( fpNTS );
950 :
951 0 : return bGotHit;
952 : }
953 :
954 : /************************************************************************/
955 : /* USGSDEMProductSetup_CDED50K() */
956 : /************************************************************************/
957 :
958 2 : static int USGSDEMProductSetup_CDED50K( USGSDEMWriteInfo *psWInfo )
959 :
960 : {
961 : /* -------------------------------------------------------------------- */
962 : /* Fetch TOPLEFT location so we know what cell we are dealing */
963 : /* with. */
964 : /* -------------------------------------------------------------------- */
965 : const char *pszNTS =
966 2 : CSLFetchNameValue( psWInfo->papszOptions, "NTS" );
967 : const char *pszTOPLEFT = CSLFetchNameValue( psWInfo->papszOptions,
968 2 : "TOPLEFT" );
969 2 : double dfULX = (psWInfo->dfULX+psWInfo->dfURX)*0.5;
970 2 : double dfULY = (psWInfo->dfULY+psWInfo->dfURY)*0.5;
971 :
972 : // Have we been given an explicit NTS mapsheet name?
973 2 : if( pszNTS != NULL )
974 : {
975 : char szTrimmedTile[7];
976 :
977 0 : strncpy( szTrimmedTile, pszNTS, 6 );
978 0 : szTrimmedTile[6] = '\0';
979 :
980 0 : if( !USGSDEM_LookupNTSByTile( szTrimmedTile, NULL, &dfULX, &dfULY ) )
981 0 : return FALSE;
982 :
983 0 : if( EQUALN(pszNTS+6,"e",1) )
984 0 : dfULX += (( dfULY < 68.1 ) ? 0.25 : ( dfULY < 80.1 ) ? 0.5 : 1);
985 : }
986 :
987 : // Try looking up TOPLEFT as a NTS mapsheet name.
988 2 : else if( pszTOPLEFT != NULL && strstr(pszTOPLEFT,",") == NULL
989 : && (strlen(pszTOPLEFT) == 6 || strlen(pszTOPLEFT) == 7) )
990 : {
991 : char szTrimmedTile[7];
992 :
993 0 : strncpy( szTrimmedTile, pszTOPLEFT, 6 );
994 0 : szTrimmedTile[6] = '\0';
995 :
996 0 : if( !USGSDEM_LookupNTSByTile( szTrimmedTile, NULL, &dfULX, &dfULY ) )
997 0 : return FALSE;
998 :
999 0 : if( EQUAL(pszTOPLEFT+6,"e") )
1000 0 : dfULX += (( dfULY < 68.1 ) ? 0.25 : ( dfULY < 80.1 ) ? 0.5 : 1);
1001 : }
1002 :
1003 : // Assume TOPLEFT is a long/lat corner.
1004 2 : else if( pszTOPLEFT != NULL )
1005 : {
1006 2 : char **papszTokens = CSLTokenizeString2( pszTOPLEFT, ",", 0 );
1007 :
1008 2 : if( CSLCount( papszTokens ) != 2 )
1009 : {
1010 0 : CSLDestroy( papszTokens );
1011 : CPLError( CE_Failure, CPLE_AppDefined,
1012 0 : "Failed to parse TOPLEFT, should have form like '138d15W,59d0N'." );
1013 0 : return FALSE;
1014 : }
1015 :
1016 2 : dfULX = CPLDMSToDec( papszTokens[0] );
1017 2 : dfULY = CPLDMSToDec( papszTokens[1] );
1018 2 : CSLDestroy( papszTokens );
1019 :
1020 2 : if( ABS(dfULX*4-floor(dfULX*4+0.00005)) > 0.0001
1021 : || ABS(dfULY*4-floor(dfULY*4+0.00005)) > 0.0001 )
1022 : {
1023 : CPLError( CE_Failure, CPLE_AppDefined,
1024 0 : "TOPLEFT must be on a 15\" boundary for CDED50K, but is not." );
1025 0 : return FALSE;
1026 : }
1027 : }
1028 0 : else if( strlen(psWInfo->pszFilename) == 12
1029 0 : && psWInfo->pszFilename[6] == '_'
1030 : && EQUAL(psWInfo->pszFilename+8,".dem") )
1031 : {
1032 : char szTrimmedTile[7];
1033 :
1034 0 : strncpy( szTrimmedTile, psWInfo->pszFilename, 6 );
1035 0 : szTrimmedTile[6] = '\0';
1036 :
1037 0 : if( !USGSDEM_LookupNTSByTile( szTrimmedTile, NULL, &dfULX, &dfULY ) )
1038 0 : return FALSE;
1039 :
1040 0 : if( EQUALN(psWInfo->pszFilename+7,"e",1) )
1041 0 : dfULX += (( dfULY < 68.1 ) ? 0.25 : ( dfULY < 80.1 ) ? 0.5 : 1);
1042 : }
1043 :
1044 0 : else if( strlen(psWInfo->pszFilename) == 14
1045 : && EQUALN(psWInfo->pszFilename+6,"DEM",3)
1046 : && EQUAL(psWInfo->pszFilename+10,".dem") )
1047 : {
1048 : char szTrimmedTile[7];
1049 :
1050 0 : strncpy( szTrimmedTile, psWInfo->pszFilename, 6 );
1051 0 : szTrimmedTile[6] = '\0';
1052 :
1053 0 : if( !USGSDEM_LookupNTSByTile( szTrimmedTile, NULL, &dfULX, &dfULY ) )
1054 0 : return FALSE;
1055 :
1056 0 : if( EQUALN(psWInfo->pszFilename+9,"e",1) )
1057 0 : dfULX += (( dfULY < 68.1 ) ? 0.25 : ( dfULY < 80.1 ) ? 0.5 : 1);
1058 : }
1059 :
1060 : /* -------------------------------------------------------------------- */
1061 : /* Set resolution and size information. */
1062 : /* -------------------------------------------------------------------- */
1063 :
1064 2 : dfULX = floor( dfULX * 4 + 0.00005 ) / 4.0;
1065 2 : dfULY = floor( dfULY * 4 + 0.00005 ) / 4.0;
1066 :
1067 2 : psWInfo->nXSize = 1201;
1068 2 : psWInfo->nYSize = 1201;
1069 2 : psWInfo->dfVertStepSize = 0.75 / 3600.0;
1070 :
1071 : /* Region A */
1072 2 : if( dfULY < 68.1 )
1073 : {
1074 2 : psWInfo->dfHorizStepSize = 0.75 / 3600.0;
1075 : }
1076 :
1077 : /* Region B */
1078 0 : else if( dfULY < 80.1 )
1079 : {
1080 0 : psWInfo->dfHorizStepSize = 1.5 / 3600.0;
1081 0 : dfULX = floor( dfULX * 2 + 0.001 ) / 2.0;
1082 : }
1083 :
1084 : /* Region C */
1085 : else
1086 : {
1087 0 : psWInfo->dfHorizStepSize = 3.0 / 3600.0;
1088 0 : dfULX = floor( dfULX + 0.001 );
1089 : }
1090 :
1091 : /* -------------------------------------------------------------------- */
1092 : /* Set bounds based on this top left anchor. */
1093 : /* -------------------------------------------------------------------- */
1094 :
1095 2 : psWInfo->dfULX = dfULX;
1096 2 : psWInfo->dfULY = dfULY;
1097 2 : psWInfo->dfLLX = dfULX;
1098 2 : psWInfo->dfLLY = dfULY - 0.25;
1099 2 : psWInfo->dfURX = dfULX + psWInfo->dfHorizStepSize * 1200.0;
1100 2 : psWInfo->dfURY = dfULY;
1101 2 : psWInfo->dfLRX = dfULX + psWInfo->dfHorizStepSize * 1200.0;
1102 2 : psWInfo->dfLRY = dfULY - 0.25;
1103 :
1104 : /* -------------------------------------------------------------------- */
1105 : /* Can we find the NTS 50k tile name that corresponds with */
1106 : /* this? */
1107 : /* -------------------------------------------------------------------- */
1108 : const char *pszINTERNAL =
1109 2 : CSLFetchNameValue( psWInfo->papszOptions, "INTERNALNAME" );
1110 : char szTile[10];
1111 2 : char chEWFlag = ' ';
1112 :
1113 2 : if( USGSDEM_LookupNTSByLoc( dfULX, dfULY, szTile, NULL ) )
1114 : {
1115 0 : chEWFlag = 'w';
1116 : }
1117 2 : else if( USGSDEM_LookupNTSByLoc( dfULX-0.25, dfULY, szTile, NULL ) )
1118 : {
1119 0 : chEWFlag = 'e';
1120 : }
1121 :
1122 2 : if( pszINTERNAL != NULL )
1123 : {
1124 2 : CPLFree( psWInfo->pszFilename );
1125 2 : psWInfo->pszFilename = CPLStrdup( pszINTERNAL );
1126 : }
1127 0 : else if( chEWFlag != ' ' )
1128 : {
1129 0 : CPLFree( psWInfo->pszFilename );
1130 : psWInfo->pszFilename =
1131 0 : CPLStrdup( CPLSPrintf("%sDEM%c", szTile, chEWFlag ) );
1132 : }
1133 : else
1134 : {
1135 0 : const char *pszBasename = CPLGetFilename( psWInfo->pszFilename);
1136 0 : if( !EQUALN(pszBasename+6,"DEM",3)
1137 : || strlen(pszBasename) != 10 )
1138 : CPLError( CE_Warning, CPLE_AppDefined,
1139 : "Internal filename required to be of 'nnnannDEMz', the output\n"
1140 : "filename is not of the required format, and the tile could not be\n"
1141 : "identified in the NTS mapsheet list (or the NTS mapsheet could not\n"
1142 0 : "be found). Correct output filename for correct CDED production." );
1143 : }
1144 :
1145 : /* -------------------------------------------------------------------- */
1146 : /* Set some specific options for CDED 50K. */
1147 : /* -------------------------------------------------------------------- */
1148 : psWInfo->papszOptions =
1149 2 : CSLSetNameValue( psWInfo->papszOptions, "DEMLevelCode", "1" );
1150 :
1151 2 : if( CSLFetchNameValue( psWInfo->papszOptions, "DataSpecVersion" ) == NULL )
1152 : psWInfo->papszOptions =
1153 : CSLSetNameValue( psWInfo->papszOptions, "DataSpecVersion",
1154 2 : "1020" );
1155 :
1156 : /* -------------------------------------------------------------------- */
1157 : /* Set the destination coordinate system. */
1158 : /* -------------------------------------------------------------------- */
1159 2 : OGRSpatialReference oSRS;
1160 2 : oSRS.SetWellKnownGeogCS( "NAD83" );
1161 2 : strncpy( psWInfo->horizdatum, "4", 2 ); //USGS DEM code for NAD83
1162 :
1163 2 : oSRS.exportToWkt( &(psWInfo->pszDstSRS) );
1164 :
1165 : /* -------------------------------------------------------------------- */
1166 : /* Cleanup. */
1167 : /* -------------------------------------------------------------------- */
1168 2 : CPLReadLine( NULL );
1169 :
1170 2 : return TRUE;
1171 : }
1172 :
1173 : /************************************************************************/
1174 : /* USGSDEMProductSetup_DEFAULT() */
1175 : /* */
1176 : /* Sets up the new DEM dataset parameters, using the source */
1177 : /* dataset's parameters. If the source dataset uses UTM or */
1178 : /* geographic coordinates, the coordinate system is carried over */
1179 : /* to the new DEM file's parameters. If the source dataset has a */
1180 : /* DEM compatible horizontal datum, the datum is carried over. */
1181 : /* Otherwise, the DEM dataset is configured to use geographic */
1182 : /* coordinates and a default datum. */
1183 : /* (Hunter Blanks, 8/31/04, hblanks@artifex.org) */
1184 : /************************************************************************/
1185 :
1186 34 : static int USGSDEMProductSetup_DEFAULT( USGSDEMWriteInfo *psWInfo )
1187 :
1188 : {
1189 :
1190 : /* -------------------------------------------------------------------- */
1191 : /* Set the destination coordinate system. */
1192 : /* -------------------------------------------------------------------- */
1193 34 : OGRSpatialReference DstoSRS;
1194 34 : OGRSpatialReference SrcoSRS;
1195 : char *sourceWkt;
1196 34 : int bNorth = TRUE;
1197 : /* XXX here we are assume (!) northern hemisphere UTM datasets */
1198 : char **readSourceWkt;
1199 : int i;
1200 34 : int numdatums = 4;
1201 34 : const char DatumCodes[4][2] = { "1", "2", "3", "4" };
1202 : char Datums[4][6] = { "NAD27", "WGS72", "WGS84",
1203 34 : "NAD83" };
1204 :
1205 : /* get the source dataset's projection */
1206 34 : sourceWkt = (char *) psWInfo->poSrcDS->GetProjectionRef();
1207 34 : readSourceWkt = &sourceWkt;
1208 34 : if (SrcoSRS.importFromWkt(readSourceWkt) != OGRERR_NONE)
1209 : {
1210 : CPLError( CE_Failure, CPLE_AppDefined,
1211 0 : "DEM Default Setup: Importing source dataset projection failed" );
1212 0 : return FALSE;
1213 : }
1214 :
1215 : /* Set the destination dataset's projection. If the source datum
1216 : * used is DEM compatible, just use it. Otherwise, default to the
1217 : * last datum in the Datums array.
1218 : */
1219 100 : for( i=0; i < numdatums; i++ )
1220 : {
1221 100 : if (DstoSRS.SetWellKnownGeogCS(Datums[i]) != OGRERR_NONE)
1222 : {
1223 : CPLError( CE_Failure, CPLE_AppDefined,
1224 0 : "DEM Default Setup: Failed to set datum of destination" );
1225 0 : return FALSE;
1226 : }
1227 : /* XXX Hopefully it's ok, to just keep changing the projection
1228 : * of our destination. If not, we'll want to reinitialize the
1229 : * OGRSpatialReference each time.
1230 : */
1231 100 : if ( DstoSRS.IsSameGeogCS( &SrcoSRS ) )
1232 : {
1233 34 : break;
1234 : }
1235 : }
1236 34 : strncpy( psWInfo->horizdatum, DatumCodes[i], 2 );
1237 :
1238 : /* get the UTM zone, if any */
1239 34 : psWInfo->utmzone = SrcoSRS.GetUTMZone(&bNorth);
1240 34 : if (psWInfo->utmzone)
1241 : {
1242 2 : if (DstoSRS.SetUTM(psWInfo->utmzone) != OGRERR_NONE)
1243 : {
1244 : CPLError( CE_Failure, CPLE_AppDefined,
1245 0 : "DEM Default Setup: Failed to set utm zone of destination" );
1246 : /* SetUTM isn't documented to return OGRERR_NONE
1247 : * on success, but it does, so, we'll check for it.
1248 : */
1249 0 : return FALSE;
1250 : }
1251 : }
1252 :
1253 : /* export the projection to sWInfo */
1254 34 : if (DstoSRS.exportToWkt( &(psWInfo->pszDstSRS) ) != OGRERR_NONE)
1255 : {
1256 : CPLError( CE_Failure, CPLE_AppDefined,
1257 0 : "UTMDEM: Failed to export destination Wkt to psWInfo" );
1258 : }
1259 34 : return TRUE;
1260 : }
1261 :
1262 : /************************************************************************/
1263 : /* USGSDEMLoadRaster() */
1264 : /* */
1265 : /* Loads the raster from the source dataset (not normally USGS */
1266 : /* DEM) into memory. If nodata is marked, a special effort is */
1267 : /* made to translate it properly into the USGS nodata value. */
1268 : /************************************************************************/
1269 :
1270 36 : static int USGSDEMLoadRaster( USGSDEMWriteInfo *psWInfo,
1271 : GDALRasterBand *poSrcBand )
1272 :
1273 : {
1274 : CPLErr eErr;
1275 : int i;
1276 :
1277 : /* -------------------------------------------------------------------- */
1278 : /* Allocate output array, and pre-initialize to NODATA value. */
1279 : /* -------------------------------------------------------------------- */
1280 : psWInfo->panData =
1281 36 : (GInt16 *) VSIMalloc3( 2, psWInfo->nXSize, psWInfo->nYSize );
1282 36 : if( psWInfo->panData == NULL )
1283 : {
1284 : CPLError( CE_Failure, CPLE_OutOfMemory,
1285 : "Out of memory allocating %d byte internal copy of DEM.",
1286 0 : 2 * psWInfo->nXSize * psWInfo->nYSize );
1287 0 : return FALSE;
1288 : }
1289 :
1290 8746368 : for( i = 0; i < psWInfo->nXSize * psWInfo->nYSize; i++ )
1291 8746332 : psWInfo->panData[i] = DEM_NODATA;
1292 :
1293 : /* -------------------------------------------------------------------- */
1294 : /* Make a "memory dataset" wrapper for this data array. */
1295 : /* -------------------------------------------------------------------- */
1296 36 : GDALDriver *poMemDriver = (GDALDriver *) GDALGetDriverByName( "MEM" );
1297 : GDALDataset *poMemDS;
1298 :
1299 36 : if( poMemDriver == NULL )
1300 : {
1301 : CPLError( CE_Failure, CPLE_AppDefined,
1302 0 : "Failed to find MEM driver." );
1303 0 : return FALSE;
1304 : }
1305 :
1306 : poMemDS =
1307 : poMemDriver->Create( "USGSDEM_temp", psWInfo->nXSize, psWInfo->nYSize,
1308 36 : 0, GDT_Int16, NULL );
1309 36 : if( poMemDS == NULL )
1310 0 : return FALSE;
1311 :
1312 : /* -------------------------------------------------------------------- */
1313 : /* Now add the array itself as a band. */
1314 : /* -------------------------------------------------------------------- */
1315 : char szDataPointer[100];
1316 36 : char *apszOptions[] = { szDataPointer, NULL };
1317 :
1318 36 : memset( szDataPointer, 0, sizeof(szDataPointer) );
1319 36 : sprintf( szDataPointer, "DATAPOINTER=" );
1320 : CPLPrintPointer( szDataPointer+strlen(szDataPointer),
1321 : psWInfo->panData,
1322 36 : sizeof(szDataPointer) - strlen(szDataPointer) );
1323 :
1324 36 : if( poMemDS->AddBand( GDT_Int16, apszOptions ) != CE_None )
1325 0 : return FALSE;
1326 :
1327 : /* -------------------------------------------------------------------- */
1328 : /* Assign geotransform and nodata indicators. */
1329 : /* -------------------------------------------------------------------- */
1330 : double adfGeoTransform[6];
1331 :
1332 36 : adfGeoTransform[0] = psWInfo->dfULX - psWInfo->dfHorizStepSize * 0.5;
1333 36 : adfGeoTransform[1] = psWInfo->dfHorizStepSize;
1334 36 : adfGeoTransform[2] = 0.0;
1335 36 : adfGeoTransform[3] = psWInfo->dfULY + psWInfo->dfVertStepSize * 0.5;
1336 36 : adfGeoTransform[4] = 0.0;
1337 36 : adfGeoTransform[5] = -psWInfo->dfVertStepSize;
1338 :
1339 36 : poMemDS->SetGeoTransform( adfGeoTransform );
1340 :
1341 : /* -------------------------------------------------------------------- */
1342 : /* Set coordinate system if we have a special one to set. */
1343 : /* -------------------------------------------------------------------- */
1344 36 : if( psWInfo->pszDstSRS )
1345 36 : poMemDS->SetProjection( psWInfo->pszDstSRS );
1346 :
1347 : /* -------------------------------------------------------------------- */
1348 : /* Establish the resampling kernel to use. */
1349 : /* -------------------------------------------------------------------- */
1350 36 : GDALResampleAlg eResampleAlg = GRA_Bilinear;
1351 : const char *pszResample = CSLFetchNameValue( psWInfo->papszOptions,
1352 36 : "RESAMPLE" );
1353 :
1354 36 : if( pszResample == NULL )
1355 : /* bilinear */;
1356 10 : else if( EQUAL(pszResample,"Nearest") )
1357 10 : eResampleAlg = GRA_NearestNeighbour;
1358 0 : else if( EQUAL(pszResample,"Bilinear") )
1359 0 : eResampleAlg = GRA_Bilinear;
1360 0 : else if( EQUAL(pszResample,"Cubic") )
1361 0 : eResampleAlg = GRA_Cubic;
1362 0 : else if( EQUAL(pszResample,"CubicSpline") )
1363 0 : eResampleAlg = GRA_CubicSpline;
1364 : else
1365 : {
1366 : CPLError( CE_Failure, CPLE_NotSupported,
1367 : "RESAMPLE=%s, not a supported resampling kernel.",
1368 0 : pszResample );
1369 0 : return FALSE;
1370 : }
1371 :
1372 : /* -------------------------------------------------------------------- */
1373 : /* Perform a warp from source dataset to destination buffer */
1374 : /* (memory dataset). */
1375 : /* -------------------------------------------------------------------- */
1376 : eErr = GDALReprojectImage( (GDALDatasetH) psWInfo->poSrcDS,
1377 36 : psWInfo->poSrcDS->GetProjectionRef(),
1378 : (GDALDatasetH) poMemDS,
1379 : psWInfo->pszDstSRS,
1380 : eResampleAlg, 0.0, 0.0, NULL, NULL,
1381 72 : NULL );
1382 :
1383 : /* -------------------------------------------------------------------- */
1384 : /* Deallocate memory wrapper for the buffer. */
1385 : /* -------------------------------------------------------------------- */
1386 36 : delete poMemDS;
1387 :
1388 36 : return eErr == CE_None;
1389 : }
1390 :
1391 :
1392 : /************************************************************************/
1393 : /* CreateCopy() */
1394 : /************************************************************************/
1395 :
1396 : GDALDataset *
1397 46 : USGSDEMCreateCopy( const char *pszFilename, GDALDataset *poSrcDS,
1398 : int bStrict, char **papszOptions,
1399 : GDALProgressFunc pfnProgress, void * pProgressData )
1400 :
1401 : {
1402 : USGSDEMWriteInfo sWInfo;
1403 :
1404 46 : if( poSrcDS->GetRasterCount() != 1 )
1405 : {
1406 : CPLError( CE_Failure, CPLE_AppDefined,
1407 10 : "Unable to create multi-band USGS DEM / CDED files." );
1408 10 : return NULL;
1409 : }
1410 :
1411 : /* -------------------------------------------------------------------- */
1412 : /* Capture some preliminary information. */
1413 : /* -------------------------------------------------------------------- */
1414 36 : memset( &sWInfo, 0, sizeof(sWInfo) );
1415 :
1416 36 : sWInfo.poSrcDS = poSrcDS;
1417 36 : sWInfo.pszFilename = CPLStrdup(pszFilename);
1418 36 : sWInfo.nXSize = poSrcDS->GetRasterXSize();
1419 36 : sWInfo.nYSize = poSrcDS->GetRasterYSize();
1420 36 : sWInfo.papszOptions = CSLDuplicate( papszOptions );
1421 36 : sWInfo.bStrict = bStrict;
1422 36 : sWInfo.utmzone = 0;
1423 36 : strncpy( sWInfo.horizdatum, "", 1 );
1424 :
1425 36 : if ( sWInfo.nXSize <= 1 || sWInfo.nYSize <= 1 )
1426 : {
1427 : CPLError( CE_Failure, CPLE_AppDefined,
1428 0 : "Source dataset dimensions must be at least 2x2." );
1429 0 : return NULL;
1430 : }
1431 :
1432 : /* -------------------------------------------------------------------- */
1433 : /* Work out corner coordinates. */
1434 : /* -------------------------------------------------------------------- */
1435 : double adfGeoTransform[6];
1436 :
1437 36 : poSrcDS->GetGeoTransform( adfGeoTransform );
1438 :
1439 36 : sWInfo.dfLLX = adfGeoTransform[0] + adfGeoTransform[1] * 0.5;
1440 36 : sWInfo.dfLLY = adfGeoTransform[3]
1441 36 : + adfGeoTransform[5] * (sWInfo.nYSize - 0.5);
1442 :
1443 36 : sWInfo.dfULX = adfGeoTransform[0] + adfGeoTransform[1] * 0.5;
1444 36 : sWInfo.dfULY = adfGeoTransform[3] + adfGeoTransform[5] * 0.5;
1445 :
1446 36 : sWInfo.dfURX = adfGeoTransform[0]
1447 36 : + adfGeoTransform[1] * (sWInfo.nXSize - 0.5);
1448 36 : sWInfo.dfURY = adfGeoTransform[3] + adfGeoTransform[5] * 0.5;
1449 :
1450 36 : sWInfo.dfLRX = adfGeoTransform[0]
1451 36 : + adfGeoTransform[1] * (sWInfo.nXSize - 0.5);
1452 36 : sWInfo.dfLRY = adfGeoTransform[3]
1453 36 : + adfGeoTransform[5] * (sWInfo.nYSize - 0.5);
1454 :
1455 36 : sWInfo.dfHorizStepSize = (sWInfo.dfURX - sWInfo.dfULX) / (sWInfo.nXSize-1);
1456 36 : sWInfo.dfVertStepSize = (sWInfo.dfURY - sWInfo.dfLRY) / (sWInfo.nYSize-1);
1457 :
1458 : /* -------------------------------------------------------------------- */
1459 : /* Allow override of z resolution, but default to 1.0. */
1460 : /* -------------------------------------------------------------------- */
1461 : const char *zResolution = CSLFetchNameValue(
1462 36 : sWInfo.papszOptions, "ZRESOLUTION" );
1463 :
1464 70 : if( zResolution == NULL || EQUAL(zResolution,"DEFAULT") )
1465 : {
1466 34 : sWInfo.dfElevStepSize = 1.0;
1467 : }
1468 : else
1469 : {
1470 : // XXX: We are using atof() here instead of CPLAtof() because
1471 : // zResolution value comes from user's input and supposed to be
1472 : // written according to user's current locale. atof() honors locale
1473 : // setting, CPLAtof() is not.
1474 2 : sWInfo.dfElevStepSize = atof( zResolution );
1475 2 : if ( sWInfo.dfElevStepSize <= 0 )
1476 : {
1477 : /* don't allow negative values */
1478 0 : sWInfo.dfElevStepSize = 1.0;
1479 : }
1480 : }
1481 :
1482 : /* -------------------------------------------------------------------- */
1483 : /* Initialize for special product configurations. */
1484 : /* -------------------------------------------------------------------- */
1485 : const char *pszProduct = CSLFetchNameValue( sWInfo.papszOptions,
1486 36 : "PRODUCT" );
1487 :
1488 70 : if( pszProduct == NULL || EQUAL(pszProduct,"DEFAULT") )
1489 : {
1490 34 : if ( !USGSDEMProductSetup_DEFAULT( &sWInfo ) )
1491 : {
1492 0 : USGSDEMWriteCleanup( &sWInfo );
1493 0 : return NULL;
1494 : }
1495 : }
1496 2 : else if( EQUAL(pszProduct,"CDED50K") )
1497 : {
1498 2 : if( !USGSDEMProductSetup_CDED50K( &sWInfo ) )
1499 : {
1500 0 : USGSDEMWriteCleanup( &sWInfo );
1501 0 : return NULL;
1502 : }
1503 : }
1504 : else
1505 : {
1506 : CPLError( CE_Failure, CPLE_NotSupported,
1507 : "DEM PRODUCT='%s' not recognised.",
1508 0 : pszProduct );
1509 0 : USGSDEMWriteCleanup( &sWInfo );
1510 0 : return NULL;
1511 : }
1512 :
1513 :
1514 : /* -------------------------------------------------------------------- */
1515 : /* Read the whole area of interest into memory. */
1516 : /* -------------------------------------------------------------------- */
1517 36 : if( !USGSDEMLoadRaster( &sWInfo, poSrcDS->GetRasterBand( 1 ) ) )
1518 : {
1519 2 : USGSDEMWriteCleanup( &sWInfo );
1520 2 : return NULL;
1521 : }
1522 :
1523 : /* -------------------------------------------------------------------- */
1524 : /* Create the output file. */
1525 : /* -------------------------------------------------------------------- */
1526 34 : sWInfo.fp = VSIFOpenL( pszFilename, "wb" );
1527 34 : if( sWInfo.fp == NULL )
1528 : {
1529 : CPLError( CE_Failure, CPLE_OpenFailed,
1530 2 : "%s", VSIStrerror( errno ) );
1531 2 : USGSDEMWriteCleanup( &sWInfo );
1532 2 : return NULL;
1533 : }
1534 :
1535 : /* -------------------------------------------------------------------- */
1536 : /* Write the A record. */
1537 : /* -------------------------------------------------------------------- */
1538 32 : if( !USGSDEMWriteARecord( &sWInfo ) )
1539 : {
1540 0 : USGSDEMWriteCleanup( &sWInfo );
1541 0 : return NULL;
1542 : }
1543 :
1544 : /* -------------------------------------------------------------------- */
1545 : /* Write profiles. */
1546 : /* -------------------------------------------------------------------- */
1547 : int iProfile;
1548 :
1549 3384 : for( iProfile = 0; iProfile < sWInfo.nXSize; iProfile++ )
1550 : {
1551 3352 : if( !USGSDEMWriteProfile( &sWInfo, iProfile ) )
1552 : {
1553 0 : USGSDEMWriteCleanup( &sWInfo );
1554 0 : return NULL;
1555 : }
1556 : }
1557 :
1558 : /* -------------------------------------------------------------------- */
1559 : /* Cleanup. */
1560 : /* -------------------------------------------------------------------- */
1561 32 : USGSDEMWriteCleanup( &sWInfo );
1562 :
1563 : /* -------------------------------------------------------------------- */
1564 : /* Re-open dataset, and copy any auxilary pam information. */
1565 : /* -------------------------------------------------------------------- */
1566 : GDALPamDataset *poDS = (GDALPamDataset *)
1567 32 : GDALOpen( pszFilename, GA_ReadOnly );
1568 :
1569 32 : if( poDS )
1570 32 : poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
1571 :
1572 32 : return poDS;
1573 : }
|
