Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Horizontal Datum Formats
4 : * Purpose: Implementation of NTv2 datum shift format used in Canada, France,
5 : * Australia and elsewhere.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : * Financial Support: i-cubed (http://www.i-cubed.com)
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2010, Frank Warmerdam
11 : * Copyright (c) 2010-2012, Even Rouault <even dot rouault at spatialys.com>
12 : *
13 : * SPDX-License-Identifier: MIT
14 : ****************************************************************************/
15 :
16 : // TODO(schwehr): There are a lot of magic numbers in this driver that should
17 : // be changed to constants and documented.
18 :
19 : #include "cpl_string.h"
20 : #include "gdal_frmts.h"
21 : #include "ogr_srs_api.h"
22 : #include "rawdataset.h"
23 :
24 : #include <algorithm>
25 :
26 : // Format documentation: https://github.com/Esri/ntv2-file-routines
27 : // Original archived specification:
28 : // https://web.archive.org/web/20091227232322/http://www.mgs.gov.on.ca/stdprodconsume/groups/content/@mgs/@iandit/documents/resourcelist/stel02_047447.pdf
29 :
30 : /**
31 : * The header for the file, and each grid consists of 11 16byte records.
32 : * The first half is an ASCII label, and the second half is the value
33 : * often in a little endian int or float.
34 : *
35 : * Example:
36 :
37 : 00000000 4e 55 4d 5f 4f 52 45 43 0b 00 00 00 00 00 00 00 |NUM_OREC........|
38 : 00000010 4e 55 4d 5f 53 52 45 43 0b 00 00 00 00 00 00 00 |NUM_SREC........|
39 : 00000020 4e 55 4d 5f 46 49 4c 45 01 00 00 00 00 00 00 00 |NUM_FILE........|
40 : 00000030 47 53 5f 54 59 50 45 20 53 45 43 4f 4e 44 53 20 |GS_TYPE SECONDS |
41 : 00000040 56 45 52 53 49 4f 4e 20 49 47 4e 30 37 5f 30 31 |VERSION IGN07_01|
42 : 00000050 53 59 53 54 45 4d 5f 46 4e 54 46 20 20 20 20 20 |SYSTEM_FNTF |
43 : 00000060 53 59 53 54 45 4d 5f 54 52 47 46 39 33 20 20 20 |SYSTEM_TRGF93 |
44 : 00000070 4d 41 4a 4f 52 5f 46 20 cd cc cc 4c c2 54 58 41 |MAJOR_F ...L.TXA|
45 : 00000080 4d 49 4e 4f 52 5f 46 20 00 00 00 c0 88 3f 58 41 |MINOR_F .....?XA|
46 : 00000090 4d 41 4a 4f 52 5f 54 20 00 00 00 40 a6 54 58 41 |MAJOR_T ...@.TXA|
47 : 000000a0 4d 49 4e 4f 52 5f 54 20 27 e0 1a 14 c4 3f 58 41 |MINOR_T '....?XA|
48 : 000000b0 53 55 42 5f 4e 41 4d 45 46 52 41 4e 43 45 20 20 |SUB_NAMEFRANCE |
49 : 000000c0 50 41 52 45 4e 54 20 20 4e 4f 4e 45 20 20 20 20 |PARENT NONE |
50 : 000000d0 43 52 45 41 54 45 44 20 33 31 2f 31 30 2f 30 37 |CREATED 31/10/07|
51 : 000000e0 55 50 44 41 54 45 44 20 20 20 20 20 20 20 20 20 |UPDATED |
52 : 000000f0 53 5f 4c 41 54 20 20 20 00 00 00 00 80 04 02 41 |S_LAT .......A|
53 : 00000100 4e 5f 4c 41 54 20 20 20 00 00 00 00 00 da 06 41 |N_LAT .......A|
54 : 00000110 45 5f 4c 4f 4e 47 20 20 00 00 00 00 00 94 e1 c0 |E_LONG ........|
55 : 00000120 57 5f 4c 4f 4e 47 20 20 00 00 00 00 00 56 d3 40 |W_LONG .....V.@|
56 : 00000130 4c 41 54 5f 49 4e 43 20 00 00 00 00 00 80 76 40 |LAT_INC ......v@|
57 : 00000140 4c 4f 4e 47 5f 49 4e 43 00 00 00 00 00 80 76 40 |LONG_INC......v@|
58 : 00000150 47 53 5f 43 4f 55 4e 54 a4 43 00 00 00 00 00 00 |GS_COUNT.C......|
59 : 00000160 94 f7 c1 3e 70 ee a3 3f 2a c7 84 3d ff 42 af 3d |...>p..?*..=.B.=|
60 :
61 : the actual grid data is a raster with 4 float32 bands (lat offset, long
62 : offset, lat error, long error). The offset values are in arc seconds.
63 : The grid is flipped in the x and y axis from our usual GDAL orientation.
64 : That is, the first pixel is the south east corner with scanlines going
65 : east to west, and rows from south to north. As a GDAL dataset we represent
66 : these both in the more conventional orientation.
67 : */
68 :
69 : constexpr size_t knREGULAR_RECORD_SIZE = 16;
70 : // This one is for velocity grids such as the NAD83(CRSR)v7 / NAD83v70VG.gvb
71 : // which is the only example I know actually of that format variant.
72 : constexpr size_t knMAX_RECORD_SIZE = 24;
73 :
74 : /************************************************************************/
75 : /* ==================================================================== */
76 : /* NTv2Dataset */
77 : /* ==================================================================== */
78 : /************************************************************************/
79 :
80 : class NTv2Dataset final : public RawDataset
81 : {
82 : public:
83 : RawRasterBand::ByteOrder m_eByteOrder = RawRasterBand::NATIVE_BYTE_ORDER;
84 : bool m_bMustSwap = false;
85 : VSILFILE *fpImage = nullptr; // image data file.
86 :
87 : size_t nRecordSize = 0;
88 : vsi_l_offset nGridOffset = 0;
89 :
90 : OGRSpatialReference m_oSRS{};
91 : double adfGeoTransform[6];
92 :
93 : void CaptureMetadataItem(const char *pszItem);
94 :
95 : bool OpenGrid(const char *pachGridHeader, vsi_l_offset nDataStart);
96 :
97 : CPL_DISALLOW_COPY_ASSIGN(NTv2Dataset)
98 :
99 : CPLErr Close() override;
100 :
101 : public:
102 : NTv2Dataset();
103 : ~NTv2Dataset() override;
104 :
105 : CPLErr SetGeoTransform(double *padfTransform) override;
106 : CPLErr GetGeoTransform(double *padfTransform) override;
107 :
108 6 : const OGRSpatialReference *GetSpatialRef() const override
109 : {
110 6 : return &m_oSRS;
111 : }
112 :
113 : CPLErr FlushCache(bool bAtClosing) override;
114 :
115 : static GDALDataset *Open(GDALOpenInfo *);
116 : static int Identify(GDALOpenInfo *);
117 : static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
118 : int nBandsIn, GDALDataType eType,
119 : char **papszOptions);
120 : };
121 :
122 : /************************************************************************/
123 : /* ==================================================================== */
124 : /* NTv2Dataset */
125 : /* ==================================================================== */
126 : /************************************************************************/
127 :
128 : /************************************************************************/
129 : /* NTv2Dataset() */
130 : /************************************************************************/
131 :
132 22 : NTv2Dataset::NTv2Dataset()
133 : {
134 22 : m_oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
135 22 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
136 :
137 22 : adfGeoTransform[0] = 0.0;
138 22 : adfGeoTransform[1] = 0.0; // TODO(schwehr): Should this be 1.0?
139 22 : adfGeoTransform[2] = 0.0;
140 22 : adfGeoTransform[3] = 0.0;
141 22 : adfGeoTransform[4] = 0.0;
142 22 : adfGeoTransform[5] = 0.0; // TODO(schwehr): Should this be 1.0?
143 22 : }
144 :
145 : /************************************************************************/
146 : /* ~NTv2Dataset() */
147 : /************************************************************************/
148 :
149 44 : NTv2Dataset::~NTv2Dataset()
150 :
151 : {
152 22 : NTv2Dataset::Close();
153 44 : }
154 :
155 : /************************************************************************/
156 : /* Close() */
157 : /************************************************************************/
158 :
159 44 : CPLErr NTv2Dataset::Close()
160 : {
161 44 : CPLErr eErr = CE_None;
162 44 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
163 : {
164 22 : if (NTv2Dataset::FlushCache(true) != CE_None)
165 0 : eErr = CE_Failure;
166 :
167 22 : if (fpImage)
168 : {
169 22 : if (VSIFCloseL(fpImage) != 0)
170 : {
171 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
172 0 : eErr = CE_Failure;
173 : }
174 : }
175 :
176 22 : if (GDALPamDataset::Close() != CE_None)
177 0 : eErr = CE_Failure;
178 : }
179 44 : return eErr;
180 : }
181 :
182 : /************************************************************************/
183 : /* SwapPtr32IfNecessary() */
184 : /************************************************************************/
185 :
186 78 : static void SwapPtr32IfNecessary(bool bMustSwap, void *ptr)
187 : {
188 78 : if (bMustSwap)
189 : {
190 35 : CPL_SWAP32PTR(static_cast<GByte *>(ptr));
191 : }
192 78 : }
193 :
194 : /************************************************************************/
195 : /* SwapPtr64IfNecessary() */
196 : /************************************************************************/
197 :
198 326 : static void SwapPtr64IfNecessary(bool bMustSwap, void *ptr)
199 : {
200 326 : if (bMustSwap)
201 : {
202 143 : CPL_SWAP64PTR(static_cast<GByte *>(ptr));
203 : }
204 326 : }
205 :
206 : /************************************************************************/
207 : /* FlushCache() */
208 : /************************************************************************/
209 :
210 23 : CPLErr NTv2Dataset::FlushCache(bool bAtClosing)
211 :
212 : {
213 : /* -------------------------------------------------------------------- */
214 : /* Nothing to do in readonly mode, or if nothing seems to have */
215 : /* changed metadata wise. */
216 : /* -------------------------------------------------------------------- */
217 23 : if (eAccess != GA_Update || !(GetPamFlags() & GPF_DIRTY))
218 : {
219 19 : return RawDataset::FlushCache(bAtClosing);
220 : }
221 :
222 : /* -------------------------------------------------------------------- */
223 : /* Load grid and file headers. */
224 : /* -------------------------------------------------------------------- */
225 4 : const int nRecords = 11;
226 4 : char achFileHeader[nRecords * knMAX_RECORD_SIZE] = {'\0'};
227 4 : char achGridHeader[nRecords * knMAX_RECORD_SIZE] = {'\0'};
228 :
229 4 : bool bOK = VSIFSeekL(fpImage, 0, SEEK_SET) == 0;
230 4 : bOK &=
231 4 : VSIFReadL(achFileHeader, nRecords, nRecordSize, fpImage) == nRecordSize;
232 :
233 4 : bOK &= VSIFSeekL(fpImage, nGridOffset, SEEK_SET) == 0;
234 4 : bOK &=
235 4 : VSIFReadL(achGridHeader, nRecords, nRecordSize, fpImage) == nRecordSize;
236 :
237 : /* -------------------------------------------------------------------- */
238 : /* Update the grid, and file headers with any available */
239 : /* metadata. If all available metadata is recognised then mark */
240 : /* things "clean" from a PAM point of view. */
241 : /* -------------------------------------------------------------------- */
242 4 : char **papszMD = GetMetadata();
243 4 : bool bSomeLeftOver = false;
244 :
245 52 : for (int i = 0; papszMD != nullptr && papszMD[i] != nullptr; i++)
246 : {
247 48 : const size_t nMinLen = 8;
248 48 : char *pszKey = nullptr;
249 48 : const char *pszValue = CPLParseNameValue(papszMD[i], &pszKey);
250 48 : if (pszKey == nullptr)
251 0 : continue;
252 :
253 48 : if (EQUAL(pszKey, "GS_TYPE"))
254 : {
255 4 : memcpy(achFileHeader + 3 * nRecordSize + 8, " ", 8);
256 4 : memcpy(achFileHeader + 3 * nRecordSize + 8, pszValue,
257 4 : std::min(nMinLen, strlen(pszValue)));
258 : }
259 44 : else if (EQUAL(pszKey, "VERSION"))
260 : {
261 4 : memcpy(achFileHeader + 4 * nRecordSize + 8, " ", 8);
262 4 : memcpy(achFileHeader + 4 * nRecordSize + 8, pszValue,
263 4 : std::min(nMinLen, strlen(pszValue)));
264 : }
265 40 : else if (EQUAL(pszKey, "SYSTEM_F"))
266 : {
267 4 : memcpy(achFileHeader + 5 * nRecordSize + 8, " ", 8);
268 4 : memcpy(achFileHeader + 5 * nRecordSize + 8, pszValue,
269 4 : std::min(nMinLen, strlen(pszValue)));
270 : }
271 36 : else if (EQUAL(pszKey, "SYSTEM_T"))
272 : {
273 4 : memcpy(achFileHeader + 6 * nRecordSize + 8, " ", 8);
274 4 : memcpy(achFileHeader + 6 * nRecordSize + 8, pszValue,
275 4 : std::min(nMinLen, strlen(pszValue)));
276 : }
277 32 : else if (EQUAL(pszKey, "MAJOR_F"))
278 : {
279 4 : double dfValue = CPLAtof(pszValue);
280 4 : SwapPtr64IfNecessary(m_bMustSwap, &dfValue);
281 4 : memcpy(achFileHeader + 7 * nRecordSize + 8, &dfValue, 8);
282 : }
283 28 : else if (EQUAL(pszKey, "MINOR_F"))
284 : {
285 4 : double dfValue = CPLAtof(pszValue);
286 4 : SwapPtr64IfNecessary(m_bMustSwap, &dfValue);
287 4 : memcpy(achFileHeader + 8 * nRecordSize + 8, &dfValue, 8);
288 : }
289 24 : else if (EQUAL(pszKey, "MAJOR_T"))
290 : {
291 4 : double dfValue = CPLAtof(pszValue);
292 4 : SwapPtr64IfNecessary(m_bMustSwap, &dfValue);
293 4 : memcpy(achFileHeader + 9 * nRecordSize + 8, &dfValue, 8);
294 : }
295 20 : else if (EQUAL(pszKey, "MINOR_T"))
296 : {
297 4 : double dfValue = CPLAtof(pszValue);
298 4 : SwapPtr64IfNecessary(m_bMustSwap, &dfValue);
299 4 : memcpy(achFileHeader + 10 * nRecordSize + 8, &dfValue, 8);
300 : }
301 16 : else if (EQUAL(pszKey, "SUB_NAME"))
302 : {
303 4 : memcpy(achGridHeader + /*0*nRecordSize+*/ 8, " ", 8);
304 4 : memcpy(achGridHeader + /*0*nRecordSize+*/ 8, pszValue,
305 4 : std::min(nMinLen, strlen(pszValue)));
306 : }
307 12 : else if (EQUAL(pszKey, "PARENT"))
308 : {
309 4 : memcpy(achGridHeader + 1 * nRecordSize + 8, " ", 8);
310 4 : memcpy(achGridHeader + 1 * nRecordSize + 8, pszValue,
311 4 : std::min(nMinLen, strlen(pszValue)));
312 : }
313 8 : else if (EQUAL(pszKey, "CREATED"))
314 : {
315 4 : memcpy(achGridHeader + 2 * nRecordSize + 8, " ", 8);
316 4 : memcpy(achGridHeader + 2 * nRecordSize + 8, pszValue,
317 4 : std::min(nMinLen, strlen(pszValue)));
318 : }
319 4 : else if (EQUAL(pszKey, "UPDATED"))
320 : {
321 4 : memcpy(achGridHeader + 3 * nRecordSize + 8, " ", 8);
322 4 : memcpy(achGridHeader + 3 * nRecordSize + 8, pszValue,
323 4 : std::min(nMinLen, strlen(pszValue)));
324 : }
325 : else
326 : {
327 0 : bSomeLeftOver = true;
328 : }
329 :
330 48 : CPLFree(pszKey);
331 : }
332 :
333 : /* -------------------------------------------------------------------- */
334 : /* Load grid and file headers. */
335 : /* -------------------------------------------------------------------- */
336 4 : bOK &= VSIFSeekL(fpImage, 0, SEEK_SET) == 0;
337 4 : bOK &= VSIFWriteL(achFileHeader, nRecords, nRecordSize, fpImage) ==
338 4 : nRecordSize;
339 :
340 4 : bOK &= VSIFSeekL(fpImage, nGridOffset, SEEK_SET) == 0;
341 4 : bOK &= VSIFWriteL(achGridHeader, nRecords, nRecordSize, fpImage) ==
342 4 : nRecordSize;
343 :
344 : /* -------------------------------------------------------------------- */
345 : /* Clear flags if we got everything, then let pam and below do */
346 : /* their flushing. */
347 : /* -------------------------------------------------------------------- */
348 4 : if (!bSomeLeftOver)
349 4 : SetPamFlags(GetPamFlags() & (~GPF_DIRTY));
350 :
351 4 : if (RawDataset::FlushCache(bAtClosing) != CE_None)
352 0 : bOK = false;
353 4 : return bOK ? CE_None : CE_Failure;
354 : }
355 :
356 : /************************************************************************/
357 : /* Identify() */
358 : /************************************************************************/
359 :
360 51483 : int NTv2Dataset::Identify(GDALOpenInfo *poOpenInfo)
361 :
362 : {
363 51483 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "NTv2:"))
364 8 : return TRUE;
365 :
366 51475 : if (poOpenInfo->nHeaderBytes < 64)
367 48068 : return FALSE;
368 :
369 3407 : const char *pszHeader =
370 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
371 3407 : if (!STARTS_WITH_CI(pszHeader + 0, "NUM_OREC"))
372 3369 : return FALSE;
373 :
374 38 : if (!STARTS_WITH_CI(pszHeader + knREGULAR_RECORD_SIZE, "NUM_SREC") &&
375 0 : !STARTS_WITH_CI(pszHeader + knMAX_RECORD_SIZE, "NUM_SREC"))
376 0 : return FALSE;
377 :
378 38 : return TRUE;
379 : }
380 :
381 : /************************************************************************/
382 : /* Open() */
383 : /************************************************************************/
384 :
385 22 : GDALDataset *NTv2Dataset::Open(GDALOpenInfo *poOpenInfo)
386 :
387 : {
388 22 : if (!Identify(poOpenInfo))
389 0 : return nullptr;
390 :
391 : /* -------------------------------------------------------------------- */
392 : /* Are we targeting a particular grid? */
393 : /* -------------------------------------------------------------------- */
394 44 : CPLString osFilename;
395 22 : int iTargetGrid = -1;
396 :
397 22 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "NTv2:"))
398 : {
399 4 : const char *pszRest = poOpenInfo->pszFilename + 5;
400 :
401 4 : iTargetGrid = atoi(pszRest);
402 8 : while (*pszRest != '\0' && *pszRest != ':')
403 4 : pszRest++;
404 :
405 4 : if (*pszRest == ':')
406 4 : pszRest++;
407 :
408 4 : osFilename = pszRest;
409 : }
410 : else
411 : {
412 18 : osFilename = poOpenInfo->pszFilename;
413 : }
414 :
415 : /* -------------------------------------------------------------------- */
416 : /* Create a corresponding GDALDataset. */
417 : /* -------------------------------------------------------------------- */
418 44 : auto poDS = std::make_unique<NTv2Dataset>();
419 22 : poDS->eAccess = poOpenInfo->eAccess;
420 :
421 : /* -------------------------------------------------------------------- */
422 : /* Open the file. */
423 : /* -------------------------------------------------------------------- */
424 22 : if (poOpenInfo->eAccess == GA_ReadOnly)
425 16 : poDS->fpImage = VSIFOpenL(osFilename, "rb");
426 : else
427 6 : poDS->fpImage = VSIFOpenL(osFilename, "rb+");
428 :
429 22 : if (poDS->fpImage == nullptr)
430 : {
431 0 : return nullptr;
432 : }
433 :
434 : /* -------------------------------------------------------------------- */
435 : /* Read the file header. */
436 : /* -------------------------------------------------------------------- */
437 22 : char achHeader[11 * knMAX_RECORD_SIZE] = {0};
438 :
439 44 : if (VSIFSeekL(poDS->fpImage, 0, SEEK_SET) != 0 ||
440 22 : VSIFReadL(achHeader, 1, 64, poDS->fpImage) != 64)
441 : {
442 0 : return nullptr;
443 : }
444 :
445 22 : poDS->nRecordSize =
446 22 : STARTS_WITH_CI(achHeader + knMAX_RECORD_SIZE, "NUM_SREC")
447 22 : ? knMAX_RECORD_SIZE
448 : : knREGULAR_RECORD_SIZE;
449 22 : if (VSIFReadL(achHeader + 64, 1, 11 * poDS->nRecordSize - 64,
450 44 : poDS->fpImage) != 11 * poDS->nRecordSize - 64)
451 : {
452 0 : return nullptr;
453 : }
454 :
455 13 : const bool bIsLE = achHeader[8] == 11 && achHeader[9] == 0 &&
456 35 : achHeader[10] == 0 && achHeader[11] == 0;
457 9 : const bool bIsBE = achHeader[8] == 0 && achHeader[9] == 0 &&
458 31 : achHeader[10] == 0 && achHeader[11] == 11;
459 22 : if (!bIsLE && !bIsBE)
460 : {
461 0 : return nullptr;
462 : }
463 22 : poDS->m_eByteOrder = bIsLE ? RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN
464 : : RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN;
465 22 : poDS->m_bMustSwap = poDS->m_eByteOrder != RawRasterBand::NATIVE_BYTE_ORDER;
466 :
467 22 : SwapPtr32IfNecessary(poDS->m_bMustSwap,
468 22 : achHeader + 2 * poDS->nRecordSize + 8);
469 22 : GInt32 nSubFileCount = 0;
470 22 : memcpy(&nSubFileCount, achHeader + 2 * poDS->nRecordSize + 8, 4);
471 22 : if (nSubFileCount <= 0 || nSubFileCount >= 1024)
472 : {
473 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for NUM_FILE : %d",
474 : nSubFileCount);
475 0 : return nullptr;
476 : }
477 :
478 22 : poDS->CaptureMetadataItem(achHeader + 3 * poDS->nRecordSize);
479 22 : poDS->CaptureMetadataItem(achHeader + 4 * poDS->nRecordSize);
480 22 : poDS->CaptureMetadataItem(achHeader + 5 * poDS->nRecordSize);
481 22 : poDS->CaptureMetadataItem(achHeader + 6 * poDS->nRecordSize);
482 :
483 22 : double dfValue = 0.0;
484 22 : memcpy(&dfValue, achHeader + 7 * poDS->nRecordSize + 8, 8);
485 22 : SwapPtr64IfNecessary(poDS->m_bMustSwap, &dfValue);
486 44 : CPLString osFValue;
487 22 : osFValue.Printf("%.15g", dfValue);
488 22 : poDS->SetMetadataItem("MAJOR_F", osFValue);
489 :
490 22 : memcpy(&dfValue, achHeader + 8 * poDS->nRecordSize + 8, 8);
491 22 : SwapPtr64IfNecessary(poDS->m_bMustSwap, &dfValue);
492 22 : osFValue.Printf("%.15g", dfValue);
493 22 : poDS->SetMetadataItem("MINOR_F", osFValue);
494 :
495 22 : memcpy(&dfValue, achHeader + 9 * poDS->nRecordSize + 8, 8);
496 22 : SwapPtr64IfNecessary(poDS->m_bMustSwap, &dfValue);
497 22 : osFValue.Printf("%.15g", dfValue);
498 22 : poDS->SetMetadataItem("MAJOR_T", osFValue);
499 :
500 22 : memcpy(&dfValue, achHeader + 10 * poDS->nRecordSize + 8, 8);
501 22 : SwapPtr64IfNecessary(poDS->m_bMustSwap, &dfValue);
502 22 : osFValue.Printf("%.15g", dfValue);
503 22 : poDS->SetMetadataItem("MINOR_T", osFValue);
504 :
505 : /* ==================================================================== */
506 : /* Loop over grids. */
507 : /* ==================================================================== */
508 22 : vsi_l_offset nGridOffset = 11 * poDS->nRecordSize;
509 :
510 50 : for (int iGrid = 0; iGrid < nSubFileCount; iGrid++)
511 : {
512 56 : if (VSIFSeekL(poDS->fpImage, nGridOffset, SEEK_SET) < 0 ||
513 28 : VSIFReadL(achHeader, 11, poDS->nRecordSize, poDS->fpImage) !=
514 28 : poDS->nRecordSize)
515 : {
516 0 : CPLError(CE_Failure, CPLE_AppDefined,
517 : "Cannot read header for subfile %d", iGrid);
518 0 : return nullptr;
519 : }
520 :
521 196 : for (int i = 4; i <= 9; i++)
522 168 : SwapPtr64IfNecessary(poDS->m_bMustSwap,
523 168 : achHeader + i * poDS->nRecordSize + 8);
524 :
525 28 : SwapPtr32IfNecessary(poDS->m_bMustSwap,
526 28 : achHeader + 10 * poDS->nRecordSize + 8);
527 :
528 28 : GUInt32 nGSCount = 0;
529 28 : memcpy(&nGSCount, achHeader + 10 * poDS->nRecordSize + 8, 4);
530 :
531 28 : CPLString osSubName;
532 28 : osSubName.assign(achHeader + 8, 8);
533 28 : osSubName.Trim();
534 :
535 : // If this is our target grid, open it as a dataset.
536 28 : if (iTargetGrid == iGrid || (iTargetGrid == -1 && iGrid == 0))
537 : {
538 22 : if (!poDS->OpenGrid(achHeader, nGridOffset))
539 : {
540 0 : return nullptr;
541 : }
542 : }
543 :
544 : // If we are opening the file as a whole, list subdatasets.
545 28 : if (iTargetGrid == -1)
546 : {
547 40 : CPLString osKey;
548 40 : CPLString osValue;
549 20 : osKey.Printf("SUBDATASET_%d_NAME", iGrid);
550 20 : osValue.Printf("NTv2:%d:%s", iGrid, osFilename.c_str());
551 20 : poDS->SetMetadataItem(osKey, osValue, "SUBDATASETS");
552 :
553 20 : osKey.Printf("SUBDATASET_%d_DESC", iGrid);
554 20 : osValue.Printf("%s", osSubName.c_str());
555 20 : poDS->SetMetadataItem(osKey, osValue, "SUBDATASETS");
556 : }
557 :
558 28 : nGridOffset +=
559 28 : (11 + static_cast<vsi_l_offset>(nGSCount)) * poDS->nRecordSize;
560 : }
561 :
562 : /* -------------------------------------------------------------------- */
563 : /* Initialize any PAM information. */
564 : /* -------------------------------------------------------------------- */
565 22 : poDS->SetDescription(poOpenInfo->pszFilename);
566 22 : poDS->TryLoadXML();
567 :
568 : /* -------------------------------------------------------------------- */
569 : /* Check for overviews. */
570 : /* -------------------------------------------------------------------- */
571 22 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
572 :
573 22 : return poDS.release();
574 : }
575 :
576 : /************************************************************************/
577 : /* OpenGrid() */
578 : /* */
579 : /* Note that the caller will already have byte swapped needed */
580 : /* portions of the header. */
581 : /************************************************************************/
582 :
583 22 : bool NTv2Dataset::OpenGrid(const char *pachHeader, vsi_l_offset nGridOffsetIn)
584 :
585 : {
586 22 : nGridOffset = nGridOffsetIn;
587 :
588 : /* -------------------------------------------------------------------- */
589 : /* Read the grid header. */
590 : /* -------------------------------------------------------------------- */
591 22 : CaptureMetadataItem(pachHeader + 0 * nRecordSize);
592 22 : CaptureMetadataItem(pachHeader + 1 * nRecordSize);
593 22 : CaptureMetadataItem(pachHeader + 2 * nRecordSize);
594 22 : CaptureMetadataItem(pachHeader + 3 * nRecordSize);
595 :
596 : double s_lat, n_lat, e_long, w_long, lat_inc, long_inc;
597 22 : memcpy(&s_lat, pachHeader + 4 * nRecordSize + 8, 8);
598 22 : memcpy(&n_lat, pachHeader + 5 * nRecordSize + 8, 8);
599 22 : memcpy(&e_long, pachHeader + 6 * nRecordSize + 8, 8);
600 22 : memcpy(&w_long, pachHeader + 7 * nRecordSize + 8, 8);
601 22 : memcpy(&lat_inc, pachHeader + 8 * nRecordSize + 8, 8);
602 22 : memcpy(&long_inc, pachHeader + 9 * nRecordSize + 8, 8);
603 :
604 22 : e_long *= -1;
605 22 : w_long *= -1;
606 :
607 22 : if (long_inc == 0.0 || lat_inc == 0.0)
608 0 : return false;
609 22 : const double dfXSize = floor((e_long - w_long) / long_inc + 1.5);
610 22 : const double dfYSize = floor((n_lat - s_lat) / lat_inc + 1.5);
611 22 : if (!(dfXSize >= 0 && dfXSize < INT_MAX) ||
612 22 : !(dfYSize >= 0 && dfYSize < INT_MAX))
613 0 : return false;
614 22 : nRasterXSize = static_cast<int>(dfXSize);
615 22 : nRasterYSize = static_cast<int>(dfYSize);
616 :
617 22 : const int l_nBands = nRecordSize == knREGULAR_RECORD_SIZE ? 4 : 6;
618 22 : const int nPixelSize = l_nBands * 4;
619 :
620 22 : if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
621 0 : return false;
622 22 : if (nRasterXSize > INT_MAX / nPixelSize)
623 0 : return false;
624 :
625 : /* -------------------------------------------------------------------- */
626 : /* Create band information object. */
627 : /* */
628 : /* We use unusual offsets to remap from bottom to top, to top */
629 : /* to bottom orientation, and also to remap east to west, to */
630 : /* west to east. */
631 : /* -------------------------------------------------------------------- */
632 110 : for (int iBand = 0; iBand < l_nBands; iBand++)
633 : {
634 : auto poBand = RawRasterBand::Create(
635 : this, iBand + 1, fpImage,
636 88 : nGridOffset + 4 * iBand + 11 * nRecordSize +
637 88 : static_cast<vsi_l_offset>(nRasterXSize - 1) * nPixelSize +
638 88 : static_cast<vsi_l_offset>(nRasterYSize - 1) * nPixelSize *
639 88 : nRasterXSize,
640 88 : -nPixelSize, -nPixelSize * nRasterXSize, GDT_Float32, m_eByteOrder,
641 88 : RawRasterBand::OwnFP::NO);
642 88 : if (!poBand)
643 0 : return false;
644 88 : SetBand(iBand + 1, std::move(poBand));
645 : }
646 :
647 22 : if (l_nBands == 4)
648 : {
649 22 : GetRasterBand(1)->SetDescription("Latitude Offset (arc seconds)");
650 22 : GetRasterBand(2)->SetDescription("Longitude Offset (arc seconds)");
651 22 : GetRasterBand(2)->SetMetadataItem("positive_value", "west");
652 22 : GetRasterBand(3)->SetDescription("Latitude Error");
653 22 : GetRasterBand(4)->SetDescription("Longitude Error");
654 : }
655 : else
656 : {
657 : // A bit surprising that the order is easting, northing here, contrary
658 : // to the classic NTv2 order.... Verified on NAD83v70VG.gvb
659 : // (https://webapp.geod.nrcan.gc.ca/geod/process/download-helper.php?file_id=NAD83v70VG)
660 : // against the TRX software
661 : // (https://webapp.geod.nrcan.gc.ca/geod/process/download-helper.php?file_id=trx)
662 : // https://webapp.geod.nrcan.gc.ca/geod/tools-outils/nad83-docs.php
663 : // Unfortunately I couldn't find an official documentation of the format
664 : // !
665 0 : GetRasterBand(1)->SetDescription("East velocity (mm/year)");
666 0 : GetRasterBand(2)->SetDescription("North velocity (mm/year)");
667 0 : GetRasterBand(3)->SetDescription("Up velocity (mm/year)");
668 0 : GetRasterBand(4)->SetDescription("East velocity Error (mm/year)");
669 0 : GetRasterBand(5)->SetDescription("North velocity Error (mm/year)");
670 0 : GetRasterBand(6)->SetDescription("Up velocity Error (mm/year)");
671 : }
672 :
673 : /* -------------------------------------------------------------------- */
674 : /* Setup georeferencing. */
675 : /* -------------------------------------------------------------------- */
676 22 : adfGeoTransform[0] = (w_long - long_inc * 0.5) / 3600.0;
677 22 : adfGeoTransform[1] = long_inc / 3600.0;
678 22 : adfGeoTransform[2] = 0.0;
679 22 : adfGeoTransform[3] = (n_lat + lat_inc * 0.5) / 3600.0;
680 22 : adfGeoTransform[4] = 0.0;
681 22 : adfGeoTransform[5] = (-1 * lat_inc) / 3600.0;
682 :
683 22 : return true;
684 : }
685 :
686 : /************************************************************************/
687 : /* CaptureMetadataItem() */
688 : /************************************************************************/
689 :
690 176 : void NTv2Dataset::CaptureMetadataItem(const char *pszItem)
691 :
692 : {
693 352 : CPLString osKey;
694 352 : CPLString osValue;
695 :
696 176 : osKey.assign(pszItem, 8);
697 176 : osValue.assign(pszItem + 8, 8);
698 :
699 176 : SetMetadataItem(osKey.Trim(), osValue.Trim());
700 176 : }
701 :
702 : /************************************************************************/
703 : /* GetGeoTransform() */
704 : /************************************************************************/
705 :
706 8 : CPLErr NTv2Dataset::GetGeoTransform(double *padfTransform)
707 :
708 : {
709 8 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
710 8 : return CE_None;
711 : }
712 :
713 : /************************************************************************/
714 : /* SetGeoTransform() */
715 : /************************************************************************/
716 :
717 4 : CPLErr NTv2Dataset::SetGeoTransform(double *padfTransform)
718 :
719 : {
720 4 : if (eAccess == GA_ReadOnly)
721 : {
722 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
723 : "Unable to update geotransform on readonly file.");
724 0 : return CE_Failure;
725 : }
726 :
727 4 : if (padfTransform[2] != 0.0 || padfTransform[4] != 0.0)
728 : {
729 0 : CPLError(CE_Failure, CPLE_AppDefined,
730 : "Rotated and sheared geotransforms not supported for NTv2.");
731 0 : return CE_Failure;
732 : }
733 :
734 4 : memcpy(adfGeoTransform, padfTransform, sizeof(double) * 6);
735 :
736 : /* -------------------------------------------------------------------- */
737 : /* Update grid header. */
738 : /* -------------------------------------------------------------------- */
739 4 : char achHeader[11 * knMAX_RECORD_SIZE] = {'\0'};
740 :
741 : // read grid header
742 4 : CPL_IGNORE_RET_VAL(VSIFSeekL(fpImage, nGridOffset, SEEK_SET));
743 4 : CPL_IGNORE_RET_VAL(VSIFReadL(achHeader, 11, nRecordSize, fpImage));
744 :
745 : // S_LAT
746 4 : double dfValue =
747 4 : 3600 * (adfGeoTransform[3] + (nRasterYSize - 0.5) * adfGeoTransform[5]);
748 4 : SwapPtr64IfNecessary(m_bMustSwap, &dfValue);
749 4 : memcpy(achHeader + 4 * nRecordSize + 8, &dfValue, 8);
750 :
751 : // N_LAT
752 4 : dfValue = 3600 * (adfGeoTransform[3] + 0.5 * adfGeoTransform[5]);
753 4 : SwapPtr64IfNecessary(m_bMustSwap, &dfValue);
754 4 : memcpy(achHeader + 5 * nRecordSize + 8, &dfValue, 8);
755 :
756 : // E_LONG
757 4 : dfValue = -3600 *
758 4 : (adfGeoTransform[0] + (nRasterXSize - 0.5) * adfGeoTransform[1]);
759 4 : SwapPtr64IfNecessary(m_bMustSwap, &dfValue);
760 4 : memcpy(achHeader + 6 * nRecordSize + 8, &dfValue, 8);
761 :
762 : // W_LONG
763 4 : dfValue = -3600 * (adfGeoTransform[0] + 0.5 * adfGeoTransform[1]);
764 4 : SwapPtr64IfNecessary(m_bMustSwap, &dfValue);
765 4 : memcpy(achHeader + 7 * nRecordSize + 8, &dfValue, 8);
766 :
767 : // LAT_INC
768 4 : dfValue = -3600 * adfGeoTransform[5];
769 4 : SwapPtr64IfNecessary(m_bMustSwap, &dfValue);
770 4 : memcpy(achHeader + 8 * nRecordSize + 8, &dfValue, 8);
771 :
772 : // LONG_INC
773 4 : dfValue = 3600 * adfGeoTransform[1];
774 4 : SwapPtr64IfNecessary(m_bMustSwap, &dfValue);
775 4 : memcpy(achHeader + 9 * nRecordSize + 8, &dfValue, 8);
776 :
777 : // write grid header.
778 4 : CPL_IGNORE_RET_VAL(VSIFSeekL(fpImage, nGridOffset, SEEK_SET));
779 4 : CPL_IGNORE_RET_VAL(VSIFWriteL(achHeader, 11, nRecordSize, fpImage));
780 :
781 4 : return CE_None;
782 : }
783 :
784 : /************************************************************************/
785 : /* Create() */
786 : /************************************************************************/
787 :
788 57 : GDALDataset *NTv2Dataset::Create(const char *pszFilename, int nXSize,
789 : int nYSize, int nBandsIn, GDALDataType eType,
790 : char **papszOptions)
791 : {
792 57 : if (eType != GDT_Float32)
793 : {
794 44 : CPLError(CE_Failure, CPLE_AppDefined,
795 : "Attempt to create NTv2 file with unsupported data type '%s'.",
796 : GDALGetDataTypeName(eType));
797 44 : return nullptr;
798 : }
799 13 : if (nBandsIn != 4)
800 : {
801 5 : CPLError(CE_Failure, CPLE_AppDefined,
802 : "Attempt to create NTv2 file with unsupported "
803 : "band number '%d'.",
804 : nBandsIn);
805 5 : return nullptr;
806 : }
807 :
808 : /* -------------------------------------------------------------------- */
809 : /* Are we extending an existing file? */
810 : /* -------------------------------------------------------------------- */
811 8 : const bool bAppend = CPLFetchBool(papszOptions, "APPEND_SUBDATASET", false);
812 :
813 : /* -------------------------------------------------------------------- */
814 : /* Try to open or create file. */
815 : /* -------------------------------------------------------------------- */
816 8 : VSILFILE *fp = nullptr;
817 8 : if (bAppend)
818 3 : fp = VSIFOpenL(pszFilename, "rb+");
819 : else
820 5 : fp = VSIFOpenL(pszFilename, "wb");
821 :
822 8 : if (fp == nullptr)
823 : {
824 2 : CPLError(CE_Failure, CPLE_OpenFailed,
825 : "Attempt to open/create file `%s' failed.\n", pszFilename);
826 2 : return nullptr;
827 : }
828 :
829 : /* -------------------------------------------------------------------- */
830 : /* Create a file level header if we are creating new. */
831 : /* -------------------------------------------------------------------- */
832 6 : char achHeader[11 * 16] = {'\0'};
833 6 : const char *pszValue = nullptr;
834 6 : GUInt32 nNumFile = 1;
835 6 : bool bMustSwap = false;
836 6 : bool bIsLE = false;
837 :
838 6 : if (!bAppend)
839 : {
840 4 : memset(achHeader, 0, sizeof(achHeader));
841 :
842 4 : bIsLE =
843 4 : EQUAL(CSLFetchNameValueDef(papszOptions, "ENDIANNESS", "LE"), "LE");
844 : #ifdef CPL_LSB
845 4 : bMustSwap = !bIsLE;
846 : #else
847 : bMustSwap = bIsLE;
848 : #endif
849 :
850 4 : memcpy(achHeader + 0 * 16, "NUM_OREC", 8);
851 4 : int nNumOrec = 11;
852 4 : SwapPtr32IfNecessary(bMustSwap, &nNumOrec);
853 4 : memcpy(achHeader + 0 * 16 + 8, &nNumOrec, 4);
854 :
855 4 : memcpy(achHeader + 1 * 16, "NUM_SREC", 8);
856 4 : int nNumSrec = 11;
857 4 : SwapPtr32IfNecessary(bMustSwap, &nNumSrec);
858 4 : memcpy(achHeader + 1 * 16 + 8, &nNumSrec, 4);
859 :
860 4 : memcpy(achHeader + 2 * 16, "NUM_FILE", 8);
861 4 : SwapPtr32IfNecessary(bMustSwap, &nNumFile);
862 4 : memcpy(achHeader + 2 * 16 + 8, &nNumFile, 4);
863 4 : SwapPtr32IfNecessary(bMustSwap, &nNumFile);
864 :
865 4 : const size_t nMinLen = 16;
866 4 : memcpy(achHeader + 3 * 16, "GS_TYPE ", 16);
867 4 : pszValue = CSLFetchNameValueDef(papszOptions, "GS_TYPE", "SECONDS");
868 4 : memcpy(achHeader + 3 * 16 + 8, pszValue,
869 4 : std::min(nMinLen, strlen(pszValue)));
870 :
871 4 : memcpy(achHeader + 4 * 16, "VERSION ", 16);
872 4 : pszValue = CSLFetchNameValueDef(papszOptions, "VERSION", "");
873 4 : memcpy(achHeader + 4 * 16 + 8, pszValue,
874 4 : std::min(nMinLen, strlen(pszValue)));
875 :
876 4 : memcpy(achHeader + 5 * 16, "SYSTEM_F ", 16);
877 4 : pszValue = CSLFetchNameValueDef(papszOptions, "SYSTEM_F", "");
878 4 : memcpy(achHeader + 5 * 16 + 8, pszValue,
879 4 : std::min(nMinLen, strlen(pszValue)));
880 :
881 4 : memcpy(achHeader + 6 * 16, "SYSTEM_T ", 16);
882 4 : pszValue = CSLFetchNameValueDef(papszOptions, "SYSTEM_T", "");
883 4 : memcpy(achHeader + 6 * 16 + 8, pszValue,
884 4 : std::min(nMinLen, strlen(pszValue)));
885 :
886 4 : memcpy(achHeader + 7 * 16, "MAJOR_F ", 8);
887 4 : memcpy(achHeader + 8 * 16, "MINOR_F ", 8);
888 4 : memcpy(achHeader + 9 * 16, "MAJOR_T ", 8);
889 4 : memcpy(achHeader + 10 * 16, "MINOR_T ", 8);
890 :
891 4 : CPL_IGNORE_RET_VAL(VSIFWriteL(achHeader, 1, sizeof(achHeader), fp));
892 : }
893 :
894 : /* -------------------------------------------------------------------- */
895 : /* Otherwise update the header with an increased subfile count, */
896 : /* and advanced to the last record of the file. */
897 : /* -------------------------------------------------------------------- */
898 : else
899 : {
900 2 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, 0, SEEK_SET));
901 2 : CPL_IGNORE_RET_VAL(VSIFReadL(achHeader, 1, 16, fp));
902 :
903 3 : bIsLE = achHeader[8] == 11 && achHeader[9] == 0 && achHeader[10] == 0 &&
904 1 : achHeader[11] == 0;
905 1 : const bool bIsBE = achHeader[8] == 0 && achHeader[9] == 0 &&
906 3 : achHeader[10] == 0 && achHeader[11] == 11;
907 2 : if (!bIsLE && !bIsBE)
908 : {
909 0 : VSIFCloseL(fp);
910 0 : return nullptr;
911 : }
912 : #ifdef CPL_LSB
913 2 : bMustSwap = bIsBE;
914 : #else
915 : bMustSwap = bIsLE;
916 : #endif
917 :
918 2 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, 2 * 16 + 8, SEEK_SET));
919 2 : CPL_IGNORE_RET_VAL(VSIFReadL(&nNumFile, 1, 4, fp));
920 2 : SwapPtr32IfNecessary(bMustSwap, &nNumFile);
921 :
922 2 : nNumFile++;
923 :
924 2 : SwapPtr32IfNecessary(bMustSwap, &nNumFile);
925 2 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, 2 * 16 + 8, SEEK_SET));
926 2 : CPL_IGNORE_RET_VAL(VSIFWriteL(&nNumFile, 1, 4, fp));
927 2 : SwapPtr32IfNecessary(bMustSwap, &nNumFile);
928 :
929 2 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, 0, SEEK_END));
930 2 : const vsi_l_offset nEnd = VSIFTellL(fp);
931 2 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, nEnd - 16, SEEK_SET));
932 : }
933 :
934 : /* -------------------------------------------------------------------- */
935 : /* Write the grid header. */
936 : /* -------------------------------------------------------------------- */
937 6 : memset(achHeader, 0, sizeof(achHeader));
938 :
939 6 : const size_t nMinLen = 16;
940 :
941 6 : memcpy(achHeader + 0 * 16, "SUB_NAME ", 16);
942 6 : pszValue = CSLFetchNameValueDef(papszOptions, "SUB_NAME", "");
943 6 : memcpy(achHeader + 0 * 16 + 8, pszValue,
944 6 : std::min(nMinLen, strlen(pszValue)));
945 :
946 6 : memcpy(achHeader + 1 * 16, "PARENT ", 16);
947 6 : pszValue = CSLFetchNameValueDef(papszOptions, "PARENT", "NONE");
948 6 : memcpy(achHeader + 1 * 16 + 8, pszValue,
949 6 : std::min(nMinLen, strlen(pszValue)));
950 :
951 6 : memcpy(achHeader + 2 * 16, "CREATED ", 16);
952 6 : pszValue = CSLFetchNameValueDef(papszOptions, "CREATED", "");
953 6 : memcpy(achHeader + 2 * 16 + 8, pszValue,
954 6 : std::min(nMinLen, strlen(pszValue)));
955 :
956 6 : memcpy(achHeader + 3 * 16, "UPDATED ", 16);
957 6 : pszValue = CSLFetchNameValueDef(papszOptions, "UPDATED", "");
958 6 : memcpy(achHeader + 3 * 16 + 8, pszValue,
959 6 : std::min(nMinLen, strlen(pszValue)));
960 :
961 : double dfValue;
962 :
963 6 : memcpy(achHeader + 4 * 16, "S_LAT ", 8);
964 6 : dfValue = 0;
965 6 : SwapPtr64IfNecessary(bMustSwap, &dfValue);
966 6 : memcpy(achHeader + 4 * 16 + 8, &dfValue, 8);
967 :
968 6 : memcpy(achHeader + 5 * 16, "N_LAT ", 8);
969 6 : dfValue = nYSize - 1;
970 6 : SwapPtr64IfNecessary(bMustSwap, &dfValue);
971 6 : memcpy(achHeader + 5 * 16 + 8, &dfValue, 8);
972 :
973 6 : memcpy(achHeader + 6 * 16, "E_LONG ", 8);
974 6 : dfValue = -1 * (nXSize - 1);
975 6 : SwapPtr64IfNecessary(bMustSwap, &dfValue);
976 6 : memcpy(achHeader + 6 * 16 + 8, &dfValue, 8);
977 :
978 6 : memcpy(achHeader + 7 * 16, "W_LONG ", 8);
979 6 : dfValue = 0;
980 6 : SwapPtr64IfNecessary(bMustSwap, &dfValue);
981 6 : memcpy(achHeader + 7 * 16 + 8, &dfValue, 8);
982 :
983 6 : memcpy(achHeader + 8 * 16, "LAT_INC ", 8);
984 6 : dfValue = 1;
985 6 : SwapPtr64IfNecessary(bMustSwap, &dfValue);
986 6 : memcpy(achHeader + 8 * 16 + 8, &dfValue, 8);
987 :
988 6 : memcpy(achHeader + 9 * 16, "LONG_INC", 8);
989 6 : memcpy(achHeader + 9 * 16 + 8, &dfValue, 8);
990 :
991 6 : memcpy(achHeader + 10 * 16, "GS_COUNT", 8);
992 6 : GUInt32 nGSCount = nXSize * nYSize;
993 6 : SwapPtr32IfNecessary(bMustSwap, &nGSCount);
994 6 : memcpy(achHeader + 10 * 16 + 8, &nGSCount, 4);
995 :
996 6 : CPL_IGNORE_RET_VAL(VSIFWriteL(achHeader, 1, sizeof(achHeader), fp));
997 :
998 : /* -------------------------------------------------------------------- */
999 : /* Write zeroed grid data. */
1000 : /* -------------------------------------------------------------------- */
1001 6 : memset(achHeader, 0, 16);
1002 :
1003 : // Use -1 (0x000080bf) as the default error value.
1004 6 : memset(achHeader + ((bIsLE) ? 10 : 9), 0x80, 1);
1005 6 : memset(achHeader + ((bIsLE) ? 11 : 8), 0xbf, 1);
1006 6 : memset(achHeader + ((bIsLE) ? 14 : 13), 0x80, 1);
1007 6 : memset(achHeader + ((bIsLE) ? 15 : 12), 0xbf, 1);
1008 :
1009 24 : for (int i = 0; i < nXSize * nYSize; i++)
1010 18 : CPL_IGNORE_RET_VAL(VSIFWriteL(achHeader, 1, 16, fp));
1011 :
1012 : /* -------------------------------------------------------------------- */
1013 : /* Write the end record. */
1014 : /* -------------------------------------------------------------------- */
1015 6 : memcpy(achHeader, "END ", 8);
1016 6 : memset(achHeader + 8, 0, 8);
1017 6 : CPL_IGNORE_RET_VAL(VSIFWriteL(achHeader, 1, 16, fp));
1018 6 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
1019 :
1020 6 : if (nNumFile == 1)
1021 4 : return GDALDataset::FromHandle(GDALOpen(pszFilename, GA_Update));
1022 :
1023 4 : CPLString osSubDSName;
1024 2 : osSubDSName.Printf("NTv2:%d:%s", nNumFile - 1, pszFilename);
1025 2 : return GDALDataset::FromHandle(GDALOpen(osSubDSName, GA_Update));
1026 : }
1027 :
1028 : /************************************************************************/
1029 : /* GDALRegister_NTv2() */
1030 : /************************************************************************/
1031 :
1032 1595 : void GDALRegister_NTv2()
1033 :
1034 : {
1035 1595 : if (GDALGetDriverByName("NTv2") != nullptr)
1036 302 : return;
1037 :
1038 1293 : GDALDriver *poDriver = new GDALDriver();
1039 :
1040 1293 : poDriver->SetDescription("NTv2");
1041 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1042 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "NTv2 Datum Grid Shift");
1043 1293 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSIONS, "gsb gvb");
1044 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1045 1293 : poDriver->SetMetadataItem(GDAL_DMD_SUBDATASETS, "YES");
1046 1293 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Float32");
1047 :
1048 1293 : poDriver->pfnOpen = NTv2Dataset::Open;
1049 1293 : poDriver->pfnIdentify = NTv2Dataset::Identify;
1050 1293 : poDriver->pfnCreate = NTv2Dataset::Create;
1051 :
1052 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1053 : }
|
