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