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 GetGeoTransform(double *padfTransform) override;
106 :
107 2 : const OGRSpatialReference *GetSpatialRef() const override
108 : {
109 2 : return &m_oSRS;
110 : }
111 :
112 : static GDALDataset *Open(GDALOpenInfo *);
113 : static int Identify(GDALOpenInfo *);
114 : };
115 :
116 : /************************************************************************/
117 : /* ==================================================================== */
118 : /* NTv2Dataset */
119 : /* ==================================================================== */
120 : /************************************************************************/
121 :
122 : /************************************************************************/
123 : /* NTv2Dataset() */
124 : /************************************************************************/
125 :
126 4 : NTv2Dataset::NTv2Dataset()
127 : {
128 4 : m_oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
129 4 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
130 :
131 4 : adfGeoTransform[0] = 0.0;
132 4 : adfGeoTransform[1] = 0.0; // TODO(schwehr): Should this be 1.0?
133 4 : adfGeoTransform[2] = 0.0;
134 4 : adfGeoTransform[3] = 0.0;
135 4 : adfGeoTransform[4] = 0.0;
136 4 : adfGeoTransform[5] = 0.0; // TODO(schwehr): Should this be 1.0?
137 4 : }
138 :
139 : /************************************************************************/
140 : /* ~NTv2Dataset() */
141 : /************************************************************************/
142 :
143 8 : NTv2Dataset::~NTv2Dataset()
144 :
145 : {
146 4 : NTv2Dataset::Close();
147 8 : }
148 :
149 : /************************************************************************/
150 : /* Close() */
151 : /************************************************************************/
152 :
153 8 : CPLErr NTv2Dataset::Close()
154 : {
155 8 : CPLErr eErr = CE_None;
156 8 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
157 : {
158 4 : if (fpImage)
159 : {
160 4 : if (VSIFCloseL(fpImage) != 0)
161 : {
162 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
163 0 : eErr = CE_Failure;
164 : }
165 : }
166 :
167 4 : if (GDALPamDataset::Close() != CE_None)
168 0 : eErr = CE_Failure;
169 : }
170 8 : return eErr;
171 : }
172 :
173 : /************************************************************************/
174 : /* SwapPtr32IfNecessary() */
175 : /************************************************************************/
176 :
177 8 : static void SwapPtr32IfNecessary(bool bMustSwap, void *ptr)
178 : {
179 8 : if (bMustSwap)
180 : {
181 4 : CPL_SWAP32PTR(static_cast<GByte *>(ptr));
182 : }
183 8 : }
184 :
185 : /************************************************************************/
186 : /* SwapPtr64IfNecessary() */
187 : /************************************************************************/
188 :
189 40 : static void SwapPtr64IfNecessary(bool bMustSwap, void *ptr)
190 : {
191 40 : if (bMustSwap)
192 : {
193 20 : CPL_SWAP64PTR(static_cast<GByte *>(ptr));
194 : }
195 40 : }
196 :
197 : /************************************************************************/
198 : /* Identify() */
199 : /************************************************************************/
200 :
201 49609 : int NTv2Dataset::Identify(GDALOpenInfo *poOpenInfo)
202 :
203 : {
204 49609 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "NTv2:"))
205 0 : return TRUE;
206 :
207 49609 : if (poOpenInfo->nHeaderBytes < 64)
208 46254 : return FALSE;
209 :
210 3355 : const char *pszHeader =
211 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
212 3355 : if (!STARTS_WITH_CI(pszHeader + 0, "NUM_OREC"))
213 3347 : return FALSE;
214 :
215 8 : if (!STARTS_WITH_CI(pszHeader + knREGULAR_RECORD_SIZE, "NUM_SREC") &&
216 0 : !STARTS_WITH_CI(pszHeader + knMAX_RECORD_SIZE, "NUM_SREC"))
217 0 : return FALSE;
218 :
219 8 : return TRUE;
220 : }
221 :
222 : /************************************************************************/
223 : /* Open() */
224 : /************************************************************************/
225 :
226 4 : GDALDataset *NTv2Dataset::Open(GDALOpenInfo *poOpenInfo)
227 :
228 : {
229 4 : if (!Identify(poOpenInfo) || poOpenInfo->eAccess == GA_Update)
230 0 : return nullptr;
231 :
232 : /* -------------------------------------------------------------------- */
233 : /* Are we targeting a particular grid? */
234 : /* -------------------------------------------------------------------- */
235 8 : CPLString osFilename;
236 4 : int iTargetGrid = -1;
237 :
238 4 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "NTv2:"))
239 : {
240 0 : const char *pszRest = poOpenInfo->pszFilename + 5;
241 :
242 0 : iTargetGrid = atoi(pszRest);
243 0 : while (*pszRest != '\0' && *pszRest != ':')
244 0 : pszRest++;
245 :
246 0 : if (*pszRest == ':')
247 0 : pszRest++;
248 :
249 0 : osFilename = pszRest;
250 : }
251 : else
252 : {
253 4 : osFilename = poOpenInfo->pszFilename;
254 : }
255 :
256 : /* -------------------------------------------------------------------- */
257 : /* Create a corresponding GDALDataset. */
258 : /* -------------------------------------------------------------------- */
259 8 : auto poDS = std::make_unique<NTv2Dataset>();
260 4 : poDS->eAccess = poOpenInfo->eAccess;
261 :
262 : /* -------------------------------------------------------------------- */
263 : /* Open the file. */
264 : /* -------------------------------------------------------------------- */
265 4 : if (poOpenInfo->eAccess == GA_ReadOnly)
266 4 : poDS->fpImage = VSIFOpenL(osFilename, "rb");
267 : else
268 0 : poDS->fpImage = VSIFOpenL(osFilename, "rb+");
269 :
270 4 : if (poDS->fpImage == nullptr)
271 : {
272 0 : return nullptr;
273 : }
274 :
275 : /* -------------------------------------------------------------------- */
276 : /* Read the file header. */
277 : /* -------------------------------------------------------------------- */
278 4 : char achHeader[11 * knMAX_RECORD_SIZE] = {0};
279 :
280 8 : if (VSIFSeekL(poDS->fpImage, 0, SEEK_SET) != 0 ||
281 4 : VSIFReadL(achHeader, 1, 64, poDS->fpImage) != 64)
282 : {
283 0 : return nullptr;
284 : }
285 :
286 4 : poDS->nRecordSize =
287 4 : STARTS_WITH_CI(achHeader + knMAX_RECORD_SIZE, "NUM_SREC")
288 4 : ? knMAX_RECORD_SIZE
289 : : knREGULAR_RECORD_SIZE;
290 4 : if (VSIFReadL(achHeader + 64, 1, 11 * poDS->nRecordSize - 64,
291 8 : poDS->fpImage) != 11 * poDS->nRecordSize - 64)
292 : {
293 0 : return nullptr;
294 : }
295 :
296 2 : const bool bIsLE = achHeader[8] == 11 && achHeader[9] == 0 &&
297 6 : achHeader[10] == 0 && achHeader[11] == 0;
298 2 : const bool bIsBE = achHeader[8] == 0 && achHeader[9] == 0 &&
299 6 : achHeader[10] == 0 && achHeader[11] == 11;
300 4 : if (!bIsLE && !bIsBE)
301 : {
302 0 : return nullptr;
303 : }
304 4 : poDS->m_eByteOrder = bIsLE ? RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN
305 : : RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN;
306 4 : poDS->m_bMustSwap = poDS->m_eByteOrder != RawRasterBand::NATIVE_BYTE_ORDER;
307 :
308 4 : SwapPtr32IfNecessary(poDS->m_bMustSwap,
309 4 : achHeader + 2 * poDS->nRecordSize + 8);
310 4 : GInt32 nSubFileCount = 0;
311 4 : memcpy(&nSubFileCount, achHeader + 2 * poDS->nRecordSize + 8, 4);
312 4 : if (nSubFileCount <= 0 || nSubFileCount >= 1024)
313 : {
314 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for NUM_FILE : %d",
315 : nSubFileCount);
316 0 : return nullptr;
317 : }
318 :
319 4 : poDS->CaptureMetadataItem(achHeader + 3 * poDS->nRecordSize);
320 4 : poDS->CaptureMetadataItem(achHeader + 4 * poDS->nRecordSize);
321 4 : poDS->CaptureMetadataItem(achHeader + 5 * poDS->nRecordSize);
322 4 : poDS->CaptureMetadataItem(achHeader + 6 * poDS->nRecordSize);
323 :
324 4 : double dfValue = 0.0;
325 4 : memcpy(&dfValue, achHeader + 7 * poDS->nRecordSize + 8, 8);
326 4 : SwapPtr64IfNecessary(poDS->m_bMustSwap, &dfValue);
327 8 : CPLString osFValue;
328 4 : osFValue.Printf("%.15g", dfValue);
329 4 : poDS->SetMetadataItem("MAJOR_F", osFValue);
330 :
331 4 : memcpy(&dfValue, achHeader + 8 * poDS->nRecordSize + 8, 8);
332 4 : SwapPtr64IfNecessary(poDS->m_bMustSwap, &dfValue);
333 4 : osFValue.Printf("%.15g", dfValue);
334 4 : poDS->SetMetadataItem("MINOR_F", osFValue);
335 :
336 4 : memcpy(&dfValue, achHeader + 9 * poDS->nRecordSize + 8, 8);
337 4 : SwapPtr64IfNecessary(poDS->m_bMustSwap, &dfValue);
338 4 : osFValue.Printf("%.15g", dfValue);
339 4 : poDS->SetMetadataItem("MAJOR_T", osFValue);
340 :
341 4 : memcpy(&dfValue, achHeader + 10 * poDS->nRecordSize + 8, 8);
342 4 : SwapPtr64IfNecessary(poDS->m_bMustSwap, &dfValue);
343 4 : osFValue.Printf("%.15g", dfValue);
344 4 : poDS->SetMetadataItem("MINOR_T", osFValue);
345 :
346 : /* ==================================================================== */
347 : /* Loop over grids. */
348 : /* ==================================================================== */
349 4 : vsi_l_offset nGridOffset = 11 * poDS->nRecordSize;
350 :
351 8 : for (int iGrid = 0; iGrid < nSubFileCount; iGrid++)
352 : {
353 8 : if (VSIFSeekL(poDS->fpImage, nGridOffset, SEEK_SET) < 0 ||
354 4 : VSIFReadL(achHeader, 11, poDS->nRecordSize, poDS->fpImage) !=
355 4 : poDS->nRecordSize)
356 : {
357 0 : CPLError(CE_Failure, CPLE_AppDefined,
358 : "Cannot read header for subfile %d", iGrid);
359 0 : return nullptr;
360 : }
361 :
362 28 : for (int i = 4; i <= 9; i++)
363 24 : SwapPtr64IfNecessary(poDS->m_bMustSwap,
364 24 : achHeader + i * poDS->nRecordSize + 8);
365 :
366 4 : SwapPtr32IfNecessary(poDS->m_bMustSwap,
367 4 : achHeader + 10 * poDS->nRecordSize + 8);
368 :
369 4 : GUInt32 nGSCount = 0;
370 4 : memcpy(&nGSCount, achHeader + 10 * poDS->nRecordSize + 8, 4);
371 :
372 4 : CPLString osSubName;
373 4 : osSubName.assign(achHeader + 8, 8);
374 4 : osSubName.Trim();
375 :
376 : // If this is our target grid, open it as a dataset.
377 4 : if (iTargetGrid == iGrid || (iTargetGrid == -1 && iGrid == 0))
378 : {
379 4 : if (!poDS->OpenGrid(achHeader, nGridOffset))
380 : {
381 0 : return nullptr;
382 : }
383 : }
384 :
385 : // If we are opening the file as a whole, list subdatasets.
386 4 : if (iTargetGrid == -1)
387 : {
388 8 : CPLString osKey;
389 8 : CPLString osValue;
390 4 : osKey.Printf("SUBDATASET_%d_NAME", iGrid);
391 4 : osValue.Printf("NTv2:%d:%s", iGrid, osFilename.c_str());
392 4 : poDS->SetMetadataItem(osKey, osValue, "SUBDATASETS");
393 :
394 4 : osKey.Printf("SUBDATASET_%d_DESC", iGrid);
395 4 : osValue.Printf("%s", osSubName.c_str());
396 4 : poDS->SetMetadataItem(osKey, osValue, "SUBDATASETS");
397 : }
398 :
399 4 : nGridOffset +=
400 4 : (11 + static_cast<vsi_l_offset>(nGSCount)) * poDS->nRecordSize;
401 : }
402 :
403 : /* -------------------------------------------------------------------- */
404 : /* Initialize any PAM information. */
405 : /* -------------------------------------------------------------------- */
406 4 : poDS->SetDescription(poOpenInfo->pszFilename);
407 4 : poDS->TryLoadXML();
408 :
409 : /* -------------------------------------------------------------------- */
410 : /* Check for overviews. */
411 : /* -------------------------------------------------------------------- */
412 4 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
413 :
414 4 : return poDS.release();
415 : }
416 :
417 : /************************************************************************/
418 : /* OpenGrid() */
419 : /* */
420 : /* Note that the caller will already have byte swapped needed */
421 : /* portions of the header. */
422 : /************************************************************************/
423 :
424 4 : bool NTv2Dataset::OpenGrid(const char *pachHeader, vsi_l_offset nGridOffsetIn)
425 :
426 : {
427 4 : nGridOffset = nGridOffsetIn;
428 :
429 : /* -------------------------------------------------------------------- */
430 : /* Read the grid header. */
431 : /* -------------------------------------------------------------------- */
432 4 : CaptureMetadataItem(pachHeader + 0 * nRecordSize);
433 4 : CaptureMetadataItem(pachHeader + 1 * nRecordSize);
434 4 : CaptureMetadataItem(pachHeader + 2 * nRecordSize);
435 4 : CaptureMetadataItem(pachHeader + 3 * nRecordSize);
436 :
437 : double s_lat, n_lat, e_long, w_long, lat_inc, long_inc;
438 4 : memcpy(&s_lat, pachHeader + 4 * nRecordSize + 8, 8);
439 4 : memcpy(&n_lat, pachHeader + 5 * nRecordSize + 8, 8);
440 4 : memcpy(&e_long, pachHeader + 6 * nRecordSize + 8, 8);
441 4 : memcpy(&w_long, pachHeader + 7 * nRecordSize + 8, 8);
442 4 : memcpy(&lat_inc, pachHeader + 8 * nRecordSize + 8, 8);
443 4 : memcpy(&long_inc, pachHeader + 9 * nRecordSize + 8, 8);
444 :
445 4 : e_long *= -1;
446 4 : w_long *= -1;
447 :
448 4 : if (long_inc == 0.0 || lat_inc == 0.0)
449 0 : return false;
450 4 : const double dfXSize = floor((e_long - w_long) / long_inc + 1.5);
451 4 : const double dfYSize = floor((n_lat - s_lat) / lat_inc + 1.5);
452 4 : if (!(dfXSize >= 0 && dfXSize < INT_MAX) ||
453 4 : !(dfYSize >= 0 && dfYSize < INT_MAX))
454 0 : return false;
455 4 : nRasterXSize = static_cast<int>(dfXSize);
456 4 : nRasterYSize = static_cast<int>(dfYSize);
457 :
458 4 : const int l_nBands = nRecordSize == knREGULAR_RECORD_SIZE ? 4 : 6;
459 4 : const int nPixelSize = l_nBands * 4;
460 :
461 4 : if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
462 0 : return false;
463 4 : if (nRasterXSize > INT_MAX / nPixelSize)
464 0 : return false;
465 :
466 : /* -------------------------------------------------------------------- */
467 : /* Create band information object. */
468 : /* */
469 : /* We use unusual offsets to remap from bottom to top, to top */
470 : /* to bottom orientation, and also to remap east to west, to */
471 : /* west to east. */
472 : /* -------------------------------------------------------------------- */
473 20 : for (int iBand = 0; iBand < l_nBands; iBand++)
474 : {
475 : auto poBand = RawRasterBand::Create(
476 : this, iBand + 1, fpImage,
477 16 : nGridOffset + 4 * iBand + 11 * nRecordSize +
478 16 : static_cast<vsi_l_offset>(nRasterXSize - 1) * nPixelSize +
479 16 : static_cast<vsi_l_offset>(nRasterYSize - 1) * nPixelSize *
480 16 : nRasterXSize,
481 16 : -nPixelSize, -nPixelSize * nRasterXSize, GDT_Float32, m_eByteOrder,
482 16 : RawRasterBand::OwnFP::NO);
483 16 : if (!poBand)
484 0 : return false;
485 16 : SetBand(iBand + 1, std::move(poBand));
486 : }
487 :
488 4 : if (l_nBands == 4)
489 : {
490 4 : GetRasterBand(1)->SetDescription("Latitude Offset (arc seconds)");
491 4 : GetRasterBand(2)->SetDescription("Longitude Offset (arc seconds)");
492 4 : GetRasterBand(2)->SetMetadataItem("positive_value", "west");
493 4 : GetRasterBand(3)->SetDescription("Latitude Error");
494 4 : GetRasterBand(4)->SetDescription("Longitude Error");
495 : }
496 : else
497 : {
498 : // A bit surprising that the order is easting, northing here, contrary
499 : // to the classic NTv2 order.... Verified on NAD83v70VG.gvb
500 : // (https://webapp.geod.nrcan.gc.ca/geod/process/download-helper.php?file_id=NAD83v70VG)
501 : // against the TRX software
502 : // (https://webapp.geod.nrcan.gc.ca/geod/process/download-helper.php?file_id=trx)
503 : // https://webapp.geod.nrcan.gc.ca/geod/tools-outils/nad83-docs.php
504 : // Unfortunately I couldn't find an official documentation of the format
505 : // !
506 0 : GetRasterBand(1)->SetDescription("East velocity (mm/year)");
507 0 : GetRasterBand(2)->SetDescription("North velocity (mm/year)");
508 0 : GetRasterBand(3)->SetDescription("Up velocity (mm/year)");
509 0 : GetRasterBand(4)->SetDescription("East velocity Error (mm/year)");
510 0 : GetRasterBand(5)->SetDescription("North velocity Error (mm/year)");
511 0 : GetRasterBand(6)->SetDescription("Up velocity Error (mm/year)");
512 : }
513 :
514 : /* -------------------------------------------------------------------- */
515 : /* Setup georeferencing. */
516 : /* -------------------------------------------------------------------- */
517 4 : adfGeoTransform[0] = (w_long - long_inc * 0.5) / 3600.0;
518 4 : adfGeoTransform[1] = long_inc / 3600.0;
519 4 : adfGeoTransform[2] = 0.0;
520 4 : adfGeoTransform[3] = (n_lat + lat_inc * 0.5) / 3600.0;
521 4 : adfGeoTransform[4] = 0.0;
522 4 : adfGeoTransform[5] = (-1 * lat_inc) / 3600.0;
523 :
524 4 : return true;
525 : }
526 :
527 : /************************************************************************/
528 : /* CaptureMetadataItem() */
529 : /************************************************************************/
530 :
531 32 : void NTv2Dataset::CaptureMetadataItem(const char *pszItem)
532 :
533 : {
534 64 : CPLString osKey;
535 64 : CPLString osValue;
536 :
537 32 : osKey.assign(pszItem, 8);
538 32 : osValue.assign(pszItem + 8, 8);
539 :
540 32 : SetMetadataItem(osKey.Trim(), osValue.Trim());
541 32 : }
542 :
543 : /************************************************************************/
544 : /* GetGeoTransform() */
545 : /************************************************************************/
546 :
547 2 : CPLErr NTv2Dataset::GetGeoTransform(double *padfTransform)
548 :
549 : {
550 2 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
551 2 : return CE_None;
552 : }
553 :
554 : /************************************************************************/
555 : /* GDALRegister_NTv2() */
556 : /************************************************************************/
557 :
558 1667 : void GDALRegister_NTv2()
559 :
560 : {
561 1667 : if (GDALGetDriverByName("NTv2") != nullptr)
562 282 : return;
563 :
564 1385 : GDALDriver *poDriver = new GDALDriver();
565 :
566 1385 : poDriver->SetDescription("NTv2");
567 1385 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
568 1385 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "NTv2 Datum Grid Shift");
569 1385 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSIONS, "gsb gvb");
570 1385 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
571 1385 : poDriver->SetMetadataItem(GDAL_DMD_SUBDATASETS, "YES");
572 :
573 1385 : poDriver->pfnOpen = NTv2Dataset::Open;
574 1385 : poDriver->pfnIdentify = NTv2Dataset::Identify;
575 :
576 1385 : GetGDALDriverManager()->RegisterDriver(poDriver);
577 : }
|
