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