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