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