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 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 68 : {
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 68 : DTEDRasterBand::DTEDRasterBand( DTEDDataset *poDS, int nBand )
102 :
103 : {
104 68 : this->poDS = poDS;
105 68 : this->nBand = nBand;
106 :
107 68 : eDataType = GDT_Int16;
108 :
109 68 : bNoDataSet = TRUE;
110 68 : 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 68 : poDS->GetRasterXSize() : 1;
119 68 : nBlockYSize = poDS->GetRasterYSize();
120 68 : }
121 :
122 : /************************************************************************/
123 : /* IReadBlock() */
124 : /************************************************************************/
125 :
126 4631 : CPLErr DTEDRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
127 : void * pImage )
128 :
129 : {
130 4631 : DTEDDataset *poDTED_DS = (DTEDDataset *) poDS;
131 4631 : int nYSize = poDTED_DS->psDTED->nYSize;
132 : GInt16 *panData;
133 :
134 : (void) nBlockXOff;
135 4631 : CPLAssert( nBlockYOff == 0 );
136 :
137 4631 : 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 4631 : panData = (GInt16 *) pImage;
165 4631 : 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 611600 : for( int i = nYSize/2; i >= 0; i-- )
174 : {
175 : GInt16 nTemp;
176 :
177 606970 : nTemp = panData[i];
178 606970 : panData[i] = panData[nYSize - i - 1];
179 606970 : panData[nYSize - i - 1] = nTemp;
180 : }
181 :
182 4630 : 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 64 : double DTEDRasterBand::GetNoDataValue( int * pbSuccess )
237 :
238 : {
239 64 : if( pbSuccess )
240 55 : *pbSuccess = bNoDataSet;
241 :
242 64 : return dfNoDataValue;
243 : }
244 :
245 : /************************************************************************/
246 : /* ~DTEDDataset() */
247 : /************************************************************************/
248 :
249 68 : DTEDDataset::DTEDDataset()
250 : {
251 68 : pszFilename = CPLStrdup("unknown");
252 68 : pszProjection = CPLStrdup("");
253 68 : bVerifyChecksum = CSLTestBoolean(CPLGetConfigOption("DTED_VERIFY_CHECKSUM", "NO"));
254 68 : }
255 :
256 : /************************************************************************/
257 : /* ~DTEDDataset() */
258 : /************************************************************************/
259 :
260 68 : DTEDDataset::~DTEDDataset()
261 :
262 : {
263 68 : FlushCache();
264 68 : CPLFree(pszFilename);
265 68 : CPLFree( pszProjection );
266 68 : if( psDTED != NULL )
267 68 : DTEDClose( psDTED );
268 68 : }
269 :
270 : /************************************************************************/
271 : /* SetFileName() */
272 : /************************************************************************/
273 :
274 68 : void DTEDDataset::SetFileName(const char* pszFilename)
275 :
276 : {
277 68 : CPLFree(this->pszFilename);
278 68 : this->pszFilename = CPLStrdup(pszFilename);
279 68 : }
280 :
281 : /************************************************************************/
282 : /* Identify() */
283 : /************************************************************************/
284 :
285 12792 : 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 12792 : if( poOpenInfo->nHeaderBytes < 240 )
295 11094 : return FALSE;
296 :
297 1698 : 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 1630 : return FALSE;
302 : }
303 :
304 68 : int bFoundUHL = FALSE;
305 137 : for(i=0;i<poOpenInfo->nHeaderBytes-3 && !bFoundUHL ;i += DTED_UHL_SIZE)
306 : {
307 69 : if( EQUALN((const char *)poOpenInfo->pabyHeader + i,"UHL", 3) )
308 : {
309 68 : bFoundUHL = TRUE;
310 : }
311 : }
312 68 : if (!bFoundUHL)
313 0 : return FALSE;
314 :
315 68 : return TRUE;
316 : }
317 :
318 : /************************************************************************/
319 : /* Open() */
320 : /************************************************************************/
321 :
322 3513 : GDALDataset *DTEDDataset::Open( GDALOpenInfo * poOpenInfo )
323 :
324 : {
325 : int i;
326 : DTEDInfo *psDTED;
327 :
328 3513 : if (!Identify(poOpenInfo))
329 3445 : return NULL;
330 :
331 : /* -------------------------------------------------------------------- */
332 : /* Try opening the dataset. */
333 : /* -------------------------------------------------------------------- */
334 68 : psDTED = DTEDOpen( poOpenInfo->pszFilename, (poOpenInfo->eAccess == GA_Update) ? "rb+" : "rb", TRUE );
335 :
336 68 : if( psDTED == NULL )
337 0 : return( NULL );
338 :
339 : /* -------------------------------------------------------------------- */
340 : /* Create a corresponding GDALDataset. */
341 : /* -------------------------------------------------------------------- */
342 : DTEDDataset *poDS;
343 :
344 68 : poDS = new DTEDDataset();
345 68 : poDS->SetFileName(poOpenInfo->pszFilename);
346 :
347 68 : poDS->eAccess = poOpenInfo->eAccess;
348 68 : poDS->psDTED = psDTED;
349 :
350 : /* -------------------------------------------------------------------- */
351 : /* Capture some information from the file that is of interest. */
352 : /* -------------------------------------------------------------------- */
353 68 : poDS->nRasterXSize = psDTED->nXSize;
354 68 : poDS->nRasterYSize = psDTED->nYSize;
355 :
356 68 : 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 68 : poDS->nBands = 1;
366 272 : for( i = 0; i < poDS->nBands; i++ )
367 68 : poDS->SetBand( i+1, new DTEDRasterBand( poDS, i+1 ) );
368 :
369 : /* -------------------------------------------------------------------- */
370 : /* Collect any metadata available. */
371 : /* -------------------------------------------------------------------- */
372 : char *pszValue;
373 :
374 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_VERTACCURACY_UHL );
375 68 : poDS->SetMetadataItem( "DTED_VerticalAccuracy_UHL", pszValue );
376 68 : free( pszValue );
377 :
378 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_VERTACCURACY_ACC );
379 68 : poDS->SetMetadataItem( "DTED_VerticalAccuracy_ACC", pszValue );
380 68 : free( pszValue );
381 :
382 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_SECURITYCODE_UHL );
383 68 : poDS->SetMetadataItem( "DTED_SecurityCode_UHL", pszValue );
384 68 : free( pszValue );
385 :
386 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_SECURITYCODE_DSI );
387 68 : poDS->SetMetadataItem( "DTED_SecurityCode_DSI", pszValue );
388 68 : free( pszValue );
389 :
390 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_UNIQUEREF_UHL );
391 68 : poDS->SetMetadataItem( "DTED_UniqueRef_UHL", pszValue );
392 68 : free( pszValue );
393 :
394 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_UNIQUEREF_DSI );
395 68 : poDS->SetMetadataItem( "DTED_UniqueRef_DSI", pszValue );
396 68 : free( pszValue );
397 :
398 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_DATA_EDITION );
399 68 : poDS->SetMetadataItem( "DTED_DataEdition", pszValue );
400 68 : free( pszValue );
401 :
402 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_MATCHMERGE_VERSION );
403 68 : poDS->SetMetadataItem( "DTED_MatchMergeVersion", pszValue );
404 68 : free( pszValue );
405 :
406 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_MAINT_DATE );
407 68 : poDS->SetMetadataItem( "DTED_MaintenanceDate", pszValue );
408 68 : free( pszValue );
409 :
410 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_MATCHMERGE_DATE );
411 68 : poDS->SetMetadataItem( "DTED_MatchMergeDate", pszValue );
412 68 : free( pszValue );
413 :
414 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_MAINT_DESCRIPTION );
415 68 : poDS->SetMetadataItem( "DTED_MaintenanceDescription", pszValue );
416 68 : free( pszValue );
417 :
418 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_PRODUCER );
419 68 : poDS->SetMetadataItem( "DTED_Producer", pszValue );
420 68 : free( pszValue );
421 :
422 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_VERTDATUM );
423 68 : poDS->SetMetadataItem( "DTED_VerticalDatum", pszValue );
424 68 : free( pszValue );
425 :
426 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_HORIZDATUM );
427 68 : poDS->SetMetadataItem( "DTED_HorizontalDatum", pszValue );
428 68 : free( pszValue );
429 :
430 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_DIGITIZING_SYS );
431 68 : poDS->SetMetadataItem( "DTED_DigitizingSystem", pszValue );
432 68 : free( pszValue );
433 :
434 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_COMPILATION_DATE );
435 68 : poDS->SetMetadataItem( "DTED_CompilationDate", pszValue );
436 68 : free( pszValue );
437 :
438 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_HORIZACCURACY );
439 68 : poDS->SetMetadataItem( "DTED_HorizontalAccuracy", pszValue );
440 68 : free( pszValue );
441 :
442 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_REL_HORIZACCURACY );
443 68 : poDS->SetMetadataItem( "DTED_RelHorizontalAccuracy", pszValue );
444 68 : free( pszValue );
445 :
446 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_REL_VERTACCURACY );
447 68 : poDS->SetMetadataItem( "DTED_RelVerticalAccuracy", pszValue );
448 68 : free( pszValue );
449 :
450 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_ORIGINLAT );
451 68 : poDS->SetMetadataItem( "DTED_OriginLatitude", pszValue );
452 68 : free( pszValue );
453 :
454 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_ORIGINLONG );
455 68 : poDS->SetMetadataItem( "DTED_OriginLongitude", pszValue );
456 68 : free( pszValue );
457 :
458 68 : pszValue = DTEDGetMetadata( psDTED, DTEDMD_NIMA_DESIGNATOR );
459 68 : poDS->SetMetadataItem( "DTED_NimaDesignator", pszValue );
460 68 : free( pszValue );
461 :
462 68 : poDS->SetMetadataItem( GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT );
463 :
464 : /* -------------------------------------------------------------------- */
465 : /* Initialize any PAM information. */
466 : /* -------------------------------------------------------------------- */
467 68 : poDS->SetDescription( poOpenInfo->pszFilename );
468 68 : poDS->TryLoadXML();
469 :
470 : // if no SR in xml, try aux
471 68 : const char* pszPrj = poDS->GDALPamDataset::GetProjectionRef();
472 68 : if( !pszPrj || strlen(pszPrj) == 0 )
473 : {
474 68 : GDALDataset* poAuxDS = GDALFindAssociatedAuxFile( poOpenInfo->pszFilename, GA_ReadOnly, poDS );
475 68 : 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 68 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
492 68 : return( poDS );
493 : }
494 :
495 : /************************************************************************/
496 : /* GetGeoTransform() */
497 : /************************************************************************/
498 :
499 74 : CPLErr DTEDDataset::GetGeoTransform( double * padfTransform )
500 :
501 : {
502 74 : padfTransform[0] = psDTED->dfULCornerX;
503 74 : padfTransform[1] = psDTED->dfPixelSizeX;
504 74 : padfTransform[2] = 0.0;
505 74 : padfTransform[3] = psDTED->dfULCornerY;
506 74 : padfTransform[4] = 0.0;
507 74 : padfTransform[5] = psDTED->dfPixelSizeY * -1;
508 :
509 74 : return( CE_None );
510 : }
511 :
512 : /************************************************************************/
513 : /* GetProjectionRef() */
514 : /************************************************************************/
515 :
516 70 : const char *DTEDDataset::GetProjectionRef()
517 :
518 : {
519 : // get xml and aux SR first
520 70 : const char* pszPrj = GDALPamDataset::GetProjectionRef();
521 70 : if(pszPrj && strlen(pszPrj) > 0)
522 0 : return pszPrj;
523 :
524 70 : if (pszProjection && strlen(pszProjection) > 0)
525 0 : return pszProjection;
526 :
527 70 : pszPrj = GetMetadataItem( "DTED_HorizontalDatum");
528 70 : if (EQUAL(pszPrj, "WGS84"))
529 : {
530 69 : 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 1 : else if (EQUAL(pszPrj, "WGS72"))
533 : {
534 : static int bWarned = FALSE;
535 1 : if (!bWarned)
536 : {
537 1 : 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 1 : "No more warnings will be issued in this session about this operation.", GetFileName() );
545 : }
546 1 : 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 21 : 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 21 : int nBands = poSrcDS->GetRasterCount();
586 21 : if (nBands == 0)
587 : {
588 : CPLError( CE_Failure, CPLE_NotSupported,
589 1 : "DTED driver does not support source dataset with zero band.\n");
590 1 : return NULL;
591 : }
592 :
593 20 : if (nBands != 1)
594 : {
595 : CPLError( (bStrict) ? CE_Failure : CE_Warning, CPLE_NotSupported,
596 4 : "DTED driver only uses the first band of the dataset.\n");
597 4 : if (bStrict)
598 4 : return NULL;
599 : }
600 :
601 16 : if( pfnProgress && !pfnProgress( 0.0, NULL, pProgressData ) )
602 0 : return NULL;
603 :
604 : /* -------------------------------------------------------------------- */
605 : /* Work out the level. */
606 : /* -------------------------------------------------------------------- */
607 : int nLevel;
608 :
609 16 : if( poSrcDS->GetRasterYSize() == 121 )
610 2 : nLevel = 0;
611 14 : else if( poSrcDS->GetRasterYSize() == 1201 )
612 3 : nLevel = 1;
613 11 : else if( poSrcDS->GetRasterYSize() == 3601 )
614 0 : nLevel = 2;
615 : else
616 : {
617 : CPLError( CE_Warning, CPLE_AppDefined,
618 11 : "The source does not appear to be a properly formatted cell." );
619 11 : nLevel = 1;
620 : }
621 :
622 : /* -------------------------------------------------------------------- */
623 : /* Checks the input SRS */
624 : /* -------------------------------------------------------------------- */
625 16 : OGRSpatialReference ogrsr_input;
626 16 : OGRSpatialReference ogrsr_wgs84;
627 16 : char* c = (char*)poSrcDS->GetProjectionRef();
628 16 : ogrsr_input.importFromWkt(&c);
629 16 : ogrsr_wgs84.SetWellKnownGeogCS( "WGS84" );
630 16 : 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 16 : poSrcDS->GetGeoTransform( adfGeoTransform );
645 :
646 : nLLOriginLat = (int)
647 16 : floor(adfGeoTransform[3]
648 16 : + poSrcDS->GetRasterYSize() * adfGeoTransform[5] + 0.5);
649 :
650 16 : nLLOriginLong = (int) floor(adfGeoTransform[0] + 0.5);
651 :
652 35 : if (fabs(nLLOriginLat - (adfGeoTransform[3]
653 16 : + (poSrcDS->GetRasterYSize() - 0.5) * adfGeoTransform[5])) > 1e-10 ||
654 3 : 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 13 : "aligned on plain latitude/longitude boundaries.");
659 : }
660 :
661 : /* -------------------------------------------------------------------- */
662 : /* Check horizontal source size. */
663 : /* -------------------------------------------------------------------- */
664 : int expectedXSize;
665 16 : if( ABS(nLLOriginLat) >= 80 )
666 0 : expectedXSize = (poSrcDS->GetRasterYSize() - 1) / 6 + 1;
667 16 : else if( ABS(nLLOriginLat) >= 75 )
668 0 : expectedXSize = (poSrcDS->GetRasterYSize() - 1) / 4 + 1;
669 16 : else if( ABS(nLLOriginLat) >= 70 )
670 0 : expectedXSize = (poSrcDS->GetRasterYSize() - 1) / 3 + 1;
671 16 : else if( ABS(nLLOriginLat) >= 50 )
672 1 : expectedXSize = (poSrcDS->GetRasterYSize() - 1) / 2 + 1;
673 : else
674 15 : expectedXSize = poSrcDS->GetRasterYSize();
675 :
676 16 : 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 16 : pszError = DTEDCreate( pszFilename, nLevel, nLLOriginLat, nLLOriginLong );
690 :
691 16 : if( pszError != NULL )
692 : {
693 : CPLError( CE_Failure, CPLE_AppDefined,
694 2 : "%s", pszError );
695 2 : return NULL;
696 : }
697 :
698 : /* -------------------------------------------------------------------- */
699 : /* Open the DTED file so we can output the data to it. */
700 : /* -------------------------------------------------------------------- */
701 : DTEDInfo *psDTED;
702 :
703 14 : psDTED = DTEDOpen( pszFilename, "rb+", FALSE );
704 14 : if( psDTED == NULL )
705 0 : return NULL;
706 :
707 : /* -------------------------------------------------------------------- */
708 : /* Read all the data in a single buffer. */
709 : /* -------------------------------------------------------------------- */
710 14 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( 1 );
711 : GInt16 *panData;
712 :
713 : panData = (GInt16 *)
714 14 : VSIMalloc(sizeof(GInt16) * psDTED->nXSize * psDTED->nYSize);
715 14 : 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 14668 : 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 14654 : GDT_Int16, 0, 0 );
727 :
728 14654 : 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 14 : double srcBandNoData = poSrcBand->GetNoDataValue(&bSrcBandHasNoData);
740 :
741 : /* -------------------------------------------------------------------- */
742 : /* Write all the profiles. */
743 : /* -------------------------------------------------------------------- */
744 : GInt16 anProfData[3601];
745 14 : int dfNodataCount=0;
746 : GByte iPartialCell;
747 :
748 14068 : for( int iProfile = 0; iProfile < psDTED->nXSize; iProfile++ )
749 : {
750 16631548 : for( int iY = 0; iY < psDTED->nYSize; iY++ )
751 : {
752 16617494 : anProfData[iY] = panData[iProfile + iY * psDTED->nXSize];
753 16617494 : if ( bSrcBandHasNoData && anProfData[iY] == srcBandNoData)
754 : {
755 0 : anProfData[iY] = DTED_NODATA_VALUE;
756 0 : dfNodataCount++;
757 : }
758 16617494 : else if ( anProfData[iY] == DTED_NODATA_VALUE )
759 0 : dfNodataCount++;
760 : }
761 14054 : DTEDWriteProfile( psDTED, iProfile, anProfData );
762 :
763 14054 : 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 14 : CPLFree( panData );
775 :
776 : /* -------------------------------------------------------------------- */
777 : /* Partial cell indicator: 0 for complete coverage; 1-99 for incomplete */
778 : /* -------------------------------------------------------------------- */
779 : char szPartialCell[3];
780 :
781 14 : if ( dfNodataCount == 0 )
782 14 : 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 14 : sprintf(szPartialCell,"%02d",iPartialCell);
791 14 : strncpy((char *) (psDTED->pachDSIRecord+289), szPartialCell, 2 );
792 :
793 : /* -------------------------------------------------------------------- */
794 : /* Try to copy any matching available metadata. */
795 : /* -------------------------------------------------------------------- */
796 14 : if( poSrcDS->GetMetadataItem( "DTED_VerticalAccuracy_UHL" ) != NULL )
797 : DTEDSetMetadata( psDTED, DTEDMD_VERTACCURACY_UHL,
798 2 : poSrcDS->GetMetadataItem( "DTED_VerticalAccuracy_UHL" ) );
799 :
800 14 : if( poSrcDS->GetMetadataItem( "DTED_VerticalAccuracy_ACC" ) != NULL )
801 : DTEDSetMetadata( psDTED, DTEDMD_VERTACCURACY_ACC,
802 2 : poSrcDS->GetMetadataItem( "DTED_VerticalAccuracy_ACC" ) );
803 :
804 14 : if( poSrcDS->GetMetadataItem( "DTED_SecurityCode_UHL" ) != NULL )
805 : DTEDSetMetadata( psDTED, DTEDMD_SECURITYCODE_UHL,
806 2 : poSrcDS->GetMetadataItem( "DTED_SecurityCode_UHL" ) );
807 :
808 14 : if( poSrcDS->GetMetadataItem( "DTED_SecurityCode_DSI" ) != NULL )
809 : DTEDSetMetadata( psDTED, DTEDMD_SECURITYCODE_DSI,
810 2 : poSrcDS->GetMetadataItem( "DTED_SecurityCode_DSI" ) );
811 :
812 14 : if( poSrcDS->GetMetadataItem( "DTED_UniqueRef_UHL" ) != NULL )
813 : DTEDSetMetadata( psDTED, DTEDMD_UNIQUEREF_UHL,
814 2 : poSrcDS->GetMetadataItem( "DTED_UniqueRef_UHL" ) );
815 :
816 14 : if( poSrcDS->GetMetadataItem( "DTED_UniqueRef_DSI" ) != NULL )
817 : DTEDSetMetadata( psDTED, DTEDMD_UNIQUEREF_DSI,
818 2 : poSrcDS->GetMetadataItem( "DTED_UniqueRef_DSI" ) );
819 :
820 14 : if( poSrcDS->GetMetadataItem( "DTED_DataEdition" ) != NULL )
821 : DTEDSetMetadata( psDTED, DTEDMD_DATA_EDITION,
822 2 : poSrcDS->GetMetadataItem( "DTED_DataEdition" ) );
823 :
824 14 : if( poSrcDS->GetMetadataItem( "DTED_MatchMergeVersion" ) != NULL )
825 : DTEDSetMetadata( psDTED, DTEDMD_MATCHMERGE_VERSION,
826 2 : poSrcDS->GetMetadataItem( "DTED_MatchMergeVersion" ) );
827 :
828 14 : if( poSrcDS->GetMetadataItem( "DTED_MaintenanceDate" ) != NULL )
829 : DTEDSetMetadata( psDTED, DTEDMD_MAINT_DATE,
830 2 : poSrcDS->GetMetadataItem( "DTED_MaintenanceDate" ) );
831 :
832 14 : if( poSrcDS->GetMetadataItem( "DTED_MatchMergeDate" ) != NULL )
833 : DTEDSetMetadata( psDTED, DTEDMD_MATCHMERGE_DATE,
834 2 : poSrcDS->GetMetadataItem( "DTED_MatchMergeDate" ) );
835 :
836 14 : if( poSrcDS->GetMetadataItem( "DTED_MaintenanceDescription" ) != NULL )
837 : DTEDSetMetadata( psDTED, DTEDMD_MAINT_DESCRIPTION,
838 2 : poSrcDS->GetMetadataItem( "DTED_MaintenanceDescription" ) );
839 :
840 14 : if( poSrcDS->GetMetadataItem( "DTED_Producer" ) != NULL )
841 : DTEDSetMetadata( psDTED, DTEDMD_PRODUCER,
842 2 : poSrcDS->GetMetadataItem( "DTED_Producer" ) );
843 :
844 14 : if( poSrcDS->GetMetadataItem( "DTED_VerticalDatum" ) != NULL )
845 : DTEDSetMetadata( psDTED, DTEDMD_VERTDATUM,
846 2 : poSrcDS->GetMetadataItem( "DTED_VerticalDatum" ) );
847 :
848 14 : if( poSrcDS->GetMetadataItem( "DTED_HorizontalDatum" ) != NULL )
849 : DTEDSetMetadata( psDTED, DTEDMD_HORIZDATUM,
850 2 : poSrcDS->GetMetadataItem( "DTED_HorizontalDatum" ) );
851 :
852 14 : if( poSrcDS->GetMetadataItem( "DTED_DigitizingSystem" ) != NULL )
853 : DTEDSetMetadata( psDTED, DTEDMD_DIGITIZING_SYS,
854 2 : poSrcDS->GetMetadataItem( "DTED_DigitizingSystem" ) );
855 :
856 14 : if( poSrcDS->GetMetadataItem( "DTED_CompilationDate" ) != NULL )
857 : DTEDSetMetadata( psDTED, DTEDMD_COMPILATION_DATE,
858 2 : poSrcDS->GetMetadataItem( "DTED_CompilationDate" ) );
859 :
860 14 : if( poSrcDS->GetMetadataItem( "DTED_HorizontalAccuracy" ) != NULL )
861 : DTEDSetMetadata( psDTED, DTEDMD_HORIZACCURACY,
862 2 : poSrcDS->GetMetadataItem( "DTED_HorizontalAccuracy" ) );
863 :
864 14 : if( poSrcDS->GetMetadataItem( "DTED_RelHorizontalAccuracy" ) != NULL )
865 : DTEDSetMetadata( psDTED, DTEDMD_REL_HORIZACCURACY,
866 2 : poSrcDS->GetMetadataItem( "DTED_RelHorizontalAccuracy" ) );
867 :
868 14 : if( poSrcDS->GetMetadataItem( "DTED_RelVerticalAccuracy" ) != NULL )
869 : DTEDSetMetadata( psDTED, DTEDMD_REL_VERTACCURACY,
870 2 : poSrcDS->GetMetadataItem( "DTED_RelVerticalAccuracy" ) );
871 :
872 : /* -------------------------------------------------------------------- */
873 : /* Try to open the resulting DTED file. */
874 : /* -------------------------------------------------------------------- */
875 14 : DTEDClose( psDTED );
876 :
877 : /* -------------------------------------------------------------------- */
878 : /* Reopen and copy missing information into a PAM file. */
879 : /* -------------------------------------------------------------------- */
880 : GDALPamDataset *poDS = (GDALPamDataset *)
881 14 : GDALOpen( pszFilename, GA_ReadOnly );
882 :
883 14 : if( poDS )
884 14 : poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
885 :
886 14 : return poDS;
887 : }
888 :
889 : /************************************************************************/
890 : /* GDALRegister_DTED() */
891 : /************************************************************************/
892 :
893 558 : void GDALRegister_DTED()
894 :
895 : {
896 : GDALDriver *poDriver;
897 :
898 558 : if( GDALGetDriverByName( "DTED" ) == NULL )
899 : {
900 537 : poDriver = new GDALDriver();
901 :
902 537 : poDriver->SetDescription( "DTED" );
903 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
904 537 : "DTED Elevation Raster" );
905 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
906 537 : "frmt_various.html#DTED" );
907 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
908 537 : "Byte Int16 UInt16" );
909 :
910 537 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
911 :
912 537 : poDriver->pfnOpen = DTEDDataset::Open;
913 537 : poDriver->pfnIdentify = DTEDDataset::Identify;
914 537 : poDriver->pfnCreateCopy = DTEDCreateCopy;
915 :
916 537 : GetGDALDriverManager()->RegisterDriver( poDriver );
917 : }
918 558 : }
919 :
|
