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