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