LCOV - code coverage report
Current view: top level - frmts/raw - ntv2dataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 188 223 84.3 %
Date: 2025-07-05 13:22:42 Functions: 13 13 100.0 %

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

Generated by: LCOV version 1.14