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