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