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