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