1 : /******************************************************************************
2 : * $Id: dteddataset.cpp 18790 2010-02-11 19:05:02Z warmerdam $
3 : *
4 : * Project: DTED Translator
5 : * Purpose: GDALDataset driver for DTED translator.
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 : #include "gdal_pam.h"
32 : #include "ogr_spatialref.h"
33 :
34 : CPL_CVSID("$Id: dteddataset.cpp 18790 2010-02-11 19:05:02Z warmerdam $");
35 :
36 : CPL_C_START
37 : void GDALRegister_DTED(void);
38 : CPL_C_END
39 :
40 : /************************************************************************/
41 : /* ==================================================================== */
42 : /* DTEDDataset */
43 : /* ==================================================================== */
44 : /************************************************************************/
45 :
46 : class DTEDRasterBand;
47 :
48 : class DTEDDataset : public GDALPamDataset
49 : {
50 : friend class DTEDRasterBand;
51 :
52 : char *pszFilename;
53 : DTEDInfo *psDTED;
54 : int bVerifyChecksum;
55 : char *pszProjection;
56 :
57 : public:
58 : DTEDDataset();
59 : virtual ~DTEDDataset();
60 :
61 : virtual const char *GetProjectionRef(void);
62 : virtual CPLErr GetGeoTransform( double * );
63 :
64 1 : const char* GetFileName() { return pszFilename; }
65 : void SetFileName(const char* pszFilename);
66 :
67 : static GDALDataset *Open( GDALOpenInfo * );
68 : };
69 :
70 : /************************************************************************/
71 : /* ==================================================================== */
72 : /* DTEDRasterBand */
73 : /* ==================================================================== */
74 : /************************************************************************/
75 :
76 : class DTEDRasterBand : public GDALPamRasterBand
77 66 : {
78 : friend class DTEDDataset;
79 :
80 : int bNoDataSet;
81 : double dfNoDataValue;
82 :
83 : public:
84 :
85 : DTEDRasterBand( DTEDDataset *, int );
86 :
87 : virtual CPLErr IReadBlock( int, int, void * );
88 : virtual CPLErr IWriteBlock( int, int, void * );
89 :
90 : virtual double GetNoDataValue( int *pbSuccess = NULL );
91 :
92 23 : virtual const char* GetUnitType() { return "m"; }
93 : };
94 :
95 :
96 : /************************************************************************/
97 : /* DTEDRasterBand() */
98 : /************************************************************************/
99 :
100 66 : DTEDRasterBand::DTEDRasterBand( DTEDDataset *poDS, int nBand )
101 :
102 : {
103 66 : this->poDS = poDS;
104 66 : this->nBand = nBand;
105 :
106 66 : eDataType = GDT_Int16;
107 :
108 66 : bNoDataSet = TRUE;
109 66 : dfNoDataValue = (double) DTED_NODATA_VALUE;
110 :
111 : /* For some applications, it may be valuable to consider the whole DTED */
112 : /* file as single block, as the column orientation doesn't fit very well */
113 : /* with some scanline oriented algorithms */
114 : /* Of course you need to have a big enough case size, particularly for DTED 2 */
115 : /* datasets */
116 : nBlockXSize = CSLTestBoolean(CPLGetConfigOption("GDAL_DTED_SINGLE_BLOCK", "NO")) ?
117 66 : poDS->GetRasterXSize() : 1;
118 66 : nBlockYSize = poDS->GetRasterYSize();
119 66 : }
120 :
121 : /************************************************************************/
122 : /* IReadBlock() */
123 : /************************************************************************/
124 :
125 : CPLErr DTEDRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
126 4510 : void * pImage )
127 :
128 : {
129 4510 : DTEDDataset *poDTED_DS = (DTEDDataset *) poDS;
130 4510 : int nYSize = poDTED_DS->psDTED->nYSize;
131 : GInt16 *panData;
132 :
133 : (void) nBlockXOff;
134 4510 : CPLAssert( nBlockYOff == 0 );
135 :
136 4510 : if (nBlockXSize != 1)
137 : {
138 0 : panData = (GInt16 *) pImage;
139 0 : GInt16* panBuffer = (GInt16*) CPLMalloc(sizeof(GInt16) * nBlockYSize);
140 : int i;
141 0 : for(i=0;i<nBlockXSize;i++)
142 : {
143 0 : if( !DTEDReadProfileEx( poDTED_DS->psDTED, i, panBuffer,
144 : poDTED_DS->bVerifyChecksum ) )
145 : {
146 0 : CPLFree(panBuffer);
147 0 : return CE_Failure;
148 : }
149 : int j;
150 0 : for(j=0;j<nBlockYSize;j++)
151 : {
152 0 : panData[j * nBlockXSize + i] = panBuffer[nYSize - j - 1];
153 : }
154 : }
155 :
156 0 : CPLFree(panBuffer);
157 0 : return CE_None;
158 : }
159 :
160 : /* -------------------------------------------------------------------- */
161 : /* Read the data. */
162 : /* -------------------------------------------------------------------- */
163 4510 : panData = (GInt16 *) pImage;
164 4510 : if( !DTEDReadProfileEx( poDTED_DS->psDTED, nBlockXOff, panData,
165 : poDTED_DS->bVerifyChecksum ) )
166 1 : return CE_Failure;
167 :
168 : /* -------------------------------------------------------------------- */
169 : /* Flip line to orient it top to bottom instead of bottom to */
170 : /* top. */
171 : /* -------------------------------------------------------------------- */
172 604098 : for( int i = nYSize/2; i >= 0; i-- )
173 : {
174 : GInt16 nTemp;
175 :
176 599589 : nTemp = panData[i];
177 599589 : panData[i] = panData[nYSize - i - 1];
178 599589 : panData[nYSize - i - 1] = nTemp;
179 : }
180 :
181 4509 : return CE_None;
182 : }
183 :
184 : /************************************************************************/
185 : /* IReadBlock() */
186 : /************************************************************************/
187 :
188 : CPLErr DTEDRasterBand::IWriteBlock( int nBlockXOff, int nBlockYOff,
189 0 : void * pImage )
190 :
191 : {
192 0 : DTEDDataset *poDTED_DS = (DTEDDataset *) poDS;
193 : GInt16 *panData;
194 :
195 : (void) nBlockXOff;
196 0 : CPLAssert( nBlockYOff == 0 );
197 :
198 0 : if (poDTED_DS->eAccess != GA_Update)
199 0 : return CE_Failure;
200 :
201 0 : if (nBlockXSize != 1)
202 : {
203 0 : panData = (GInt16 *) pImage;
204 0 : GInt16* panBuffer = (GInt16*) CPLMalloc(sizeof(GInt16) * nBlockYSize);
205 : int i;
206 0 : for(i=0;i<nBlockXSize;i++)
207 : {
208 : int j;
209 0 : for(j=0;j<nBlockYSize;j++)
210 : {
211 0 : panBuffer[j] = panData[j * nBlockXSize + i];
212 : }
213 0 : if( !DTEDWriteProfile( poDTED_DS->psDTED, i, panBuffer) )
214 : {
215 0 : CPLFree(panBuffer);
216 0 : return CE_Failure;
217 : }
218 : }
219 :
220 0 : CPLFree(panBuffer);
221 0 : return CE_None;
222 : }
223 :
224 0 : panData = (GInt16 *) pImage;
225 0 : if( !DTEDWriteProfile( poDTED_DS->psDTED, nBlockXOff, panData) )
226 0 : return CE_Failure;
227 :
228 0 : return CE_None;
229 : }
230 :
231 : /************************************************************************/
232 : /* GetNoDataValue() */
233 : /************************************************************************/
234 :
235 61 : double DTEDRasterBand::GetNoDataValue( int * pbSuccess )
236 :
237 : {
238 61 : if( pbSuccess )
239 53 : *pbSuccess = bNoDataSet;
240 :
241 61 : return dfNoDataValue;
242 : }
243 :
244 : /************************************************************************/
245 : /* ~DTEDDataset() */
246 : /************************************************************************/
247 :
248 66 : DTEDDataset::DTEDDataset()
249 : {
250 66 : pszFilename = CPLStrdup("unknown");
251 66 : pszProjection = CPLStrdup("");
252 66 : bVerifyChecksum = CSLTestBoolean(CPLGetConfigOption("DTED_VERIFY_CHECKSUM", "NO"));
253 66 : }
254 :
255 : /************************************************************************/
256 : /* ~DTEDDataset() */
257 : /************************************************************************/
258 :
259 66 : DTEDDataset::~DTEDDataset()
260 :
261 : {
262 66 : FlushCache();
263 66 : CPLFree(pszFilename);
264 66 : CPLFree( pszProjection );
265 66 : if( psDTED != NULL )
266 66 : DTEDClose( psDTED );
267 66 : }
268 :
269 : /************************************************************************/
270 : /* SetFileName() */
271 : /************************************************************************/
272 :
273 66 : void DTEDDataset::SetFileName(const char* pszFilename)
274 :
275 : {
276 66 : CPLFree(this->pszFilename);
277 66 : this->pszFilename = CPLStrdup(pszFilename);
278 66 : }
279 :
280 : /************************************************************************/
281 : /* Open() */
282 : /************************************************************************/
283 :
284 11272 : GDALDataset *DTEDDataset::Open( GDALOpenInfo * poOpenInfo )
285 :
286 : {
287 : int i;
288 : DTEDInfo *psDTED;
289 :
290 : /* -------------------------------------------------------------------- */
291 : /* Does the file start with one of the possible DTED header */
292 : /* record types, and do we have a UHL marker? */
293 : /* -------------------------------------------------------------------- */
294 11272 : if( poOpenInfo->nHeaderBytes < 240 )
295 9963 : return NULL;
296 :
297 1309 : if( !EQUALN((const char *)poOpenInfo->pabyHeader,"VOL",3)
298 : && !EQUALN((const char *)poOpenInfo->pabyHeader,"HDR",3)
299 : && !EQUALN((const char *)poOpenInfo->pabyHeader,"UHL",3) )
300 : {
301 1243 : return NULL;
302 : }
303 :
304 66 : int bFoundUHL = FALSE;
305 133 : for(i=0;i<poOpenInfo->nHeaderBytes-3 && !bFoundUHL ;i += DTED_UHL_SIZE)
306 : {
307 67 : if( EQUALN((const char *)poOpenInfo->pabyHeader + i,"UHL", 3) )
308 : {
309 66 : bFoundUHL = TRUE;
310 : }
311 : }
312 66 : if (!bFoundUHL)
313 0 : return NULL;
314 :
315 : /* -------------------------------------------------------------------- */
316 : /* Try opening the dataset. */
317 : /* -------------------------------------------------------------------- */
318 66 : psDTED = DTEDOpen( poOpenInfo->pszFilename, (poOpenInfo->eAccess == GA_Update) ? "rb+" : "rb", TRUE );
319 :
320 66 : if( psDTED == NULL )
321 0 : return( NULL );
322 :
323 : /* -------------------------------------------------------------------- */
324 : /* Create a corresponding GDALDataset. */
325 : /* -------------------------------------------------------------------- */
326 : DTEDDataset *poDS;
327 :
328 66 : poDS = new DTEDDataset();
329 66 : poDS->SetFileName(poOpenInfo->pszFilename);
330 :
331 66 : poDS->eAccess = poOpenInfo->eAccess;
332 66 : poDS->psDTED = psDTED;
333 :
334 : /* -------------------------------------------------------------------- */
335 : /* Capture some information from the file that is of interest. */
336 : /* -------------------------------------------------------------------- */
337 66 : poDS->nRasterXSize = psDTED->nXSize;
338 66 : poDS->nRasterYSize = psDTED->nYSize;
339 :
340 66 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
341 : {
342 0 : delete poDS;
343 0 : return NULL;
344 : }
345 :
346 : /* -------------------------------------------------------------------- */
347 : /* Create band information objects. */
348 : /* -------------------------------------------------------------------- */
349 66 : poDS->nBands = 1;;
350 264 : for( i = 0; i < poDS->nBands; i++ )
351 66 : poDS->SetBand( i+1, new DTEDRasterBand( poDS, i+1 ) );
352 :
353 : /* -------------------------------------------------------------------- */
354 : /* Collect any metadata available. */
355 : /* -------------------------------------------------------------------- */
356 : char *pszValue;
357 :
358 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_VERTACCURACY_UHL );
359 66 : poDS->SetMetadataItem( "DTED_VerticalAccuracy_UHL", pszValue );
360 66 : free( pszValue );
361 :
362 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_VERTACCURACY_ACC );
363 66 : poDS->SetMetadataItem( "DTED_VerticalAccuracy_ACC", pszValue );
364 66 : free( pszValue );
365 :
366 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_SECURITYCODE_UHL );
367 66 : poDS->SetMetadataItem( "DTED_SecurityCode_UHL", pszValue );
368 66 : free( pszValue );
369 :
370 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_SECURITYCODE_DSI );
371 66 : poDS->SetMetadataItem( "DTED_SecurityCode_DSI", pszValue );
372 66 : free( pszValue );
373 :
374 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_UNIQUEREF_UHL );
375 66 : poDS->SetMetadataItem( "DTED_UniqueRef_UHL", pszValue );
376 66 : free( pszValue );
377 :
378 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_UNIQUEREF_DSI );
379 66 : poDS->SetMetadataItem( "DTED_UniqueRef_DSI", pszValue );
380 66 : free( pszValue );
381 :
382 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_DATA_EDITION );
383 66 : poDS->SetMetadataItem( "DTED_DataEdition", pszValue );
384 66 : free( pszValue );
385 :
386 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_MATCHMERGE_VERSION );
387 66 : poDS->SetMetadataItem( "DTED_MatchMergeVersion", pszValue );
388 66 : free( pszValue );
389 :
390 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_MAINT_DATE );
391 66 : poDS->SetMetadataItem( "DTED_MaintenanceDate", pszValue );
392 66 : free( pszValue );
393 :
394 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_MATCHMERGE_DATE );
395 66 : poDS->SetMetadataItem( "DTED_MatchMergeDate", pszValue );
396 66 : free( pszValue );
397 :
398 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_MAINT_DESCRIPTION );
399 66 : poDS->SetMetadataItem( "DTED_MaintenanceDescription", pszValue );
400 66 : free( pszValue );
401 :
402 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_PRODUCER );
403 66 : poDS->SetMetadataItem( "DTED_Producer", pszValue );
404 66 : free( pszValue );
405 :
406 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_VERTDATUM );
407 66 : poDS->SetMetadataItem( "DTED_VerticalDatum", pszValue );
408 66 : free( pszValue );
409 :
410 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_HORIZDATUM );
411 66 : poDS->SetMetadataItem( "DTED_HorizontalDatum", pszValue );
412 66 : free( pszValue );
413 :
414 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_DIGITIZING_SYS );
415 66 : poDS->SetMetadataItem( "DTED_DigitizingSystem", pszValue );
416 66 : free( pszValue );
417 :
418 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_COMPILATION_DATE );
419 66 : poDS->SetMetadataItem( "DTED_CompilationDate", pszValue );
420 66 : free( pszValue );
421 :
422 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_HORIZACCURACY );
423 66 : poDS->SetMetadataItem( "DTED_HorizontalAccuracy", pszValue );
424 66 : free( pszValue );
425 :
426 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_REL_HORIZACCURACY );
427 66 : poDS->SetMetadataItem( "DTED_RelHorizontalAccuracy", pszValue );
428 66 : free( pszValue );
429 :
430 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_REL_VERTACCURACY );
431 66 : poDS->SetMetadataItem( "DTED_RelVerticalAccuracy", pszValue );
432 66 : free( pszValue );
433 :
434 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_ORIGINLAT );
435 66 : poDS->SetMetadataItem( "DTED_OriginLatitude", pszValue );
436 66 : free( pszValue );
437 :
438 66 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_ORIGINLONG );
439 66 : poDS->SetMetadataItem( "DTED_OriginLongitude", pszValue );
440 66 : free( pszValue );
441 :
442 66 : poDS->SetMetadataItem( GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT );
443 :
444 : /* -------------------------------------------------------------------- */
445 : /* Initialize any PAM information. */
446 : /* -------------------------------------------------------------------- */
447 66 : poDS->SetDescription( poOpenInfo->pszFilename );
448 66 : poDS->TryLoadXML();
449 :
450 : // if no SR in xml, try aux
451 66 : const char* pszPrj = poDS->GDALPamDataset::GetProjectionRef();
452 66 : if( !pszPrj || strlen(pszPrj) == 0 )
453 : {
454 66 : GDALDataset* poAuxDS = GDALFindAssociatedAuxFile( poOpenInfo->pszFilename, GA_ReadOnly, poDS );
455 66 : if( poAuxDS )
456 : {
457 0 : pszPrj = poAuxDS->GetProjectionRef();
458 0 : if( pszPrj && strlen(pszPrj) > 0 )
459 : {
460 0 : CPLFree( poDS->pszProjection );
461 0 : poDS->pszProjection = CPLStrdup(pszPrj);
462 : }
463 :
464 0 : GDALClose( poAuxDS );
465 : }
466 : }
467 :
468 : /* -------------------------------------------------------------------- */
469 : /* Support overviews. */
470 : /* -------------------------------------------------------------------- */
471 66 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
472 66 : return( poDS );
473 : }
474 :
475 : /************************************************************************/
476 : /* GetGeoTransform() */
477 : /************************************************************************/
478 :
479 72 : CPLErr DTEDDataset::GetGeoTransform( double * padfTransform )
480 :
481 : {
482 72 : padfTransform[0] = psDTED->dfULCornerX;
483 72 : padfTransform[1] = psDTED->dfPixelSizeX;
484 72 : padfTransform[2] = 0.0;
485 72 : padfTransform[3] = psDTED->dfULCornerY;
486 72 : padfTransform[4] = 0.0;
487 72 : padfTransform[5] = psDTED->dfPixelSizeY * -1;
488 :
489 72 : return( CE_None );
490 : }
491 :
492 : /************************************************************************/
493 : /* GetProjectionRef() */
494 : /************************************************************************/
495 :
496 69 : const char *DTEDDataset::GetProjectionRef()
497 :
498 : {
499 : // get xml and aux SR first
500 69 : const char* pszPrj = GDALPamDataset::GetProjectionRef();
501 69 : if(pszPrj && strlen(pszPrj) > 0)
502 0 : return pszPrj;
503 :
504 69 : if (pszProjection && strlen(pszProjection) > 0)
505 0 : return pszProjection;
506 :
507 69 : pszPrj = GetMetadataItem( "DTED_HorizontalDatum");
508 69 : if (EQUAL(pszPrj, "WGS84"))
509 : {
510 68 : return( "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]" );
511 : }
512 1 : else if (EQUAL(pszPrj, "WGS72"))
513 : {
514 : static int bWarned = FALSE;
515 1 : if (!bWarned)
516 : {
517 1 : bWarned = TRUE;
518 : CPLError( CE_Warning, CPLE_AppDefined,
519 : "The DTED file %s indicates WGS72 as horizontal datum. \n"
520 : "As this is outdated nowadays, you should contact your data producer to get data georeferenced in WGS84.\n"
521 : "In some cases, WGS72 is a wrong indication and the georeferencing is really WGS84. In that case\n"
522 : "you might consider doing 'gdal_translate -of DTED -mo \"DTED_HorizontalDatum=WGS84\" src.dtX dst.dtX' to\n"
523 : "fix the DTED file.\n"
524 1 : "No more warnings will be issued in this session about this operation.", GetFileName() );
525 : }
526 1 : return "GEOGCS[\"WGS 72\",DATUM[\"WGS_1972\",SPHEROID[\"WGS 72\",6378135,298.26]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4322\"]]";
527 : }
528 : else
529 : {
530 : static int bWarned = FALSE;
531 0 : if (!bWarned)
532 : {
533 0 : bWarned = TRUE;
534 : CPLError( CE_Warning, CPLE_AppDefined,
535 : "The DTED file %s indicates %s as horizontal datum, which is not recognized by the DTED driver. \n"
536 : "The DTED driver is going to consider it as WGS84.\n"
537 0 : "No more warnings will be issued in this session about this operation.", GetFileName(), pszProjection );
538 : }
539 0 : return( "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]" );
540 : }
541 : }
542 :
543 : /************************************************************************/
544 : /* DTEDCreateCopy() */
545 : /* */
546 : /* For now we will assume the input is exactly one proper */
547 : /* cell. */
548 : /************************************************************************/
549 :
550 : static GDALDataset *
551 : DTEDCreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
552 : int bStrict, char ** papszOptions,
553 19 : GDALProgressFunc pfnProgress, void * pProgressData )
554 :
555 : {
556 : (void) pProgressData;
557 : (void) pfnProgress;
558 : (void) papszOptions;
559 : (void) bStrict;
560 :
561 :
562 : /* -------------------------------------------------------------------- */
563 : /* Some some rudimentary checks */
564 : /* -------------------------------------------------------------------- */
565 19 : int nBands = poSrcDS->GetRasterCount();
566 19 : if (nBands == 0)
567 : {
568 : CPLError( CE_Failure, CPLE_NotSupported,
569 1 : "DTED driver does not support source dataset with zero band.\n");
570 1 : return NULL;
571 : }
572 :
573 18 : if (nBands != 1)
574 : {
575 : CPLError( (bStrict) ? CE_Failure : CE_Warning, CPLE_NotSupported,
576 4 : "DTED driver only uses the first band of the dataset.\n");
577 4 : if (bStrict)
578 4 : return NULL;
579 : }
580 :
581 14 : if( pfnProgress && !pfnProgress( 0.0, NULL, pProgressData ) )
582 0 : return NULL;
583 :
584 : /* -------------------------------------------------------------------- */
585 : /* Work out the level. */
586 : /* -------------------------------------------------------------------- */
587 : int nLevel;
588 :
589 14 : if( poSrcDS->GetRasterYSize() == 121 )
590 2 : nLevel = 0;
591 12 : else if( poSrcDS->GetRasterYSize() == 1201 )
592 1 : nLevel = 1;
593 11 : else if( poSrcDS->GetRasterYSize() == 3601 )
594 0 : nLevel = 2;
595 : else
596 : {
597 : CPLError( CE_Warning, CPLE_AppDefined,
598 11 : "The source does not appear to be a properly formatted cell." );
599 11 : nLevel = 1;
600 : }
601 :
602 : /* -------------------------------------------------------------------- */
603 : /* Checks the input SRS */
604 : /* -------------------------------------------------------------------- */
605 14 : OGRSpatialReference ogrsr_input;
606 14 : OGRSpatialReference ogrsr_wgs84;
607 14 : char* c = (char*)poSrcDS->GetProjectionRef();
608 14 : ogrsr_input.importFromWkt(&c);
609 14 : ogrsr_wgs84.SetWellKnownGeogCS( "WGS84" );
610 14 : if ( ogrsr_input.IsSameGeogCS(&ogrsr_wgs84) == FALSE)
611 : {
612 : CPLError( CE_Warning, CPLE_AppDefined,
613 : "The source projection coordinate system is %s. Only WGS 84 is supported.\n"
614 : "The DTED driver will generate a file as if the source was WGS 84 projection coordinate system.",
615 0 : poSrcDS->GetProjectionRef() );
616 : }
617 :
618 : /* -------------------------------------------------------------------- */
619 : /* Work out the LL origin. */
620 : /* -------------------------------------------------------------------- */
621 : int nLLOriginLat, nLLOriginLong;
622 : double adfGeoTransform[6];
623 :
624 14 : poSrcDS->GetGeoTransform( adfGeoTransform );
625 :
626 : nLLOriginLat = (int)
627 : floor(adfGeoTransform[3]
628 14 : + poSrcDS->GetRasterYSize() * adfGeoTransform[5] + 0.5);
629 :
630 14 : nLLOriginLong = (int) floor(adfGeoTransform[0] + 0.5);
631 :
632 14 : if (fabs(nLLOriginLat - (adfGeoTransform[3]
633 : + (poSrcDS->GetRasterYSize() - 0.5) * adfGeoTransform[5])) > 1e-10 ||
634 : fabs(nLLOriginLong - (adfGeoTransform[0] + 0.5 * adfGeoTransform[1])) > 1e-10)
635 : {
636 : CPLError( CE_Warning, CPLE_AppDefined,
637 : "The corner coordinates of the source are not properly "
638 11 : "aligned on plain latitude/longitude boundaries.");
639 : }
640 :
641 : /* -------------------------------------------------------------------- */
642 : /* Check horizontal source size. */
643 : /* -------------------------------------------------------------------- */
644 : int expectedXSize;
645 14 : if( ABS(nLLOriginLat) >= 80 )
646 0 : expectedXSize = (poSrcDS->GetRasterYSize() - 1) / 6 + 1;
647 14 : else if( ABS(nLLOriginLat) >= 75 )
648 0 : expectedXSize = (poSrcDS->GetRasterYSize() - 1) / 4 + 1;
649 14 : else if( ABS(nLLOriginLat) >= 70 )
650 0 : expectedXSize = (poSrcDS->GetRasterYSize() - 1) / 3 + 1;
651 14 : else if( ABS(nLLOriginLat) >= 50 )
652 1 : expectedXSize = (poSrcDS->GetRasterYSize() - 1) / 2 + 1;
653 : else
654 13 : expectedXSize = poSrcDS->GetRasterYSize();
655 :
656 14 : if (poSrcDS->GetRasterXSize() != expectedXSize)
657 : {
658 : CPLError( CE_Warning, CPLE_AppDefined,
659 : "The horizontal source size is not conformant with the one "
660 : "expected by DTED Level %d at this latitude (%d pixels found instead of %d).", nLevel,
661 0 : poSrcDS->GetRasterXSize(), expectedXSize);
662 : }
663 :
664 : /* -------------------------------------------------------------------- */
665 : /* Create the output dted file. */
666 : /* -------------------------------------------------------------------- */
667 : const char *pszError;
668 :
669 14 : pszError = DTEDCreate( pszFilename, nLevel, nLLOriginLat, nLLOriginLong );
670 :
671 14 : if( pszError != NULL )
672 : {
673 : CPLError( CE_Failure, CPLE_AppDefined,
674 0 : "%s", pszError );
675 14 : return NULL;
676 : }
677 :
678 : /* -------------------------------------------------------------------- */
679 : /* Open the DTED file so we can output the data to it. */
680 : /* -------------------------------------------------------------------- */
681 : DTEDInfo *psDTED;
682 :
683 14 : psDTED = DTEDOpen( pszFilename, "rb+", FALSE );
684 14 : if( psDTED == NULL )
685 0 : return NULL;
686 :
687 : /* -------------------------------------------------------------------- */
688 : /* Read all the data in a single buffer. */
689 : /* -------------------------------------------------------------------- */
690 14 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( 1 );
691 : GInt16 *panData;
692 :
693 : panData = (GInt16 *)
694 14 : CPLMalloc(sizeof(GInt16) * psDTED->nXSize * psDTED->nYSize);
695 :
696 14668 : for( int iY = 0; iY < psDTED->nYSize; iY++ )
697 : {
698 : poSrcBand->RasterIO( GF_Read, 0, iY, psDTED->nXSize, 1,
699 : (void *) (panData + iY * psDTED->nXSize), psDTED->nXSize, 1,
700 14654 : GDT_Int16, 0, 0 );
701 :
702 14654 : if( pfnProgress && !pfnProgress(0.5 * (iY+1) / (double) psDTED->nYSize, NULL, pProgressData ) )
703 : {
704 : CPLError( CE_Failure, CPLE_UserInterrupt,
705 0 : "User terminated CreateCopy()" );
706 0 : DTEDClose( psDTED );
707 0 : CPLFree( panData );
708 0 : return NULL;
709 : }
710 : }
711 :
712 : int bSrcBandHasNoData;
713 14 : double srcBandNoData = poSrcBand->GetNoDataValue(&bSrcBandHasNoData);
714 :
715 : /* -------------------------------------------------------------------- */
716 : /* Write all the profiles. */
717 : /* -------------------------------------------------------------------- */
718 : GInt16 anProfData[3601];
719 14 : int dfNodataCount=0;
720 : GByte iPartialCell;
721 :
722 14068 : for( int iProfile = 0; iProfile < psDTED->nXSize; iProfile++ )
723 : {
724 16631548 : for( int iY = 0; iY < psDTED->nYSize; iY++ )
725 : {
726 16617494 : anProfData[iY] = panData[iProfile + iY * psDTED->nXSize];
727 16617494 : if ( bSrcBandHasNoData && anProfData[iY] == srcBandNoData)
728 : {
729 0 : anProfData[iY] = DTED_NODATA_VALUE;
730 0 : dfNodataCount++;
731 : }
732 16617494 : else if ( anProfData[iY] == DTED_NODATA_VALUE )
733 0 : dfNodataCount++;
734 : }
735 14054 : DTEDWriteProfile( psDTED, iProfile, anProfData );
736 :
737 14054 : if( pfnProgress
738 : && !pfnProgress( 0.5 + 0.5 * (iProfile+1) / (double) psDTED->nXSize,
739 : NULL, pProgressData ) )
740 : {
741 : CPLError( CE_Failure, CPLE_UserInterrupt,
742 0 : "User terminated CreateCopy()" );
743 0 : DTEDClose( psDTED );
744 0 : CPLFree( panData );
745 0 : return NULL;
746 : }
747 : }
748 14 : CPLFree( panData );
749 :
750 : /* -------------------------------------------------------------------- */
751 : /* Partial cell indicator: 0 for complete coverage; 1-99 for incomplete */
752 : /* -------------------------------------------------------------------- */
753 : char szPartialCell[3];
754 :
755 14 : if ( dfNodataCount == 0 )
756 14 : iPartialCell = 0;
757 : else
758 : {
759 : iPartialCell = (GByte)int(floor(100.0 -
760 0 : (dfNodataCount*100.0/(psDTED->nXSize * psDTED->nYSize))));
761 0 : if (iPartialCell < 1)
762 0 : iPartialCell=1;
763 : }
764 14 : sprintf(szPartialCell,"%02d",iPartialCell);
765 14 : strncpy((char *) (psDTED->pachDSIRecord+289), szPartialCell, 2 );
766 :
767 : /* -------------------------------------------------------------------- */
768 : /* Try to copy any matching available metadata. */
769 : /* -------------------------------------------------------------------- */
770 14 : if( poSrcDS->GetMetadataItem( "DTED_VerticalAccuracy_UHL" ) != NULL )
771 : DTEDSetMetadata( psDTED, DTEDMD_VERTACCURACY_UHL,
772 2 : poSrcDS->GetMetadataItem( "DTED_VerticalAccuracy_UHL" ) );
773 :
774 14 : if( poSrcDS->GetMetadataItem( "DTED_VerticalAccuracy_ACC" ) != NULL )
775 : DTEDSetMetadata( psDTED, DTEDMD_VERTACCURACY_ACC,
776 2 : poSrcDS->GetMetadataItem( "DTED_VerticalAccuracy_ACC" ) );
777 :
778 14 : if( poSrcDS->GetMetadataItem( "DTED_SecurityCode_UHL" ) != NULL )
779 : DTEDSetMetadata( psDTED, DTEDMD_SECURITYCODE_UHL,
780 2 : poSrcDS->GetMetadataItem( "DTED_SecurityCode_UHL" ) );
781 :
782 14 : if( poSrcDS->GetMetadataItem( "DTED_SecurityCode_DSI" ) != NULL )
783 : DTEDSetMetadata( psDTED, DTEDMD_SECURITYCODE_DSI,
784 2 : poSrcDS->GetMetadataItem( "DTED_SecurityCode_DSI" ) );
785 :
786 14 : if( poSrcDS->GetMetadataItem( "DTED_UniqueRef_UHL" ) != NULL )
787 : DTEDSetMetadata( psDTED, DTEDMD_UNIQUEREF_UHL,
788 2 : poSrcDS->GetMetadataItem( "DTED_UniqueRef_UHL" ) );
789 :
790 14 : if( poSrcDS->GetMetadataItem( "DTED_UniqueRef_DSI" ) != NULL )
791 : DTEDSetMetadata( psDTED, DTEDMD_UNIQUEREF_DSI,
792 2 : poSrcDS->GetMetadataItem( "DTED_UniqueRef_DSI" ) );
793 :
794 14 : if( poSrcDS->GetMetadataItem( "DTED_DataEdition" ) != NULL )
795 : DTEDSetMetadata( psDTED, DTEDMD_DATA_EDITION,
796 2 : poSrcDS->GetMetadataItem( "DTED_DataEdition" ) );
797 :
798 14 : if( poSrcDS->GetMetadataItem( "DTED_MatchMergeVersion" ) != NULL )
799 : DTEDSetMetadata( psDTED, DTEDMD_MATCHMERGE_VERSION,
800 2 : poSrcDS->GetMetadataItem( "DTED_MatchMergeVersion" ) );
801 :
802 14 : if( poSrcDS->GetMetadataItem( "DTED_MaintenanceDate" ) != NULL )
803 : DTEDSetMetadata( psDTED, DTEDMD_MAINT_DATE,
804 2 : poSrcDS->GetMetadataItem( "DTED_MaintenanceDate" ) );
805 :
806 14 : if( poSrcDS->GetMetadataItem( "DTED_MatchMergeDate" ) != NULL )
807 : DTEDSetMetadata( psDTED, DTEDMD_MATCHMERGE_DATE,
808 2 : poSrcDS->GetMetadataItem( "DTED_MatchMergeDate" ) );
809 :
810 14 : if( poSrcDS->GetMetadataItem( "DTED_MaintenanceDescription" ) != NULL )
811 : DTEDSetMetadata( psDTED, DTEDMD_MAINT_DESCRIPTION,
812 2 : poSrcDS->GetMetadataItem( "DTED_MaintenanceDescription" ) );
813 :
814 14 : if( poSrcDS->GetMetadataItem( "DTED_Producer" ) != NULL )
815 : DTEDSetMetadata( psDTED, DTEDMD_PRODUCER,
816 2 : poSrcDS->GetMetadataItem( "DTED_Producer" ) );
817 :
818 14 : if( poSrcDS->GetMetadataItem( "DTED_VerticalDatum" ) != NULL )
819 : DTEDSetMetadata( psDTED, DTEDMD_VERTDATUM,
820 2 : poSrcDS->GetMetadataItem( "DTED_VerticalDatum" ) );
821 :
822 14 : if( poSrcDS->GetMetadataItem( "DTED_HorizontalDatum" ) != NULL )
823 : DTEDSetMetadata( psDTED, DTEDMD_HORIZDATUM,
824 2 : poSrcDS->GetMetadataItem( "DTED_HorizontalDatum" ) );
825 :
826 14 : if( poSrcDS->GetMetadataItem( "DTED_DigitizingSystem" ) != NULL )
827 : DTEDSetMetadata( psDTED, DTEDMD_DIGITIZING_SYS,
828 2 : poSrcDS->GetMetadataItem( "DTED_DigitizingSystem" ) );
829 :
830 14 : if( poSrcDS->GetMetadataItem( "DTED_CompilationDate" ) != NULL )
831 : DTEDSetMetadata( psDTED, DTEDMD_COMPILATION_DATE,
832 2 : poSrcDS->GetMetadataItem( "DTED_CompilationDate" ) );
833 :
834 14 : if( poSrcDS->GetMetadataItem( "DTED_HorizontalAccuracy" ) != NULL )
835 : DTEDSetMetadata( psDTED, DTEDMD_HORIZACCURACY,
836 2 : poSrcDS->GetMetadataItem( "DTED_HorizontalAccuracy" ) );
837 :
838 14 : if( poSrcDS->GetMetadataItem( "DTED_RelHorizontalAccuracy" ) != NULL )
839 : DTEDSetMetadata( psDTED, DTEDMD_REL_HORIZACCURACY,
840 2 : poSrcDS->GetMetadataItem( "DTED_RelHorizontalAccuracy" ) );
841 :
842 14 : if( poSrcDS->GetMetadataItem( "DTED_RelVerticalAccuracy" ) != NULL )
843 : DTEDSetMetadata( psDTED, DTEDMD_REL_VERTACCURACY,
844 2 : poSrcDS->GetMetadataItem( "DTED_RelVerticalAccuracy" ) );
845 :
846 : /* -------------------------------------------------------------------- */
847 : /* Try to open the resulting DTED file. */
848 : /* -------------------------------------------------------------------- */
849 14 : DTEDClose( psDTED );
850 :
851 : /* -------------------------------------------------------------------- */
852 : /* Reopen and copy missing information into a PAM file. */
853 : /* -------------------------------------------------------------------- */
854 : GDALPamDataset *poDS = (GDALPamDataset *)
855 14 : GDALOpen( pszFilename, GA_ReadOnly );
856 :
857 14 : if( poDS )
858 14 : poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
859 :
860 14 : return poDS;
861 : }
862 :
863 : /************************************************************************/
864 : /* GDALRegister_DTED() */
865 : /************************************************************************/
866 :
867 409 : void GDALRegister_DTED()
868 :
869 : {
870 : GDALDriver *poDriver;
871 :
872 409 : if( GDALGetDriverByName( "DTED" ) == NULL )
873 : {
874 392 : poDriver = new GDALDriver();
875 :
876 392 : poDriver->SetDescription( "DTED" );
877 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
878 392 : "DTED Elevation Raster" );
879 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
880 392 : "frmt_various.html#DTED" );
881 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
882 392 : "Byte Int16 UInt16" );
883 :
884 392 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
885 :
886 392 : poDriver->pfnOpen = DTEDDataset::Open;
887 392 : poDriver->pfnCreateCopy = DTEDCreateCopy;
888 :
889 392 : GetGDALDriverManager()->RegisterDriver( poDriver );
890 : }
891 409 : }
892 :
|
