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 "ogr_spatialref.h"
18 :
19 : #include <cstdlib>
20 : #include <algorithm>
21 :
22 : /************************************************************************/
23 : /* ==================================================================== */
24 : /* DTEDDataset */
25 : /* ==================================================================== */
26 : /************************************************************************/
27 :
28 : class DTEDRasterBand;
29 :
30 : class DTEDDataset final : public GDALPamDataset
31 : {
32 : friend class DTEDRasterBand;
33 :
34 : char *pszFilename;
35 : DTEDInfo *psDTED;
36 : int bVerifyChecksum;
37 : mutable OGRSpatialReference m_oSRS{};
38 :
39 : public:
40 : DTEDDataset();
41 : ~DTEDDataset() override;
42 :
43 : const OGRSpatialReference *GetSpatialRef() const override;
44 : CPLErr GetGeoTransform(double *) override;
45 :
46 1 : const char *GetFileName() const
47 : {
48 1 : return pszFilename;
49 : }
50 :
51 : void SetFileName(const char *pszFilename);
52 :
53 : static GDALDataset *Open(GDALOpenInfo *);
54 : static int Identify(GDALOpenInfo *);
55 : };
56 :
57 : /************************************************************************/
58 : /* ==================================================================== */
59 : /* DTEDRasterBand */
60 : /* ==================================================================== */
61 : /************************************************************************/
62 :
63 : class DTEDRasterBand final : public GDALPamRasterBand
64 : {
65 : friend class DTEDDataset;
66 :
67 : int bNoDataSet;
68 : double dfNoDataValue;
69 :
70 : public:
71 : DTEDRasterBand(DTEDDataset *, int);
72 :
73 : CPLErr IReadBlock(int, int, void *) override;
74 : CPLErr IWriteBlock(int, int, void *) override;
75 :
76 : double GetNoDataValue(int *pbSuccess = nullptr) override;
77 :
78 26 : const char *GetUnitType() override
79 : {
80 26 : return "m";
81 : }
82 : };
83 :
84 : /************************************************************************/
85 : /* DTEDRasterBand() */
86 : /************************************************************************/
87 :
88 51 : DTEDRasterBand::DTEDRasterBand(DTEDDataset *poDSIn, int nBandIn)
89 51 : : bNoDataSet(TRUE), dfNoDataValue(static_cast<double>(DTED_NODATA_VALUE))
90 : {
91 51 : poDS = poDSIn;
92 51 : nBand = nBandIn;
93 :
94 51 : eDataType = GDT_Int16;
95 :
96 : /* For some applications, it may be valuable to consider the whole DTED */
97 : /* file as single block, as the column orientation doesn't fit very well */
98 : /* with some scanline oriented algorithms */
99 : /* Of course you need to have a big enough case size, particularly for DTED
100 : * 2 */
101 : /* datasets */
102 51 : nBlockXSize =
103 51 : CPLTestBool(CPLGetConfigOption("GDAL_DTED_SINGLE_BLOCK", "NO"))
104 51 : ? poDS->GetRasterXSize()
105 : : 1;
106 51 : nBlockYSize = poDS->GetRasterYSize();
107 51 : }
108 :
109 : /************************************************************************/
110 : /* IReadBlock() */
111 : /************************************************************************/
112 :
113 3911 : CPLErr DTEDRasterBand::IReadBlock(int nBlockXOff, CPL_UNUSED int nBlockYOff,
114 : void *pImage)
115 : {
116 3911 : DTEDDataset *poDTED_DS = (DTEDDataset *)poDS;
117 3911 : int nYSize = poDTED_DS->psDTED->nYSize;
118 : GInt16 *panData;
119 :
120 : (void)nBlockXOff;
121 3911 : CPLAssert(nBlockYOff == 0);
122 :
123 3911 : if (nBlockXSize != 1)
124 : {
125 1 : const int cbs = 32; // optimize for 64 byte cache line size
126 1 : const int bsy = (nBlockYSize + cbs - 1) / cbs * cbs;
127 1 : panData = (GInt16 *)pImage;
128 1 : GInt16 *panBuffer = (GInt16 *)CPLMalloc(sizeof(GInt16) * cbs * bsy);
129 5 : for (int i = 0; i < nBlockXSize; i += cbs)
130 : {
131 4 : int n = std::min(cbs, nBlockXSize - i);
132 125 : for (int j = 0; j < n; ++j)
133 : {
134 121 : if (!DTEDReadProfileEx(poDTED_DS->psDTED, i + j,
135 121 : panBuffer + j * bsy,
136 : poDTED_DS->bVerifyChecksum))
137 : {
138 0 : CPLFree(panBuffer);
139 0 : return CE_Failure;
140 : }
141 : }
142 488 : for (int y = 0; y < nBlockYSize; ++y)
143 : {
144 484 : GInt16 *dst = panData + i + (nYSize - y - 1) * nBlockXSize;
145 484 : GInt16 *src = panBuffer + y;
146 15125 : for (int j = 0; j < n; ++j)
147 : {
148 14641 : dst[j] = src[j * bsy];
149 : }
150 : }
151 : }
152 :
153 1 : CPLFree(panBuffer);
154 1 : return CE_None;
155 : }
156 :
157 : /* -------------------------------------------------------------------- */
158 : /* Read the data. */
159 : /* -------------------------------------------------------------------- */
160 3910 : panData = (GInt16 *)pImage;
161 3910 : if (!DTEDReadProfileEx(poDTED_DS->psDTED, nBlockXOff, panData,
162 : poDTED_DS->bVerifyChecksum))
163 1 : return CE_Failure;
164 :
165 : /* -------------------------------------------------------------------- */
166 : /* Flip line to orient it top to bottom instead of bottom to */
167 : /* top. */
168 : /* -------------------------------------------------------------------- */
169 566898 : for (int i = nYSize / 2; i >= 0; i--)
170 : {
171 562989 : std::swap(panData[i], panData[nYSize - i - 1]);
172 : }
173 :
174 3909 : return CE_None;
175 : }
176 :
177 : /************************************************************************/
178 : /* IReadBlock() */
179 : /************************************************************************/
180 :
181 0 : CPLErr DTEDRasterBand::IWriteBlock(int nBlockXOff, CPL_UNUSED int nBlockYOff,
182 : void *pImage)
183 : {
184 0 : DTEDDataset *poDTED_DS = (DTEDDataset *)poDS;
185 : GInt16 *panData;
186 :
187 : (void)nBlockXOff;
188 0 : CPLAssert(nBlockYOff == 0);
189 :
190 0 : if (poDTED_DS->eAccess != GA_Update)
191 0 : return CE_Failure;
192 :
193 0 : if (nBlockXSize != 1)
194 : {
195 0 : panData = (GInt16 *)pImage;
196 0 : GInt16 *panBuffer = (GInt16 *)CPLMalloc(sizeof(GInt16) * nBlockYSize);
197 0 : for (int i = 0; i < nBlockXSize; i++)
198 : {
199 0 : for (int j = 0; j < nBlockYSize; j++)
200 : {
201 0 : panBuffer[j] = panData[j * nBlockXSize + i];
202 : }
203 0 : if (!DTEDWriteProfile(poDTED_DS->psDTED, i, panBuffer))
204 : {
205 0 : CPLFree(panBuffer);
206 0 : return CE_Failure;
207 : }
208 : }
209 :
210 0 : CPLFree(panBuffer);
211 0 : return CE_None;
212 : }
213 :
214 0 : panData = (GInt16 *)pImage;
215 0 : if (!DTEDWriteProfile(poDTED_DS->psDTED, nBlockXOff, panData))
216 0 : return CE_Failure;
217 :
218 0 : return CE_None;
219 : }
220 :
221 : /************************************************************************/
222 : /* GetNoDataValue() */
223 : /************************************************************************/
224 :
225 64 : double DTEDRasterBand::GetNoDataValue(int *pbSuccess)
226 :
227 : {
228 64 : if (pbSuccess)
229 53 : *pbSuccess = bNoDataSet;
230 :
231 64 : return dfNoDataValue;
232 : }
233 :
234 : /************************************************************************/
235 : /* ~DTEDDataset() */
236 : /************************************************************************/
237 :
238 51 : DTEDDataset::DTEDDataset()
239 102 : : pszFilename(CPLStrdup("unknown")), psDTED(nullptr),
240 : bVerifyChecksum(
241 51 : CPLTestBool(CPLGetConfigOption("DTED_VERIFY_CHECKSUM", "NO")))
242 : {
243 51 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
244 51 : }
245 :
246 : /************************************************************************/
247 : /* ~DTEDDataset() */
248 : /************************************************************************/
249 :
250 102 : DTEDDataset::~DTEDDataset()
251 :
252 : {
253 51 : FlushCache(true);
254 51 : CPLFree(pszFilename);
255 51 : if (psDTED != nullptr)
256 51 : DTEDClose(psDTED);
257 102 : }
258 :
259 : /************************************************************************/
260 : /* SetFileName() */
261 : /************************************************************************/
262 :
263 51 : void DTEDDataset::SetFileName(const char *pszFilenameIn)
264 :
265 : {
266 51 : CPLFree(this->pszFilename);
267 51 : this->pszFilename = CPLStrdup(pszFilenameIn);
268 51 : }
269 :
270 : /************************************************************************/
271 : /* Identify() */
272 : /************************************************************************/
273 :
274 56881 : int DTEDDataset::Identify(GDALOpenInfo *poOpenInfo)
275 :
276 : {
277 : /* -------------------------------------------------------------------- */
278 : /* Does the file start with one of the possible DTED header */
279 : /* record types, and do we have a UHL marker? */
280 : /* -------------------------------------------------------------------- */
281 56881 : if (poOpenInfo->nHeaderBytes < 240)
282 50496 : return FALSE;
283 :
284 6385 : if (!STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "VOL") &&
285 6385 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "HDR") &&
286 6383 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "UHL"))
287 : {
288 6283 : return FALSE;
289 : }
290 :
291 102 : bool bFoundUHL = false;
292 206 : for (int i = 0; i < poOpenInfo->nHeaderBytes - 3 && !bFoundUHL;
293 104 : i += DTED_UHL_SIZE)
294 : {
295 104 : if (STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + i, "UHL"))
296 : {
297 102 : bFoundUHL = true;
298 : }
299 : }
300 102 : if (!bFoundUHL)
301 0 : return FALSE;
302 :
303 102 : return TRUE;
304 : }
305 :
306 : /************************************************************************/
307 : /* Open() */
308 : /************************************************************************/
309 :
310 51 : GDALDataset *DTEDDataset::Open(GDALOpenInfo *poOpenInfo)
311 :
312 : {
313 51 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
314 0 : return nullptr;
315 :
316 : /* -------------------------------------------------------------------- */
317 : /* Try opening the dataset. */
318 : /* -------------------------------------------------------------------- */
319 51 : VSILFILE *fp = poOpenInfo->fpL;
320 51 : poOpenInfo->fpL = nullptr;
321 : DTEDInfo *psDTED =
322 51 : DTEDOpenEx(fp, poOpenInfo->pszFilename,
323 51 : (poOpenInfo->eAccess == GA_Update) ? "rb+" : "rb", TRUE);
324 :
325 51 : if (psDTED == nullptr)
326 0 : return nullptr;
327 :
328 : /* -------------------------------------------------------------------- */
329 : /* Create a corresponding GDALDataset. */
330 : /* -------------------------------------------------------------------- */
331 51 : DTEDDataset *poDS = new DTEDDataset();
332 51 : poDS->SetFileName(poOpenInfo->pszFilename);
333 :
334 51 : poDS->eAccess = poOpenInfo->eAccess;
335 51 : poDS->psDTED = psDTED;
336 :
337 : /* -------------------------------------------------------------------- */
338 : /* Capture some information from the file that is of interest. */
339 : /* -------------------------------------------------------------------- */
340 51 : poDS->nRasterXSize = psDTED->nXSize;
341 51 : poDS->nRasterYSize = psDTED->nYSize;
342 :
343 51 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
344 : {
345 0 : delete poDS;
346 0 : return nullptr;
347 : }
348 :
349 : /* -------------------------------------------------------------------- */
350 : /* Create band information objects. */
351 : /* -------------------------------------------------------------------- */
352 51 : poDS->nBands = 1;
353 102 : for (int i = 0; i < poDS->nBands; i++)
354 51 : poDS->SetBand(i + 1, new DTEDRasterBand(poDS, i + 1));
355 :
356 : /* -------------------------------------------------------------------- */
357 : /* Collect any metadata available. */
358 : /* -------------------------------------------------------------------- */
359 51 : char *pszValue = DTEDGetMetadata(psDTED, DTEDMD_VERTACCURACY_UHL);
360 51 : poDS->SetMetadataItem("DTED_VerticalAccuracy_UHL", pszValue);
361 51 : CPLFree(pszValue);
362 :
363 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_VERTACCURACY_ACC);
364 51 : poDS->SetMetadataItem("DTED_VerticalAccuracy_ACC", pszValue);
365 51 : CPLFree(pszValue);
366 :
367 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_SECURITYCODE_UHL);
368 51 : poDS->SetMetadataItem("DTED_SecurityCode_UHL", pszValue);
369 51 : CPLFree(pszValue);
370 :
371 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_SECURITYCODE_DSI);
372 51 : poDS->SetMetadataItem("DTED_SecurityCode_DSI", pszValue);
373 51 : CPLFree(pszValue);
374 :
375 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_UNIQUEREF_UHL);
376 51 : poDS->SetMetadataItem("DTED_UniqueRef_UHL", pszValue);
377 51 : CPLFree(pszValue);
378 :
379 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_UNIQUEREF_DSI);
380 51 : poDS->SetMetadataItem("DTED_UniqueRef_DSI", pszValue);
381 51 : CPLFree(pszValue);
382 :
383 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_DATA_EDITION);
384 51 : poDS->SetMetadataItem("DTED_DataEdition", pszValue);
385 51 : CPLFree(pszValue);
386 :
387 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_MATCHMERGE_VERSION);
388 51 : poDS->SetMetadataItem("DTED_MatchMergeVersion", pszValue);
389 51 : CPLFree(pszValue);
390 :
391 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_MAINT_DATE);
392 51 : poDS->SetMetadataItem("DTED_MaintenanceDate", pszValue);
393 51 : CPLFree(pszValue);
394 :
395 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_MATCHMERGE_DATE);
396 51 : poDS->SetMetadataItem("DTED_MatchMergeDate", pszValue);
397 51 : CPLFree(pszValue);
398 :
399 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_MAINT_DESCRIPTION);
400 51 : poDS->SetMetadataItem("DTED_MaintenanceDescription", pszValue);
401 51 : CPLFree(pszValue);
402 :
403 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_PRODUCER);
404 51 : poDS->SetMetadataItem("DTED_Producer", pszValue);
405 51 : CPLFree(pszValue);
406 :
407 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_VERTDATUM);
408 51 : poDS->SetMetadataItem("DTED_VerticalDatum", pszValue);
409 51 : CPLFree(pszValue);
410 :
411 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_HORIZDATUM);
412 51 : poDS->SetMetadataItem("DTED_HorizontalDatum", pszValue);
413 51 : CPLFree(pszValue);
414 :
415 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_DIGITIZING_SYS);
416 51 : poDS->SetMetadataItem("DTED_DigitizingSystem", pszValue);
417 51 : CPLFree(pszValue);
418 :
419 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_COMPILATION_DATE);
420 51 : poDS->SetMetadataItem("DTED_CompilationDate", pszValue);
421 51 : CPLFree(pszValue);
422 :
423 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_HORIZACCURACY);
424 51 : poDS->SetMetadataItem("DTED_HorizontalAccuracy", pszValue);
425 51 : CPLFree(pszValue);
426 :
427 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_REL_HORIZACCURACY);
428 51 : poDS->SetMetadataItem("DTED_RelHorizontalAccuracy", pszValue);
429 51 : CPLFree(pszValue);
430 :
431 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_REL_VERTACCURACY);
432 51 : poDS->SetMetadataItem("DTED_RelVerticalAccuracy", pszValue);
433 51 : CPLFree(pszValue);
434 :
435 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_ORIGINLAT);
436 51 : poDS->SetMetadataItem("DTED_OriginLatitude", pszValue);
437 51 : CPLFree(pszValue);
438 :
439 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_ORIGINLONG);
440 51 : poDS->SetMetadataItem("DTED_OriginLongitude", pszValue);
441 51 : CPLFree(pszValue);
442 :
443 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_NIMA_DESIGNATOR);
444 51 : poDS->SetMetadataItem("DTED_NimaDesignator", pszValue);
445 51 : CPLFree(pszValue);
446 :
447 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_PARTIALCELL_DSI);
448 51 : poDS->SetMetadataItem("DTED_PartialCellIndicator", pszValue);
449 51 : CPLFree(pszValue);
450 :
451 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_SECURITYCONTROL);
452 51 : poDS->SetMetadataItem("DTED_SecurityControl", pszValue);
453 51 : CPLFree(pszValue);
454 :
455 51 : pszValue = DTEDGetMetadata(psDTED, DTEDMD_SECURITYHANDLING);
456 51 : poDS->SetMetadataItem("DTED_SecurityHandling", pszValue);
457 51 : CPLFree(pszValue);
458 :
459 51 : poDS->SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
460 :
461 : /* -------------------------------------------------------------------- */
462 : /* Initialize any PAM information. */
463 : /* -------------------------------------------------------------------- */
464 51 : poDS->SetDescription(poOpenInfo->pszFilename);
465 51 : poDS->TryLoadXML(poOpenInfo->GetSiblingFiles());
466 :
467 : // if no SR in xml, try aux
468 51 : if (poDS->GDALPamDataset::GetSpatialRef() == nullptr)
469 : {
470 51 : int bTryAux = TRUE;
471 51 : if (poOpenInfo->GetSiblingFiles() != nullptr &&
472 51 : CSLFindString(poOpenInfo->GetSiblingFiles(),
473 : CPLResetExtension(
474 51 : CPLGetFilename(poOpenInfo->pszFilename), "aux")) <
475 102 : 0 &&
476 102 : CSLFindString(
477 51 : poOpenInfo->GetSiblingFiles(),
478 51 : CPLSPrintf("%s.aux", CPLGetFilename(poOpenInfo->pszFilename))) <
479 : 0)
480 51 : bTryAux = FALSE;
481 51 : if (bTryAux)
482 : {
483 0 : GDALDataset *poAuxDS = GDALFindAssociatedAuxFile(
484 0 : poOpenInfo->pszFilename, GA_ReadOnly, poDS);
485 0 : if (poAuxDS)
486 : {
487 0 : const auto poSRS = poAuxDS->GetSpatialRef();
488 0 : if (poSRS)
489 : {
490 0 : poDS->m_oSRS = *poSRS;
491 : }
492 :
493 0 : GDALClose(poAuxDS);
494 : }
495 : }
496 : }
497 :
498 : /* -------------------------------------------------------------------- */
499 : /* Support overviews. */
500 : /* -------------------------------------------------------------------- */
501 102 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename,
502 51 : poOpenInfo->GetSiblingFiles());
503 51 : return poDS;
504 : }
505 :
506 : /************************************************************************/
507 : /* GetGeoTransform() */
508 : /************************************************************************/
509 :
510 43 : CPLErr DTEDDataset::GetGeoTransform(double *padfTransform)
511 :
512 : {
513 :
514 : bool bApplyPixelIsPoint =
515 43 : CPLTestBool(CPLGetConfigOption("DTED_APPLY_PIXEL_IS_POINT", "FALSE"));
516 43 : if (!bApplyPixelIsPoint)
517 : {
518 42 : padfTransform[0] = psDTED->dfULCornerX;
519 42 : padfTransform[1] = psDTED->dfPixelSizeX;
520 42 : padfTransform[2] = 0.0;
521 42 : padfTransform[3] = psDTED->dfULCornerY;
522 42 : padfTransform[4] = 0.0;
523 42 : padfTransform[5] = psDTED->dfPixelSizeY * -1;
524 :
525 42 : return CE_None;
526 : }
527 : else
528 : {
529 1 : padfTransform[0] = psDTED->dfULCornerX + (0.5 * psDTED->dfPixelSizeX);
530 1 : padfTransform[1] = psDTED->dfPixelSizeX;
531 1 : padfTransform[2] = 0.0;
532 1 : padfTransform[3] = psDTED->dfULCornerY - (0.5 * psDTED->dfPixelSizeY);
533 1 : padfTransform[4] = 0.0;
534 1 : padfTransform[5] = psDTED->dfPixelSizeY * -1;
535 :
536 1 : return CE_None;
537 : }
538 : }
539 :
540 : /************************************************************************/
541 : /* GetSpatialRef() */
542 : /************************************************************************/
543 :
544 39 : const OGRSpatialReference *DTEDDataset::GetSpatialRef() const
545 :
546 : {
547 39 : if (!m_oSRS.IsEmpty())
548 20 : return &m_oSRS;
549 :
550 : // get xml and aux SR first
551 19 : const auto poSRS = GDALPamDataset::GetSpatialRef();
552 19 : if (poSRS)
553 : {
554 0 : m_oSRS = *poSRS;
555 0 : return &m_oSRS;
556 : }
557 :
558 : const char *pszVertDatum;
559 38 : const char *pszPrj = const_cast<DTEDDataset *>(this)->GetMetadataItem(
560 19 : "DTED_HorizontalDatum");
561 19 : if (EQUAL(pszPrj, "WGS84"))
562 : {
563 :
564 36 : pszVertDatum = const_cast<DTEDDataset *>(this)->GetMetadataItem(
565 18 : "DTED_VerticalDatum");
566 36 : if ((EQUAL(pszVertDatum, "MSL") || EQUAL(pszVertDatum, "E96")) &&
567 18 : CPLTestBool(CPLGetConfigOption("REPORT_COMPD_CS", "NO")))
568 : {
569 0 : m_oSRS.importFromWkt(
570 : "COMPD_CS[\"WGS 84 + EGM96 geoid height\", GEOGCS[\"WGS 84\", "
571 : "DATUM[\"WGS_1984\", SPHEROID[\"WGS "
572 : "84\",6378137,298.257223563, AUTHORITY[\"EPSG\",\"7030\"]], "
573 : "AUTHORITY[\"EPSG\",\"6326\"]], PRIMEM[\"Greenwich\",0, "
574 : "AUTHORITY[\"EPSG\",\"8901\"]], "
575 : "UNIT[\"degree\",0.0174532925199433, "
576 : "AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS["
577 : "\"Longitude\",EAST], AUTHORITY[\"EPSG\",\"4326\"]], "
578 : "VERT_CS[\"EGM96 geoid height\", VERT_DATUM[\"EGM96 "
579 : "geoid\",2005, AUTHORITY[\"EPSG\",\"5171\"]], "
580 : "UNIT[\"metre\",1, AUTHORITY[\"EPSG\",\"9001\"]], "
581 : "AXIS[\"Up\",UP], AUTHORITY[\"EPSG\",\"5773\"]]]");
582 : }
583 : // Support DTED with EGM08 vertical datum reference
584 18 : else if ((EQUAL(pszVertDatum, "E08")) &&
585 0 : CPLTestBool(CPLGetConfigOption("REPORT_COMPD_CS", "NO")))
586 : {
587 0 : m_oSRS.importFromWkt(
588 : "COMPD_CS[\"WGS 84 + EGM2008 height\",GEOGCS[\"WGS "
589 : "84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS "
590 : "84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],"
591 : "AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,"
592 : "AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0."
593 : "0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY["
594 : "\"EPSG\",\"4326\"]],VERT_CS[\"EGM2008 "
595 : "height\",VERT_DATUM[\"EGM2008 "
596 : "geoid\",2005,AUTHORITY[\"EPSG\",\"1027\"]],UNIT[\"metre\",1,"
597 : "AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Gravity-related "
598 : "height\",UP],AUTHORITY[\"EPSG\",\"3855\"]]]");
599 : }
600 : else
601 : {
602 18 : m_oSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
603 : }
604 : }
605 1 : else if (EQUAL(pszPrj, "WGS72"))
606 : {
607 : static bool bWarned = false;
608 1 : if (!bWarned)
609 : {
610 1 : bWarned = true;
611 1 : CPLError(CE_Warning, CPLE_AppDefined,
612 : "The DTED file %s indicates WGS72 as horizontal datum. \n"
613 : "As this is outdated nowadays, you should contact your "
614 : "data producer to get data georeferenced in WGS84.\n"
615 : "In some cases, WGS72 is a wrong indication and the "
616 : "georeferencing is really WGS84. In that case\n"
617 : "you might consider doing 'gdal_translate -of DTED -mo "
618 : "\"DTED_HorizontalDatum=WGS84\" src.dtX dst.dtX' to\n"
619 : "fix the DTED file.\n"
620 : "No more warnings will be issued in this session about "
621 : "this operation.",
622 : GetFileName());
623 : }
624 1 : m_oSRS.importFromWkt(
625 : "GEOGCS[\"WGS 72\",DATUM[\"WGS_1972\",SPHEROID[\"WGS "
626 : "72\",6378135,298.26]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0."
627 : "0174532925199433],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\","
628 : "EAST],AUTHORITY[\"EPSG\",\"4322\"]]");
629 : }
630 : else
631 : {
632 : static bool bWarned = false;
633 0 : if (!bWarned)
634 : {
635 0 : bWarned = true;
636 0 : CPLError(CE_Warning, CPLE_AppDefined,
637 : "The DTED file %s indicates %s as horizontal datum, which "
638 : "is not recognized by the DTED driver. \n"
639 : "The DTED driver is going to consider it as WGS84.\n"
640 : "No more warnings will be issued in this session about "
641 : "this operation.",
642 : GetFileName(), pszPrj);
643 : }
644 0 : m_oSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
645 : }
646 19 : return &m_oSRS;
647 : }
648 :
649 : /************************************************************************/
650 : /* DTEDCreateCopy() */
651 : /* */
652 : /* For now we will assume the input is exactly one proper */
653 : /* cell. */
654 : /************************************************************************/
655 :
656 22 : static GDALDataset *DTEDCreateCopy(const char *pszFilename,
657 : GDALDataset *poSrcDS, int bStrict,
658 : char ** /* papszOptions */,
659 : GDALProgressFunc pfnProgress,
660 : void *pProgressData)
661 :
662 : {
663 : /* -------------------------------------------------------------------- */
664 : /* Some some rudimentary checks */
665 : /* -------------------------------------------------------------------- */
666 22 : const int nBands = poSrcDS->GetRasterCount();
667 22 : if (nBands == 0)
668 : {
669 1 : CPLError(
670 : CE_Failure, CPLE_NotSupported,
671 : "DTED driver does not support source dataset with zero band.\n");
672 1 : return nullptr;
673 : }
674 :
675 21 : if (nBands != 1)
676 : {
677 4 : CPLError((bStrict) ? CE_Failure : CE_Warning, CPLE_NotSupported,
678 : "DTED driver only uses the first band of the dataset.\n");
679 4 : if (bStrict)
680 4 : return nullptr;
681 : }
682 :
683 17 : if (pfnProgress && !pfnProgress(0.0, nullptr, pProgressData))
684 0 : return nullptr;
685 :
686 : /* -------------------------------------------------------------------- */
687 : /* Work out the level. */
688 : /* -------------------------------------------------------------------- */
689 : int nLevel;
690 :
691 17 : if (poSrcDS->GetRasterYSize() == 121)
692 3 : nLevel = 0;
693 14 : else if (poSrcDS->GetRasterYSize() == 1201)
694 3 : nLevel = 1;
695 11 : else if (poSrcDS->GetRasterYSize() == 3601)
696 0 : nLevel = 2;
697 : else
698 : {
699 11 : CPLError(CE_Warning, CPLE_AppDefined,
700 : "The source does not appear to be a properly formatted cell.");
701 11 : nLevel = 1;
702 : }
703 :
704 : /* -------------------------------------------------------------------- */
705 : /* Checks the input SRS */
706 : /* -------------------------------------------------------------------- */
707 34 : OGRSpatialReference ogrsr_input;
708 17 : ogrsr_input.importFromWkt(poSrcDS->GetProjectionRef());
709 34 : OGRSpatialReference ogrsr_wgs84;
710 17 : ogrsr_wgs84.SetWellKnownGeogCS("WGS84");
711 17 : if (ogrsr_input.IsSameGeogCS(&ogrsr_wgs84) == FALSE)
712 : {
713 0 : CPLError(CE_Warning, CPLE_AppDefined,
714 : "The source projection coordinate system is %s. Only WGS 84 "
715 : "is supported.\n"
716 : "The DTED driver will generate a file as if the source was "
717 : "WGS 84 projection coordinate system.",
718 : poSrcDS->GetProjectionRef());
719 : }
720 :
721 : /* -------------------------------------------------------------------- */
722 : /* Work out the LL origin. */
723 : /* -------------------------------------------------------------------- */
724 : double adfGeoTransform[6];
725 :
726 17 : poSrcDS->GetGeoTransform(adfGeoTransform);
727 :
728 : int nLLOriginLat =
729 34 : (int)floor(adfGeoTransform[3] +
730 17 : poSrcDS->GetRasterYSize() * adfGeoTransform[5] + 0.5);
731 :
732 17 : int nLLOriginLong = (int)floor(adfGeoTransform[0] + 0.5);
733 :
734 51 : if (fabs(nLLOriginLat -
735 17 : (adfGeoTransform[3] + (poSrcDS->GetRasterYSize() - 0.5) *
736 21 : adfGeoTransform[5])) > 1e-10 ||
737 4 : fabs(nLLOriginLong - (adfGeoTransform[0] + 0.5 * adfGeoTransform[1])) >
738 : 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 1595 : void GDALRegister_DTED()
982 :
983 : {
984 1595 : if (GDALGetDriverByName("DTED") != nullptr)
985 302 : return;
986 :
987 1293 : GDALDriver *poDriver = new GDALDriver();
988 :
989 1293 : poDriver->SetDescription("DTED");
990 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
991 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "DTED Elevation Raster");
992 1293 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSIONS, "dt0 dt1 dt2");
993 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/dted.html");
994 1293 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Byte Int16 UInt16");
995 :
996 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
997 :
998 1293 : poDriver->pfnOpen = DTEDDataset::Open;
999 1293 : poDriver->pfnIdentify = DTEDDataset::Identify;
1000 1293 : poDriver->pfnCreateCopy = DTEDCreateCopy;
1001 :
1002 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1003 : }
|
