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