Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Vertical Datum Transformation
4 : * Purpose: Implementation of NOAA .gtx vertical datum shift file format.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2010, Frank Warmerdam
9 : * Copyright (c) 2010-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_string.h"
15 : #include "gdal_frmts.h"
16 : #include "gdal_priv.h"
17 : #include "ogr_srs_api.h"
18 : #include "rawdataset.h"
19 :
20 : #include <algorithm>
21 : #include <limits>
22 :
23 : /**
24 :
25 : NOAA .GTX Vertical Datum Grid Shift Format
26 :
27 : All values are bigendian
28 :
29 : Header
30 : ------
31 :
32 : float64 latitude_of_origin
33 : float64 longitude_of_origin (0-360)
34 : float64 cell size (x?y?)
35 : float64 cell size (x?y?)
36 : int32 length in pixels
37 : int32 width in pixels
38 :
39 : Data
40 : ----
41 :
42 : float32 * width in pixels * length in pixels
43 :
44 : Values are an offset in meters between two vertical datums.
45 :
46 : **/
47 :
48 : /************************************************************************/
49 : /* ==================================================================== */
50 : /* GTXDataset */
51 : /* ==================================================================== */
52 : /************************************************************************/
53 :
54 : class GTXDataset final : public RawDataset
55 : {
56 : VSILFILE *fpImage = nullptr; // image data file.
57 :
58 : OGRSpatialReference m_oSRS{};
59 : GDALGeoTransform m_gt{};
60 :
61 : CPL_DISALLOW_COPY_ASSIGN(GTXDataset)
62 :
63 : CPLErr Close() override;
64 :
65 : public:
66 47 : GTXDataset()
67 47 : {
68 47 : m_oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
69 47 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
70 47 : }
71 :
72 : ~GTXDataset() override;
73 :
74 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
75 : CPLErr SetGeoTransform(const GDALGeoTransform >) override;
76 :
77 1 : const OGRSpatialReference *GetSpatialRef() const override
78 : {
79 1 : return &m_oSRS;
80 : }
81 :
82 : static GDALDataset *Open(GDALOpenInfo *);
83 : static int Identify(GDALOpenInfo *);
84 : static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
85 : int nBands, GDALDataType eType,
86 : char **papszOptions);
87 : };
88 :
89 : /************************************************************************/
90 : /* ==================================================================== */
91 : /* GTXRasterBand */
92 : /* ==================================================================== */
93 : /************************************************************************/
94 :
95 : class GTXRasterBand final : public RawRasterBand
96 : {
97 : CPL_DISALLOW_COPY_ASSIGN(GTXRasterBand)
98 :
99 : public:
100 : GTXRasterBand(GDALDataset *poDS, int nBand, VSILFILE *fpRaw,
101 : vsi_l_offset nImgOffset, int nPixelOffset, int nLineOffset,
102 : GDALDataType eDataType, int bNativeOrder);
103 :
104 : ~GTXRasterBand() override;
105 :
106 : double GetNoDataValue(int *pbSuccess = nullptr) override;
107 : };
108 :
109 : /************************************************************************/
110 : /* GTXRasterBand() */
111 : /************************************************************************/
112 :
113 47 : GTXRasterBand::GTXRasterBand(GDALDataset *poDSIn, int nBandIn,
114 : VSILFILE *fpRawIn, vsi_l_offset nImgOffsetIn,
115 : int nPixelOffsetIn, int nLineOffsetIn,
116 47 : GDALDataType eDataTypeIn, int bNativeOrderIn)
117 : : RawRasterBand(poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
118 : nLineOffsetIn, eDataTypeIn, bNativeOrderIn,
119 47 : RawRasterBand::OwnFP::NO)
120 : {
121 47 : }
122 :
123 : /************************************************************************/
124 : /* ~GTXRasterBand() */
125 : /************************************************************************/
126 :
127 94 : GTXRasterBand::~GTXRasterBand()
128 : {
129 94 : }
130 :
131 : /************************************************************************/
132 : /* GetNoDataValue() */
133 : /************************************************************************/
134 :
135 0 : double GTXRasterBand::GetNoDataValue(int *pbSuccess)
136 : {
137 0 : if (pbSuccess)
138 0 : *pbSuccess = TRUE;
139 0 : int bSuccess = FALSE;
140 0 : double dfNoData = GDALPamRasterBand::GetNoDataValue(&bSuccess);
141 0 : if (bSuccess)
142 : {
143 0 : return dfNoData;
144 : }
145 0 : return -88.8888;
146 : }
147 :
148 : /************************************************************************/
149 : /* ==================================================================== */
150 : /* GTXDataset */
151 : /* ==================================================================== */
152 : /************************************************************************/
153 :
154 : /************************************************************************/
155 : /* ~GTXDataset() */
156 : /************************************************************************/
157 :
158 94 : GTXDataset::~GTXDataset()
159 :
160 : {
161 47 : GTXDataset::Close();
162 94 : }
163 :
164 : /************************************************************************/
165 : /* Close() */
166 : /************************************************************************/
167 :
168 85 : CPLErr GTXDataset::Close()
169 : {
170 85 : CPLErr eErr = CE_None;
171 85 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
172 : {
173 47 : if (GTXDataset::FlushCache(true) != CE_None)
174 9 : eErr = CE_Failure;
175 :
176 47 : if (fpImage)
177 : {
178 47 : if (VSIFCloseL(fpImage) != 0)
179 : {
180 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
181 0 : eErr = CE_Failure;
182 : }
183 : }
184 :
185 47 : if (GDALPamDataset::Close() != CE_None)
186 0 : eErr = CE_Failure;
187 : }
188 85 : return eErr;
189 : }
190 :
191 : /************************************************************************/
192 : /* Identify() */
193 : /************************************************************************/
194 :
195 58623 : int GTXDataset::Identify(GDALOpenInfo *poOpenInfo)
196 :
197 : {
198 58623 : if (poOpenInfo->nHeaderBytes < 40)
199 54855 : return FALSE;
200 :
201 3768 : if (!poOpenInfo->IsExtensionEqualToCI("gtx"))
202 3675 : return FALSE;
203 :
204 94 : return TRUE;
205 : }
206 :
207 : /************************************************************************/
208 : /* Open() */
209 : /************************************************************************/
210 :
211 47 : GDALDataset *GTXDataset::Open(GDALOpenInfo *poOpenInfo)
212 :
213 : {
214 47 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
215 0 : return nullptr;
216 :
217 : /* -------------------------------------------------------------------- */
218 : /* Create a corresponding GDALDataset. */
219 : /* -------------------------------------------------------------------- */
220 94 : auto poDS = std::make_unique<GTXDataset>();
221 :
222 47 : poDS->eAccess = poOpenInfo->eAccess;
223 47 : std::swap(poDS->fpImage, poOpenInfo->fpL);
224 :
225 : /* -------------------------------------------------------------------- */
226 : /* Read the header. */
227 : /* -------------------------------------------------------------------- */
228 47 : double gt[6] = {0};
229 47 : CPL_IGNORE_RET_VAL(VSIFReadL(>[3], 8, 1, poDS->fpImage));
230 47 : CPL_IGNORE_RET_VAL(VSIFReadL(>[0], 8, 1, poDS->fpImage));
231 47 : CPL_IGNORE_RET_VAL(VSIFReadL(>[5], 8, 1, poDS->fpImage));
232 47 : CPL_IGNORE_RET_VAL(VSIFReadL(>[1], 8, 1, poDS->fpImage));
233 :
234 47 : CPL_IGNORE_RET_VAL(VSIFReadL(&(poDS->nRasterYSize), 4, 1, poDS->fpImage));
235 47 : CPL_IGNORE_RET_VAL(VSIFReadL(&(poDS->nRasterXSize), 4, 1, poDS->fpImage));
236 :
237 47 : CPL_MSBPTR32(&(poDS->nRasterYSize));
238 47 : CPL_MSBPTR32(&(poDS->nRasterXSize));
239 :
240 47 : CPL_MSBPTR64(>[0]);
241 47 : CPL_MSBPTR64(>[1]);
242 47 : CPL_MSBPTR64(>[3]);
243 47 : CPL_MSBPTR64(>[5]);
244 :
245 47 : poDS->m_gt = GDALGeoTransform(gt);
246 47 : poDS->m_gt[3] +=
247 47 : poDS->m_gt[5] * (static_cast<double>(poDS->nRasterYSize) - 1);
248 :
249 47 : poDS->m_gt[0] -= poDS->m_gt[1] * 0.5;
250 47 : poDS->m_gt[3] += poDS->m_gt[5] * 0.5;
251 :
252 47 : poDS->m_gt[5] *= -1;
253 :
254 47 : if (CPLFetchBool(poOpenInfo->papszOpenOptions,
255 : "SHIFT_ORIGIN_IN_MINUS_180_PLUS_180", false))
256 : {
257 0 : if (poDS->m_gt[0] < -180.0 - poDS->m_gt[1])
258 0 : poDS->m_gt[0] += 360.0;
259 0 : else if (poDS->m_gt[0] > 180.0)
260 0 : poDS->m_gt[0] -= 360.0;
261 : }
262 :
263 94 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
264 47 : static_cast<vsi_l_offset>(poDS->nRasterXSize) * poDS->nRasterYSize >
265 47 : std::numeric_limits<vsi_l_offset>::max() / sizeof(double))
266 : {
267 0 : return nullptr;
268 : }
269 :
270 : /* -------------------------------------------------------------------- */
271 : /* Guess the data type. Since October 1, 2009, it should be */
272 : /* Float32. Before it was double. */
273 : /* -------------------------------------------------------------------- */
274 47 : CPL_IGNORE_RET_VAL(VSIFSeekL(poDS->fpImage, 0, SEEK_END));
275 47 : const vsi_l_offset nSize = VSIFTellL(poDS->fpImage);
276 :
277 47 : GDALDataType eDT = GDT_Float32;
278 94 : if (nSize - 40 == sizeof(double) *
279 47 : static_cast<vsi_l_offset>(poDS->nRasterXSize) *
280 47 : poDS->nRasterYSize)
281 0 : eDT = GDT_Float64;
282 47 : const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
283 47 : if (nDTSize <= 0 || poDS->nRasterXSize > INT_MAX / nDTSize)
284 : {
285 0 : return nullptr;
286 : }
287 :
288 : /* -------------------------------------------------------------------- */
289 : /* Create band information object. */
290 : /* -------------------------------------------------------------------- */
291 : auto poBand = std::make_unique<GTXRasterBand>(
292 47 : poDS.get(), 1, poDS->fpImage,
293 47 : static_cast<vsi_l_offset>(poDS->nRasterYSize - 1) * poDS->nRasterXSize *
294 47 : nDTSize +
295 : 40,
296 141 : nDTSize, poDS->nRasterXSize * -nDTSize, eDT, !CPL_IS_LSB);
297 47 : if (!poBand->IsValid())
298 0 : return nullptr;
299 47 : poDS->SetBand(1, std::move(poBand));
300 :
301 : /* -------------------------------------------------------------------- */
302 : /* Initialize any PAM information. */
303 : /* -------------------------------------------------------------------- */
304 47 : poDS->SetDescription(poOpenInfo->pszFilename);
305 47 : poDS->TryLoadXML();
306 :
307 : /* -------------------------------------------------------------------- */
308 : /* Check for overviews. */
309 : /* -------------------------------------------------------------------- */
310 47 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
311 :
312 47 : return poDS.release();
313 : }
314 :
315 : /************************************************************************/
316 : /* GetGeoTransform() */
317 : /************************************************************************/
318 :
319 3 : CPLErr GTXDataset::GetGeoTransform(GDALGeoTransform >) const
320 :
321 : {
322 3 : gt = m_gt;
323 3 : return CE_None;
324 : }
325 :
326 : /************************************************************************/
327 : /* SetGeoTransform() */
328 : /************************************************************************/
329 :
330 33 : CPLErr GTXDataset::SetGeoTransform(const GDALGeoTransform >)
331 : {
332 33 : if (gt[2] != 0.0 || gt[4] != 0.0)
333 : {
334 0 : CPLError(CE_Failure, CPLE_AppDefined,
335 : "Attempt to write skewed or rotated geotransform to gtx.");
336 0 : return CE_Failure;
337 : }
338 :
339 33 : m_gt = gt;
340 :
341 33 : const double dfXOrigin = m_gt[0] + 0.5 * m_gt[1];
342 33 : const double dfYOrigin = m_gt[3] + (nRasterYSize - 0.5) * m_gt[5];
343 33 : const double dfWidth = m_gt[1];
344 33 : const double dfHeight = -m_gt[5];
345 :
346 33 : unsigned char header[32] = {'\0'};
347 33 : memcpy(header + 0, &dfYOrigin, 8);
348 33 : CPL_MSBPTR64(header + 0);
349 :
350 33 : memcpy(header + 8, &dfXOrigin, 8);
351 33 : CPL_MSBPTR64(header + 8);
352 :
353 33 : memcpy(header + 16, &dfHeight, 8);
354 33 : CPL_MSBPTR64(header + 16);
355 :
356 33 : memcpy(header + 24, &dfWidth, 8);
357 33 : CPL_MSBPTR64(header + 24);
358 :
359 66 : if (VSIFSeekL(fpImage, SEEK_SET, 0) != 0 ||
360 33 : VSIFWriteL(header, 32, 1, fpImage) != 1)
361 : {
362 0 : CPLError(CE_Failure, CPLE_AppDefined,
363 : "Attempt to write geotransform header to GTX failed.");
364 0 : return CE_Failure;
365 : }
366 :
367 33 : return CE_None;
368 : }
369 :
370 : /************************************************************************/
371 : /* Create() */
372 : /************************************************************************/
373 :
374 83 : GDALDataset *GTXDataset::Create(const char *pszFilename, int nXSize, int nYSize,
375 : int /* nBands */, GDALDataType eType,
376 : char ** /* papszOptions */)
377 : {
378 83 : if (eType != GDT_Float32)
379 : {
380 46 : CPLError(CE_Failure, CPLE_AppDefined,
381 : "Attempt to create gtx file with unsupported data type '%s'.",
382 : GDALGetDataTypeName(eType));
383 46 : return nullptr;
384 : }
385 :
386 37 : if (!EQUAL(CPLGetExtensionSafe(pszFilename).c_str(), "gtx"))
387 : {
388 0 : CPLError(CE_Failure, CPLE_AppDefined,
389 : "Attempt to create gtx file with extension other than gtx.");
390 0 : return nullptr;
391 : }
392 :
393 : /* -------------------------------------------------------------------- */
394 : /* Try to create the file. */
395 : /* -------------------------------------------------------------------- */
396 37 : VSILFILE *fp = VSIFOpenL(pszFilename, "wb");
397 37 : if (fp == nullptr)
398 : {
399 3 : CPLError(CE_Failure, CPLE_OpenFailed,
400 : "Attempt to create file `%s' failed.\n", pszFilename);
401 3 : return nullptr;
402 : }
403 :
404 : /* -------------------------------------------------------------------- */
405 : /* Write out the header with stub georeferencing. */
406 : /* -------------------------------------------------------------------- */
407 :
408 34 : unsigned char header[40] = {'\0'};
409 34 : double dfYOrigin = 0.0;
410 34 : memcpy(header + 0, &dfYOrigin, 8);
411 34 : CPL_MSBPTR64(header + 0);
412 :
413 34 : double dfXOrigin = 0.0;
414 34 : memcpy(header + 8, &dfXOrigin, 8);
415 34 : CPL_MSBPTR64(header + 8);
416 :
417 34 : double dfYSize = 0.01;
418 34 : memcpy(header + 16, &dfYSize, 8);
419 34 : CPL_MSBPTR64(header + 16);
420 :
421 34 : double dfXSize = 0.01;
422 34 : memcpy(header + 24, &dfXSize, 8);
423 34 : CPL_MSBPTR64(header + 24);
424 :
425 34 : GInt32 nYSize32 = nYSize;
426 34 : memcpy(header + 32, &nYSize32, 4);
427 34 : CPL_MSBPTR32(header + 32);
428 :
429 34 : GInt32 nXSize32 = nXSize;
430 34 : memcpy(header + 36, &nXSize32, 4);
431 34 : CPL_MSBPTR32(header + 36);
432 :
433 34 : CPL_IGNORE_RET_VAL(VSIFWriteL(header, 40, 1, fp));
434 34 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
435 :
436 34 : return GDALDataset::FromHandle(GDALOpen(pszFilename, GA_Update));
437 : }
438 :
439 : /************************************************************************/
440 : /* GDALRegister_GTX() */
441 : /************************************************************************/
442 :
443 2033 : void GDALRegister_GTX()
444 :
445 : {
446 2033 : if (GDALGetDriverByName("GTX") != nullptr)
447 283 : return;
448 :
449 1750 : GDALDriver *poDriver = new GDALDriver();
450 :
451 1750 : poDriver->SetDescription("GTX");
452 1750 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
453 1750 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "NOAA Vertical Datum .GTX");
454 1750 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "gtx");
455 1750 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
456 : // poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
457 : // "frmt_various.html#GTX" );
458 1750 : poDriver->SetMetadataItem(
459 : GDAL_DMD_OPENOPTIONLIST,
460 : "<OpenOptionList>"
461 : " <Option name='SHIFT_ORIGIN_IN_MINUS_180_PLUS_180' type='boolean' "
462 : "description='Whether to apply a +/-360 deg shift to the longitude of "
463 : "the top left corner so that it is in the [-180,180] range' "
464 : "default='NO'/>"
465 1750 : "</OpenOptionList>");
466 :
467 1750 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Float32");
468 :
469 1750 : poDriver->pfnOpen = GTXDataset::Open;
470 1750 : poDriver->pfnIdentify = GTXDataset::Identify;
471 1750 : poDriver->pfnCreate = GTXDataset::Create;
472 :
473 1750 : GetGDALDriverManager()->RegisterDriver(poDriver);
474 : }
|
