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