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