1 : /******************************************************************************
2 : * $Id: dted_api.c 24603 2012-06-20 18:21:29Z rouault $
3 : *
4 : * Project: DTED Translator
5 : * Purpose: Implementation of DTED/CDED access functions.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 1999, Frank Warmerdam
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "dted_api.h"
31 :
32 : #ifndef AVOID_CPL
33 : CPL_CVSID("$Id: dted_api.c 24603 2012-06-20 18:21:29Z rouault $");
34 : #endif
35 :
36 : static int bWarnedTwoComplement = FALSE;
37 :
38 : static void DTEDDetectVariantWithMissingColumns(DTEDInfo* psDInfo);
39 :
40 : /************************************************************************/
41 : /* DTEDGetField() */
42 : /* */
43 : /* Extract a field as a zero terminated string. Address is */
44 : /* deliberately 1 based so the getfield arguments will be the */
45 : /* same as the numbers in the file format specification. */
46 : /************************************************************************/
47 :
48 : static
49 860 : char *DTEDGetField( char szResult[81], const char *pachRecord, int nStart, int nSize )
50 :
51 : {
52 860 : CPLAssert( nSize < 81 );
53 860 : memcpy( szResult, pachRecord + nStart - 1, nSize );
54 860 : szResult[nSize] = '\0';
55 :
56 860 : return szResult;
57 : }
58 :
59 : /************************************************************************/
60 : /* StripLeadingZeros() */
61 : /* */
62 : /* Return a pointer to the first non-zero character in BUF. */
63 : /* BUF must be null terminated. */
64 : /* If buff is all zeros, then it will point to the last non-zero */
65 : /************************************************************************/
66 :
67 516 : static const char* stripLeadingZeros(const char* buf)
68 : {
69 516 : const char* ptr = buf;
70 :
71 : /* Go until we run out of characters or hit something non-zero */
72 :
73 1569 : while( *ptr == '0' && *(ptr+1) != '\0' )
74 : {
75 537 : ptr++;
76 : }
77 :
78 516 : return ptr;
79 : }
80 :
81 : /************************************************************************/
82 : /* DTEDOpen() */
83 : /************************************************************************/
84 :
85 86 : DTEDInfo * DTEDOpen( const char * pszFilename,
86 : const char * pszAccess,
87 : int bTestOpen )
88 :
89 : {
90 : VSILFILE *fp;
91 : char achRecord[DTED_UHL_SIZE];
92 86 : DTEDInfo *psDInfo = NULL;
93 : double dfLLOriginX, dfLLOriginY;
94 86 : int deg = 0;
95 86 : int min = 0;
96 86 : int sec = 0;
97 86 : int bSwapLatLong = FALSE;
98 : char szResult[81];
99 : int bIsWeirdDTED;
100 : char chHemisphere;
101 : /* -------------------------------------------------------------------- */
102 : /* Open the physical file. */
103 : /* -------------------------------------------------------------------- */
104 158 : if( EQUAL(pszAccess,"r") || EQUAL(pszAccess,"rb") )
105 72 : pszAccess = "rb";
106 : else
107 14 : pszAccess = "r+b";
108 :
109 86 : fp = VSIFOpenL( pszFilename, pszAccess );
110 :
111 86 : if( fp == NULL )
112 : {
113 0 : if( !bTestOpen )
114 : {
115 : #ifndef AVOID_CPL
116 0 : CPLError( CE_Failure, CPLE_OpenFailed,
117 : #else
118 : fprintf( stderr,
119 : #endif
120 : "Failed to open file %s.",
121 : pszFilename );
122 : }
123 :
124 0 : return NULL;
125 : }
126 :
127 : /* -------------------------------------------------------------------- */
128 : /* Read, trying to find the UHL record. Skip VOL or HDR */
129 : /* records if they are encountered. */
130 : /* -------------------------------------------------------------------- */
131 : do
132 : {
133 87 : if( VSIFReadL( achRecord, 1, DTED_UHL_SIZE, fp ) != DTED_UHL_SIZE )
134 : {
135 0 : if( !bTestOpen )
136 : {
137 : #ifndef AVOID_CPL
138 0 : CPLError( CE_Failure, CPLE_OpenFailed,
139 : #else
140 : fprintf( stderr,
141 : #endif
142 : "Unable to read header, %s is not DTED.",
143 : pszFilename );
144 : }
145 0 : VSIFCloseL( fp );
146 0 : return NULL;
147 : }
148 :
149 87 : } while( EQUALN(achRecord,"VOL",3) || EQUALN(achRecord,"HDR",3) );
150 :
151 86 : if( !EQUALN(achRecord,"UHL",3) )
152 : {
153 0 : if( !bTestOpen )
154 : {
155 : #ifndef AVOID_CPL
156 0 : CPLError( CE_Failure, CPLE_OpenFailed,
157 : #else
158 : fprintf( stderr,
159 : #endif
160 : "No UHL record. %s is not a DTED file.",
161 : pszFilename );
162 : }
163 0 : VSIFCloseL( fp );
164 0 : return NULL;
165 : }
166 :
167 : /* -------------------------------------------------------------------- */
168 : /* Create and initialize the DTEDInfo structure. */
169 : /* -------------------------------------------------------------------- */
170 86 : psDInfo = (DTEDInfo *) CPLCalloc(1,sizeof(DTEDInfo));
171 :
172 86 : psDInfo->fp = fp;
173 :
174 86 : psDInfo->bUpdate = EQUAL(pszAccess,"r+b");
175 86 : psDInfo->bRewriteHeaders = FALSE;
176 :
177 86 : psDInfo->nUHLOffset = (int)VSIFTellL( fp ) - DTED_UHL_SIZE;
178 86 : psDInfo->pachUHLRecord = (char *) CPLMalloc(DTED_UHL_SIZE);
179 86 : memcpy( psDInfo->pachUHLRecord, achRecord, DTED_UHL_SIZE );
180 :
181 86 : psDInfo->nDSIOffset = (int)VSIFTellL( fp );
182 86 : psDInfo->pachDSIRecord = (char *) CPLMalloc(DTED_DSI_SIZE);
183 86 : VSIFReadL( psDInfo->pachDSIRecord, 1, DTED_DSI_SIZE, fp );
184 :
185 86 : psDInfo->nACCOffset = (int)VSIFTellL( fp );
186 86 : psDInfo->pachACCRecord = (char *) CPLMalloc(DTED_ACC_SIZE);
187 86 : VSIFReadL( psDInfo->pachACCRecord, 1, DTED_ACC_SIZE, fp );
188 :
189 172 : if( !EQUALN(psDInfo->pachDSIRecord,"DSI",3)
190 86 : || !EQUALN(psDInfo->pachACCRecord,"ACC",3) )
191 : {
192 : #ifndef AVOID_CPL
193 0 : CPLError( CE_Failure, CPLE_OpenFailed,
194 : #else
195 : fprintf( stderr,
196 : #endif
197 : "DSI or ACC record missing. DTED access to\n%s failed.",
198 : pszFilename );
199 :
200 0 : DTEDClose(psDInfo);
201 0 : return NULL;
202 : }
203 :
204 86 : psDInfo->nDataOffset = (int)VSIFTellL( fp );
205 :
206 : /* DTED3 file from http://www.falconview.org/trac/FalconView/downloads/20 */
207 : /* (co_elevation.zip) has really weird offsets that don't comply with the 89020B specification */
208 86 : bIsWeirdDTED = achRecord[4] == ' ';
209 :
210 : /* -------------------------------------------------------------------- */
211 : /* Parse out position information. Note that we are extracting */
212 : /* the top left corner of the top left pixel area, not the */
213 : /* center of the area. */
214 : /* -------------------------------------------------------------------- */
215 86 : if (!bIsWeirdDTED)
216 : {
217 86 : psDInfo->dfPixelSizeX =
218 86 : atoi(DTEDGetField(szResult,achRecord,21,4)) / 36000.0;
219 :
220 86 : psDInfo->dfPixelSizeY =
221 86 : atoi(DTEDGetField(szResult,achRecord,25,4)) / 36000.0;
222 :
223 86 : psDInfo->nXSize = atoi(DTEDGetField(szResult,achRecord,48,4));
224 86 : psDInfo->nYSize = atoi(DTEDGetField(szResult,achRecord,52,4));
225 : }
226 : else
227 : {
228 0 : psDInfo->dfPixelSizeX =
229 0 : atoi(DTEDGetField(szResult,achRecord,41,4)) / 36000.0;
230 :
231 0 : psDInfo->dfPixelSizeY =
232 0 : atoi(DTEDGetField(szResult,achRecord,45,4)) / 36000.0;
233 :
234 0 : psDInfo->nXSize = atoi(DTEDGetField(szResult,psDInfo->pachDSIRecord,563,4));
235 0 : psDInfo->nYSize = atoi(DTEDGetField(szResult,psDInfo->pachDSIRecord,567,4));
236 : }
237 :
238 86 : if (psDInfo->nXSize <= 0 || psDInfo->nYSize <= 0)
239 : {
240 : #ifndef AVOID_CPL
241 0 : CPLError( CE_Failure, CPLE_OpenFailed,
242 : #else
243 : fprintf( stderr,
244 : #endif
245 : "Invalid dimensions : %d x %d. DTED access to\n%s failed.",
246 : psDInfo->nXSize, psDInfo->nYSize, pszFilename );
247 :
248 0 : DTEDClose(psDInfo);
249 0 : return NULL;
250 : }
251 :
252 : /* create a scope so I don't need to declare these up top */
253 86 : if (!bIsWeirdDTED)
254 : {
255 86 : deg = atoi(stripLeadingZeros(DTEDGetField(szResult,achRecord,5,3)));
256 86 : min = atoi(stripLeadingZeros(DTEDGetField(szResult,achRecord,8,2)));
257 86 : sec = atoi(stripLeadingZeros(DTEDGetField(szResult,achRecord,10,2)));
258 86 : chHemisphere = achRecord[11];
259 : }
260 : else
261 : {
262 0 : deg = atoi(stripLeadingZeros(DTEDGetField(szResult,achRecord,9,3)));
263 0 : min = atoi(stripLeadingZeros(DTEDGetField(szResult,achRecord,12,2)));
264 0 : sec = atoi(stripLeadingZeros(DTEDGetField(szResult,achRecord,14,2)));
265 0 : chHemisphere = achRecord[15];
266 : }
267 :
268 : /* NOTE : The first version of MIL-D-89020 was buggy.
269 : The latitude and longitude of the LL cornder of the UHF record was inverted.
270 : This was fixed in MIL-D-89020 Amendement 1, but some products may be affected.
271 : We detect this situation by looking at N/S in the longitude field and
272 : E/W in the latitude one
273 : */
274 :
275 86 : dfLLOriginX = deg + min / 60.0 + sec / 3600.0;
276 86 : if( chHemisphere == 'W' )
277 63 : dfLLOriginX *= -1;
278 23 : else if ( chHemisphere == 'N' )
279 1 : bSwapLatLong = TRUE;
280 22 : else if ( chHemisphere == 'S' )
281 : {
282 0 : dfLLOriginX *= -1;
283 0 : bSwapLatLong = TRUE;
284 : }
285 :
286 86 : if (!bIsWeirdDTED)
287 : {
288 86 : deg = atoi(stripLeadingZeros(DTEDGetField(szResult,achRecord,13,3)));
289 86 : min = atoi(stripLeadingZeros(DTEDGetField(szResult,achRecord,16,2)));
290 86 : sec = atoi(stripLeadingZeros(DTEDGetField(szResult,achRecord,18,2)));
291 86 : chHemisphere = achRecord[19];
292 : }
293 : else
294 : {
295 0 : deg = atoi(stripLeadingZeros(DTEDGetField(szResult,achRecord,25,3)));
296 0 : min = atoi(stripLeadingZeros(DTEDGetField(szResult,achRecord,28,2)));
297 0 : sec = atoi(stripLeadingZeros(DTEDGetField(szResult,achRecord,30,2)));
298 0 : chHemisphere = achRecord[31];
299 : }
300 :
301 86 : dfLLOriginY = deg + min / 60.0 + sec / 3600.0;
302 86 : if( chHemisphere == 'S' || (bSwapLatLong && chHemisphere == 'W'))
303 1 : dfLLOriginY *= -1;
304 :
305 86 : if (bSwapLatLong)
306 : {
307 1 : double dfTmp = dfLLOriginX;
308 1 : dfLLOriginX = dfLLOriginY;
309 1 : dfLLOriginY = dfTmp;
310 : }
311 :
312 86 : psDInfo->dfULCornerX = dfLLOriginX - 0.5 * psDInfo->dfPixelSizeX;
313 86 : psDInfo->dfULCornerY = dfLLOriginY - 0.5 * psDInfo->dfPixelSizeY
314 86 : + psDInfo->nYSize * psDInfo->dfPixelSizeY;
315 :
316 86 : DTEDDetectVariantWithMissingColumns(psDInfo);
317 :
318 86 : return psDInfo;
319 : }
320 :
321 : /************************************************************************/
322 : /* DTEDDetectVariantWithMissingColumns() */
323 : /************************************************************************/
324 :
325 86 : static void DTEDDetectVariantWithMissingColumns(DTEDInfo* psDInfo)
326 : {
327 : /* -------------------------------------------------------------------- */
328 : /* Some DTED files have only a subset of all possible columns. */
329 : /* They can declare for example 3601 columns, but in the file, */
330 : /* there are just columns 100->500. Detect that situation. */
331 : /* -------------------------------------------------------------------- */
332 :
333 : GByte pabyRecordHeader[8];
334 : int nFirstDataBlockCount, nFirstLongitudeCount;
335 : int nLastDataBlockCount, nLastLongitudeCount;
336 : int nSize;
337 86 : int nColByteSize = 12 + psDInfo->nYSize*2;
338 :
339 86 : VSIFSeekL(psDInfo->fp, psDInfo->nDataOffset, SEEK_SET);
340 170 : if (VSIFReadL(pabyRecordHeader, 1, 8, psDInfo->fp) != 8 ||
341 84 : pabyRecordHeader[0] != 0252)
342 : {
343 2 : CPLDebug("DTED", "Cannot find signature of first column");
344 2 : return;
345 : }
346 :
347 84 : nFirstDataBlockCount = (pabyRecordHeader[2] << 8) | pabyRecordHeader[3];
348 84 : nFirstLongitudeCount = (pabyRecordHeader[4] << 8) | pabyRecordHeader[5];
349 :
350 84 : VSIFSeekL(psDInfo->fp, 0, SEEK_END);
351 84 : nSize = (int)VSIFTellL(psDInfo->fp);
352 84 : if (nSize < 12 + psDInfo->nYSize*2)
353 : {
354 0 : CPLDebug("DTED", "File too short");
355 0 : return;
356 : }
357 :
358 84 : VSIFSeekL(psDInfo->fp, nSize - nColByteSize, SEEK_SET);
359 168 : if (VSIFReadL(pabyRecordHeader, 1, 8, psDInfo->fp) != 8 ||
360 84 : pabyRecordHeader[0] != 0252)
361 : {
362 0 : CPLDebug("DTED", "Cannot find signature of last column");
363 0 : return;
364 : }
365 :
366 84 : nLastDataBlockCount = (pabyRecordHeader[2] << 8) | pabyRecordHeader[3];
367 84 : nLastLongitudeCount = (pabyRecordHeader[4] << 8) | pabyRecordHeader[5];
368 :
369 330 : if (nFirstDataBlockCount == 0 &&
370 : nFirstLongitudeCount == 0 &&
371 82 : nLastDataBlockCount == psDInfo->nXSize - 1 &&
372 82 : nLastLongitudeCount == psDInfo->nXSize - 1 &&
373 82 : nSize - psDInfo->nDataOffset == psDInfo->nXSize * nColByteSize)
374 : {
375 : /* This is the most standard form of DTED. Return happily now. */
376 82 : return;
377 : }
378 :
379 : /* Well, we have an odd DTED file at that point */
380 :
381 2 : psDInfo->panMapLogicalColsToOffsets =
382 2 : (int*)CPLMalloc(psDInfo->nXSize * sizeof(int));
383 :
384 6 : if (nFirstDataBlockCount == 0 &&
385 2 : nLastLongitudeCount - nFirstLongitudeCount ==
386 2 : nLastDataBlockCount - nFirstDataBlockCount &&
387 1 : nSize - psDInfo->nDataOffset ==
388 1 : (nLastLongitudeCount - nFirstLongitudeCount + 1) * nColByteSize)
389 : {
390 : int i;
391 :
392 : /* Case seen in a real-world file */
393 :
394 1 : CPLDebug("DTED", "The file only contains data from column %d to column %d.",
395 : nFirstLongitudeCount, nLastLongitudeCount);
396 :
397 122 : for(i = 0; i < psDInfo->nXSize; i++)
398 : {
399 121 : if (i < nFirstLongitudeCount)
400 2 : psDInfo->panMapLogicalColsToOffsets[i] = -1;
401 119 : else if (i <= nLastLongitudeCount)
402 4 : psDInfo->panMapLogicalColsToOffsets[i] =
403 2 : psDInfo->nDataOffset + (i - nFirstLongitudeCount) * nColByteSize;
404 : else
405 117 : psDInfo->panMapLogicalColsToOffsets[i] = -1;
406 : }
407 : }
408 : else
409 : {
410 1 : int nPhysicalCols = (nSize - psDInfo->nDataOffset) / nColByteSize;
411 : int i;
412 :
413 : /* Theoretical case for now... */
414 :
415 1 : CPLDebug("DTED", "There columns appear to be in non sequential order. "
416 : "Scanning the whole file.");
417 :
418 122 : for(i = 0; i < psDInfo->nXSize; i++)
419 : {
420 121 : psDInfo->panMapLogicalColsToOffsets[i] = -1;
421 : }
422 :
423 3 : for(i = 0; i < nPhysicalCols; i++)
424 : {
425 : int nDataBlockCount, nLongitudeCount;
426 :
427 2 : VSIFSeekL(psDInfo->fp, psDInfo->nDataOffset + i * nColByteSize, SEEK_SET);
428 4 : if (VSIFReadL(pabyRecordHeader, 1, 8, psDInfo->fp) != 8 ||
429 2 : pabyRecordHeader[0] != 0252)
430 : {
431 0 : CPLDebug("DTED", "Cannot find signature of physical column %d", i);
432 0 : return;
433 : }
434 :
435 2 : nDataBlockCount = (pabyRecordHeader[2] << 8) | pabyRecordHeader[3];
436 2 : if (nDataBlockCount != i)
437 : {
438 0 : CPLDebug("DTED", "Unexpected block count(%d) at physical column %d. "
439 : "Ignoring that and going on...",
440 : nDataBlockCount, i);
441 : }
442 :
443 2 : nLongitudeCount = (pabyRecordHeader[4] << 8) | pabyRecordHeader[5];
444 2 : if (nLongitudeCount < 0 || nLongitudeCount >= psDInfo->nXSize)
445 : {
446 0 : CPLDebug("DTED", "Invalid longitude count (%d) at physical column %d",
447 : nLongitudeCount, i);
448 0 : return;
449 : }
450 :
451 4 : psDInfo->panMapLogicalColsToOffsets[nLongitudeCount] =
452 2 : psDInfo->nDataOffset + i * nColByteSize;
453 : }
454 : }
455 : }
456 :
457 : /************************************************************************/
458 : /* DTEDReadPoint() */
459 : /* */
460 : /* Read one single sample. The coordinates are given from the */
461 : /* top-left corner of the file (contrary to the internal */
462 : /* organisation or a DTED file) */
463 : /************************************************************************/
464 :
465 0 : int DTEDReadPoint( DTEDInfo * psDInfo, int nXOff, int nYOff, GInt16* panVal)
466 : {
467 : int nOffset;
468 : GByte pabyData[2];
469 :
470 0 : if (nYOff < 0 || nXOff < 0 || nYOff >= psDInfo->nYSize || nXOff >= psDInfo->nXSize)
471 : {
472 : #ifndef AVOID_CPL
473 0 : CPLError( CE_Failure, CPLE_AppDefined,
474 : #else
475 : fprintf( stderr,
476 : #endif
477 : "Invalid raster coordinates (%d,%d) in DTED file.\n", nXOff, nYOff);
478 0 : return FALSE;
479 : }
480 :
481 0 : if (psDInfo->panMapLogicalColsToOffsets != NULL)
482 : {
483 0 : nOffset = psDInfo->panMapLogicalColsToOffsets[nXOff];
484 0 : if( nOffset < 0 )
485 : {
486 0 : *panVal = DTED_NODATA_VALUE;
487 0 : return TRUE;
488 : }
489 : }
490 : else
491 0 : nOffset = psDInfo->nDataOffset + nXOff * (12+psDInfo->nYSize*2);
492 0 : nOffset += 8 + 2 * (psDInfo->nYSize-1-nYOff);
493 :
494 0 : if( VSIFSeekL( psDInfo->fp, nOffset, SEEK_SET ) != 0
495 0 : || VSIFReadL( pabyData, 2, 1, psDInfo->fp ) != 1)
496 : {
497 : #ifndef AVOID_CPL
498 0 : CPLError( CE_Failure, CPLE_FileIO,
499 : #else
500 : fprintf( stderr,
501 : #endif
502 : "Failed to seek to, or read (%d,%d) at offset %d\n"
503 : "in DTED file.\n",
504 : nXOff, nYOff, nOffset );
505 0 : return FALSE;
506 : }
507 :
508 0 : *panVal = ((pabyData[0] & 0x7f) << 8) | pabyData[1];
509 :
510 0 : if( pabyData[0] & 0x80 )
511 : {
512 0 : *panVal *= -1;
513 :
514 : /*
515 : ** It seems that some files are improperly generated in twos
516 : ** complement form for negatives. For these, redo the job
517 : ** in twos complement. eg. w_069_s50.dt0
518 : */
519 0 : if(( *panVal < -16000 ) && (*panVal != DTED_NODATA_VALUE))
520 : {
521 0 : *panVal = (pabyData[0] << 8) | pabyData[1];
522 :
523 0 : if( !bWarnedTwoComplement )
524 : {
525 0 : bWarnedTwoComplement = TRUE;
526 : #ifndef AVOID_CPL
527 0 : CPLError( CE_Warning, CPLE_AppDefined,
528 : #else
529 : fprintf( stderr,
530 : #endif
531 : "The DTED driver found values less than -16000, and has adjusted\n"
532 : "them assuming they are improperly two-complemented. No more warnings\n"
533 : "will be issued in this session about this operation." );
534 : }
535 : }
536 : }
537 :
538 0 : return TRUE;
539 : }
540 :
541 : /************************************************************************/
542 : /* DTEDReadProfile() */
543 : /* */
544 : /* Read one profile line. These are organized in bottom to top */
545 : /* order starting from the leftmost column (0). */
546 : /************************************************************************/
547 :
548 0 : int DTEDReadProfile( DTEDInfo * psDInfo, int nColumnOffset,
549 : GInt16 * panData )
550 : {
551 0 : return DTEDReadProfileEx( psDInfo, nColumnOffset, panData, FALSE);
552 : }
553 :
554 4994 : int DTEDReadProfileEx( DTEDInfo * psDInfo, int nColumnOffset,
555 : GInt16 * panData, int bVerifyChecksum )
556 : {
557 : int nOffset;
558 : int i;
559 : GByte *pabyRecord;
560 : int nLongitudeCount;
561 :
562 : /* -------------------------------------------------------------------- */
563 : /* Read data record from disk. */
564 : /* -------------------------------------------------------------------- */
565 4994 : if (psDInfo->panMapLogicalColsToOffsets != NULL)
566 : {
567 242 : nOffset = psDInfo->panMapLogicalColsToOffsets[nColumnOffset];
568 242 : if( nOffset < 0 )
569 : {
570 29036 : for( i = 0; i < psDInfo->nYSize; i++ )
571 : {
572 28798 : panData[i] = DTED_NODATA_VALUE;
573 : }
574 238 : return TRUE;
575 : }
576 : }
577 : else
578 4752 : nOffset = psDInfo->nDataOffset + nColumnOffset * (12+psDInfo->nYSize*2);
579 :
580 4756 : pabyRecord = (GByte *) CPLMalloc(12 + psDInfo->nYSize*2);
581 :
582 9512 : if( VSIFSeekL( psDInfo->fp, nOffset, SEEK_SET ) != 0
583 9512 : || VSIFReadL( pabyRecord, (12+psDInfo->nYSize*2), 1, psDInfo->fp ) != 1)
584 : {
585 : #ifndef AVOID_CPL
586 0 : CPLError( CE_Failure, CPLE_FileIO,
587 : #else
588 : fprintf( stderr,
589 : #endif
590 : "Failed to seek to, or read profile %d at offset %d\n"
591 : "in DTED file.\n",
592 : nColumnOffset, nOffset );
593 0 : CPLFree( pabyRecord );
594 0 : return FALSE;
595 : }
596 :
597 4756 : nLongitudeCount = (pabyRecord[4] << 8) | pabyRecord[5];
598 4756 : if( nLongitudeCount != nColumnOffset )
599 : {
600 : #ifndef AVOID_CPL
601 0 : CPLError( CE_Warning, CPLE_AppDefined,
602 : #else
603 : fprintf( stderr,
604 : #endif
605 : "Longitude count (%d) of column %d doesn't match expected value.\n",
606 : nLongitudeCount, nColumnOffset );
607 : }
608 :
609 : /* -------------------------------------------------------------------- */
610 : /* Translate data values from "signed magnitude" to standard */
611 : /* binary. */
612 : /* -------------------------------------------------------------------- */
613 1229312 : for( i = 0; i < psDInfo->nYSize; i++ )
614 : {
615 1224556 : panData[i] = ((pabyRecord[8+i*2] & 0x7f) << 8) | pabyRecord[8+i*2+1];
616 :
617 1224556 : if( pabyRecord[8+i*2] & 0x80 )
618 : {
619 0 : panData[i] *= -1;
620 :
621 : /*
622 : ** It seems that some files are improperly generated in twos
623 : ** complement form for negatives. For these, redo the job
624 : ** in twos complement. eg. w_069_s50.dt0
625 : */
626 0 : if(( panData[i] < -16000 ) && (panData[i] != DTED_NODATA_VALUE))
627 : {
628 0 : panData[i] = (pabyRecord[8+i*2] << 8) | pabyRecord[8+i*2+1];
629 :
630 0 : if( !bWarnedTwoComplement )
631 : {
632 0 : bWarnedTwoComplement = TRUE;
633 : #ifndef AVOID_CPL
634 0 : CPLError( CE_Warning, CPLE_AppDefined,
635 : #else
636 : fprintf( stderr,
637 : #endif
638 : "The DTED driver found values less than -16000, and has adjusted\n"
639 : "them assuming they are improperly two-complemented. No more warnings\n"
640 : "will be issued in this session about this operation." );
641 : }
642 : }
643 : }
644 : }
645 :
646 4756 : if (bVerifyChecksum)
647 : {
648 1 : unsigned int nCheckSum = 0;
649 : unsigned int fileCheckSum;
650 :
651 : /* -------------------------------------------------------------------- */
652 : /* Verify the checksum. */
653 : /* -------------------------------------------------------------------- */
654 :
655 251 : for( i = 0; i < psDInfo->nYSize*2 + 8; i++ )
656 250 : nCheckSum += pabyRecord[i];
657 :
658 2 : fileCheckSum = (pabyRecord[8+psDInfo->nYSize*2+0] << 24) |
659 1 : (pabyRecord[8+psDInfo->nYSize*2+1] << 16) |
660 1 : (pabyRecord[8+psDInfo->nYSize*2+2] << 8) |
661 1 : pabyRecord[8+psDInfo->nYSize*2+3];
662 :
663 1 : if ((GIntBig)fileCheckSum > (GIntBig)(0xff * (8+psDInfo->nYSize*2)))
664 : {
665 : static int bWarned = FALSE;
666 0 : if (! bWarned)
667 : {
668 0 : bWarned = TRUE;
669 : #ifndef AVOID_CPL
670 0 : CPLError( CE_Warning, CPLE_AppDefined,
671 : #else
672 : fprintf( stderr,
673 : #endif
674 : "The DTED driver has read from the file a checksum "
675 : "with an impossible value (0x%X) at column %d.\n"
676 : "Check with your file producer.\n"
677 : "No more warnings will be issued in this session about this operation.",
678 : fileCheckSum, nColumnOffset);
679 : }
680 : }
681 1 : else if (fileCheckSum != nCheckSum)
682 : {
683 : #ifndef AVOID_CPL
684 1 : CPLError( CE_Warning, CPLE_AppDefined,
685 : #else
686 : fprintf( stderr,
687 : #endif
688 : "The DTED driver has found a computed and read checksum "
689 : "that do not match at column %d. Computed 0x%X, read 0x%X\n",
690 : nColumnOffset, nCheckSum, fileCheckSum);
691 1 : CPLFree( pabyRecord );
692 1 : return FALSE;
693 : }
694 : }
695 :
696 4755 : CPLFree( pabyRecord );
697 :
698 4755 : return TRUE;
699 : }
700 :
701 : /************************************************************************/
702 : /* DTEDWriteProfile() */
703 : /************************************************************************/
704 :
705 14054 : int DTEDWriteProfile( DTEDInfo * psDInfo, int nColumnOffset,
706 : GInt16 * panData )
707 :
708 : {
709 : int nOffset;
710 14054 : int i, nCheckSum = 0;
711 : GByte *pabyRecord;
712 :
713 14054 : if (psDInfo->panMapLogicalColsToOffsets != NULL)
714 : {
715 : #ifndef AVOID_CPL
716 0 : CPLError( CE_Failure, CPLE_NotSupported,
717 : #else
718 : fprintf( stderr,
719 : #endif
720 : "Write to partial file not supported.\n");
721 0 : return FALSE;
722 : }
723 :
724 : /* -------------------------------------------------------------------- */
725 : /* Format the data record. */
726 : /* -------------------------------------------------------------------- */
727 14054 : pabyRecord = (GByte *) CPLMalloc(12 + psDInfo->nYSize*2);
728 :
729 16631548 : for( i = 0; i < psDInfo->nYSize; i++ )
730 : {
731 16617494 : int nABSVal = ABS(panData[psDInfo->nYSize-i-1]);
732 16617494 : pabyRecord[8+i*2] = (GByte) ((nABSVal >> 8) & 0x7f);
733 16617494 : pabyRecord[8+i*2+1] = (GByte) (nABSVal & 0xff);
734 :
735 16617494 : if( panData[psDInfo->nYSize-i-1] < 0 )
736 2836848 : pabyRecord[8+i*2] |= 0x80;
737 : }
738 :
739 14054 : pabyRecord[0] = 0xaa;
740 14054 : pabyRecord[1] = 0;
741 14054 : pabyRecord[2] = (GByte) (nColumnOffset / 256);
742 14054 : pabyRecord[3] = (GByte) (nColumnOffset % 256);
743 14054 : pabyRecord[4] = (GByte) (nColumnOffset / 256);
744 14054 : pabyRecord[5] = (GByte) (nColumnOffset % 256);
745 14054 : pabyRecord[6] = 0;
746 14054 : pabyRecord[7] = 0;
747 :
748 : /* -------------------------------------------------------------------- */
749 : /* Compute the checksum. */
750 : /* -------------------------------------------------------------------- */
751 33361474 : for( i = 0; i < psDInfo->nYSize*2 + 8; i++ )
752 33347420 : nCheckSum += pabyRecord[i];
753 :
754 14054 : pabyRecord[8+psDInfo->nYSize*2+0] = (GByte) ((nCheckSum >> 24) & 0xff);
755 14054 : pabyRecord[8+psDInfo->nYSize*2+1] = (GByte) ((nCheckSum >> 16) & 0xff);
756 14054 : pabyRecord[8+psDInfo->nYSize*2+2] = (GByte) ((nCheckSum >> 8) & 0xff);
757 14054 : pabyRecord[8+psDInfo->nYSize*2+3] = (GByte) (nCheckSum & 0xff);
758 :
759 : /* -------------------------------------------------------------------- */
760 : /* Write the record. */
761 : /* -------------------------------------------------------------------- */
762 14054 : nOffset = psDInfo->nDataOffset + nColumnOffset * (12+psDInfo->nYSize*2);
763 :
764 28108 : if( VSIFSeekL( psDInfo->fp, nOffset, SEEK_SET ) != 0
765 28108 : || VSIFWriteL( pabyRecord,(12+psDInfo->nYSize*2),1,psDInfo->fp ) != 1)
766 : {
767 : #ifndef AVOID_CPL
768 0 : CPLError( CE_Failure, CPLE_FileIO,
769 : #else
770 : fprintf( stderr,
771 : #endif
772 : "Failed to seek to, or write profile %d at offset %d\n"
773 : "in DTED file.\n",
774 : nColumnOffset, nOffset );
775 0 : CPLFree( pabyRecord );
776 0 : return FALSE;
777 : }
778 :
779 14054 : CPLFree( pabyRecord );
780 :
781 14054 : return TRUE;
782 :
783 : }
784 :
785 : /************************************************************************/
786 : /* DTEDGetMetadataLocation() */
787 : /************************************************************************/
788 :
789 1708 : static void DTEDGetMetadataLocation( DTEDInfo *psDInfo,
790 : DTEDMetaDataCode eCode,
791 : char **ppszLocation, int *pnLength )
792 : {
793 1708 : int bIsWeirdDTED = psDInfo->pachUHLRecord[4] == ' ';
794 :
795 1708 : switch( eCode )
796 : {
797 : case DTEDMD_ORIGINLONG:
798 72 : if (bIsWeirdDTED)
799 0 : *ppszLocation = psDInfo->pachUHLRecord + 8;
800 : else
801 72 : *ppszLocation = psDInfo->pachUHLRecord + 4;
802 72 : *pnLength = 8;
803 72 : break;
804 :
805 : case DTEDMD_ORIGINLAT:
806 72 : if (bIsWeirdDTED)
807 0 : *ppszLocation = psDInfo->pachUHLRecord + 24;
808 : else
809 72 : *ppszLocation = psDInfo->pachUHLRecord + 12;
810 72 : *pnLength = 8;
811 72 : break;
812 :
813 : case DTEDMD_VERTACCURACY_UHL:
814 74 : if (bIsWeirdDTED)
815 0 : *ppszLocation = psDInfo->pachUHLRecord + 56;
816 : else
817 74 : *ppszLocation = psDInfo->pachUHLRecord + 28;
818 74 : *pnLength = 4;
819 74 : break;
820 :
821 : case DTEDMD_SECURITYCODE_UHL:
822 74 : if (bIsWeirdDTED)
823 0 : *ppszLocation = psDInfo->pachUHLRecord + 60;
824 : else
825 74 : *ppszLocation = psDInfo->pachUHLRecord + 32;
826 74 : *pnLength = 3;
827 74 : break;
828 :
829 : case DTEDMD_UNIQUEREF_UHL:
830 74 : if (bIsWeirdDTED)
831 0 : *ppszLocation = NULL;
832 : else
833 74 : *ppszLocation = psDInfo->pachUHLRecord + 35;
834 74 : *pnLength = 12;
835 74 : break;
836 :
837 : case DTEDMD_DATA_EDITION:
838 74 : if (bIsWeirdDTED)
839 0 : *ppszLocation = psDInfo->pachDSIRecord + 174;
840 : else
841 74 : *ppszLocation = psDInfo->pachDSIRecord + 87;
842 74 : *pnLength = 2;
843 74 : break;
844 :
845 : case DTEDMD_MATCHMERGE_VERSION:
846 74 : if (bIsWeirdDTED)
847 0 : *ppszLocation = psDInfo->pachDSIRecord + 176;
848 : else
849 74 : *ppszLocation = psDInfo->pachDSIRecord + 89;
850 74 : *pnLength = 1;
851 74 : break;
852 :
853 : case DTEDMD_MAINT_DATE:
854 74 : if (bIsWeirdDTED)
855 0 : *ppszLocation = psDInfo->pachDSIRecord + 177;
856 : else
857 74 : *ppszLocation = psDInfo->pachDSIRecord + 90;
858 74 : *pnLength = 4;
859 74 : break;
860 :
861 : case DTEDMD_MATCHMERGE_DATE:
862 74 : if (bIsWeirdDTED)
863 0 : *ppszLocation = psDInfo->pachDSIRecord + 181;
864 : else
865 74 : *ppszLocation = psDInfo->pachDSIRecord + 94;
866 74 : *pnLength = 4;
867 74 : break;
868 :
869 : case DTEDMD_MAINT_DESCRIPTION:
870 74 : if (bIsWeirdDTED)
871 0 : *ppszLocation = psDInfo->pachDSIRecord + 185;
872 : else
873 74 : *ppszLocation = psDInfo->pachDSIRecord + 98;
874 74 : *pnLength = 4;
875 74 : break;
876 :
877 : case DTEDMD_PRODUCER:
878 74 : if (bIsWeirdDTED)
879 0 : *ppszLocation = psDInfo->pachDSIRecord + 189;
880 : else
881 74 : *ppszLocation = psDInfo->pachDSIRecord + 102;
882 74 : *pnLength = 8;
883 74 : break;
884 :
885 : case DTEDMD_VERTDATUM:
886 74 : if (bIsWeirdDTED)
887 0 : *ppszLocation = psDInfo->pachDSIRecord + 267;
888 : else
889 74 : *ppszLocation = psDInfo->pachDSIRecord + 141;
890 74 : *pnLength = 3;
891 74 : break;
892 :
893 : case DTEDMD_HORIZDATUM:
894 74 : if (bIsWeirdDTED)
895 0 : *ppszLocation = psDInfo->pachDSIRecord + 270;
896 : else
897 74 : *ppszLocation = psDInfo->pachDSIRecord + 144;
898 74 : *pnLength = 5;
899 74 : break;
900 :
901 : case DTEDMD_DIGITIZING_SYS:
902 74 : if (bIsWeirdDTED)
903 0 : *ppszLocation = NULL;
904 : else
905 74 : *ppszLocation = psDInfo->pachDSIRecord + 149;
906 74 : *pnLength = 10;
907 74 : break;
908 :
909 : case DTEDMD_COMPILATION_DATE:
910 74 : if (bIsWeirdDTED)
911 0 : *ppszLocation = NULL;
912 : else
913 74 : *ppszLocation = psDInfo->pachDSIRecord + 159;
914 74 : *pnLength = 4;
915 74 : break;
916 :
917 : case DTEDMD_HORIZACCURACY:
918 74 : *ppszLocation = psDInfo->pachACCRecord + 3;
919 74 : *pnLength = 4;
920 74 : break;
921 :
922 : case DTEDMD_REL_HORIZACCURACY:
923 74 : *ppszLocation = psDInfo->pachACCRecord + 11;
924 74 : *pnLength = 4;
925 74 : break;
926 :
927 : case DTEDMD_REL_VERTACCURACY:
928 74 : *ppszLocation = psDInfo->pachACCRecord + 15;
929 74 : *pnLength = 4;
930 74 : break;
931 :
932 : case DTEDMD_VERTACCURACY_ACC:
933 74 : *ppszLocation = psDInfo->pachACCRecord + 7;
934 74 : *pnLength = 4;
935 74 : break;
936 :
937 : case DTEDMD_SECURITYCODE_DSI:
938 74 : *ppszLocation = psDInfo->pachDSIRecord + 3;
939 74 : *pnLength = 1;
940 74 : break;
941 :
942 : case DTEDMD_UNIQUEREF_DSI:
943 74 : if (bIsWeirdDTED)
944 0 : *ppszLocation = NULL;
945 : else
946 74 : *ppszLocation = psDInfo->pachDSIRecord + 64;
947 74 : *pnLength = 15;
948 74 : break;
949 :
950 : case DTEDMD_NIMA_DESIGNATOR:
951 72 : if (bIsWeirdDTED)
952 0 : *ppszLocation = psDInfo->pachDSIRecord + 118;
953 : else
954 72 : *ppszLocation = psDInfo->pachDSIRecord + 59;
955 72 : *pnLength = 5;
956 72 : break;
957 :
958 : case DTEDMD_PARTIALCELL_DSI:
959 86 : if (bIsWeirdDTED)
960 0 : *ppszLocation = NULL;
961 : else
962 86 : *ppszLocation = psDInfo->pachDSIRecord + 289;
963 86 : *pnLength = 2;
964 86 : break;
965 :
966 : default:
967 0 : *ppszLocation = NULL;
968 0 : *pnLength = 0;
969 : }
970 1708 : }
971 :
972 : /************************************************************************/
973 : /* DTEDGetMetadata() */
974 : /************************************************************************/
975 :
976 1656 : char *DTEDGetMetadata( DTEDInfo *psDInfo, DTEDMetaDataCode eCode )
977 :
978 : {
979 : int nFieldLen;
980 : char *pszFieldSrc;
981 : char *pszResult;
982 :
983 1656 : DTEDGetMetadataLocation( psDInfo, eCode, &pszFieldSrc, &nFieldLen );
984 1656 : if( pszFieldSrc == NULL )
985 0 : return strdup( "" );
986 :
987 1656 : pszResult = (char *) malloc(nFieldLen+1);
988 1656 : strncpy( pszResult, pszFieldSrc, nFieldLen );
989 1656 : pszResult[nFieldLen] = '\0';
990 :
991 1656 : return pszResult;
992 : }
993 :
994 : /************************************************************************/
995 : /* DTEDSetMetadata() */
996 : /************************************************************************/
997 :
998 52 : int DTEDSetMetadata( DTEDInfo *psDInfo, DTEDMetaDataCode eCode,
999 : const char *pszNewValue )
1000 :
1001 : {
1002 : int nFieldLen;
1003 : char *pszFieldSrc;
1004 :
1005 52 : if( !psDInfo->bUpdate )
1006 0 : return FALSE;
1007 :
1008 : /* -------------------------------------------------------------------- */
1009 : /* Get the location in the headers to update. */
1010 : /* -------------------------------------------------------------------- */
1011 52 : DTEDGetMetadataLocation( psDInfo, eCode, &pszFieldSrc, &nFieldLen );
1012 52 : if( pszFieldSrc == NULL )
1013 0 : return FALSE;
1014 :
1015 : /* -------------------------------------------------------------------- */
1016 : /* Update it, padding with spaces. */
1017 : /* -------------------------------------------------------------------- */
1018 52 : memset( pszFieldSrc, ' ', nFieldLen );
1019 52 : strncpy( pszFieldSrc, pszNewValue,
1020 52 : MIN((size_t)nFieldLen,strlen(pszNewValue)) );
1021 :
1022 : /* Turn the flag on, so that the headers are rewritten at file */
1023 : /* closing */
1024 52 : psDInfo->bRewriteHeaders = TRUE;
1025 :
1026 52 : return TRUE;
1027 : }
1028 :
1029 : /************************************************************************/
1030 : /* DTEDClose() */
1031 : /************************************************************************/
1032 :
1033 86 : void DTEDClose( DTEDInfo * psDInfo )
1034 :
1035 : {
1036 86 : if( psDInfo->bRewriteHeaders )
1037 : {
1038 : /* -------------------------------------------------------------------- */
1039 : /* Write all headers back to disk. */
1040 : /* -------------------------------------------------------------------- */
1041 14 : VSIFSeekL( psDInfo->fp, psDInfo->nUHLOffset, SEEK_SET );
1042 14 : VSIFWriteL( psDInfo->pachUHLRecord, 1, DTED_UHL_SIZE, psDInfo->fp );
1043 :
1044 14 : VSIFSeekL( psDInfo->fp, psDInfo->nDSIOffset, SEEK_SET );
1045 14 : VSIFWriteL( psDInfo->pachDSIRecord, 1, DTED_DSI_SIZE, psDInfo->fp );
1046 :
1047 14 : VSIFSeekL( psDInfo->fp, psDInfo->nACCOffset, SEEK_SET );
1048 14 : VSIFWriteL( psDInfo->pachACCRecord, 1, DTED_ACC_SIZE, psDInfo->fp );
1049 : }
1050 :
1051 86 : VSIFCloseL( psDInfo->fp );
1052 :
1053 86 : CPLFree( psDInfo->pachUHLRecord );
1054 86 : CPLFree( psDInfo->pachDSIRecord );
1055 86 : CPLFree( psDInfo->pachACCRecord );
1056 :
1057 86 : CPLFree( psDInfo->panMapLogicalColsToOffsets );
1058 :
1059 86 : CPLFree( psDInfo );
1060 86 : }
|
