Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Horizontal Datum Formats
       4             :  * Purpose:  Implementation of the CTable2 format, a PROJ.4 specific format
       5             :  *           that is more compact (due to a lack of unused error band) than NTv2
       6             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       7             :  *
       8             :  ******************************************************************************
       9             :  * Copyright (c) 2012, Frank Warmerdam
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "cpl_string.h"
      15             : #include "gdal_frmts.h"
      16             : #include "ogr_srs_api.h"
      17             : #include "rawdataset.h"
      18             : 
      19             : #include <algorithm>
      20             : 
      21             : /************************************************************************/
      22             : /* ==================================================================== */
      23             : /*                              CTable2Dataset                          */
      24             : /* ==================================================================== */
      25             : /************************************************************************/
      26             : 
      27             : class CTable2Dataset final : public RawDataset
      28             : {
      29             :     VSILFILE *fpImage;  // image data file.
      30             : 
      31             :     double adfGeoTransform[6];
      32             :     OGRSpatialReference m_oSRS{};
      33             : 
      34             :     CPL_DISALLOW_COPY_ASSIGN(CTable2Dataset)
      35             : 
      36             :     CPLErr Close() override;
      37             : 
      38             :   public:
      39             :     CTable2Dataset();
      40             :     ~CTable2Dataset() override;
      41             : 
      42             :     CPLErr SetGeoTransform(double *padfTransform) override;
      43             :     CPLErr GetGeoTransform(double *padfTransform) override;
      44             : 
      45           0 :     const OGRSpatialReference *GetSpatialRef() const override
      46             :     {
      47           0 :         return &m_oSRS;
      48             :     }
      49             : 
      50             :     static GDALDataset *Open(GDALOpenInfo *);
      51             :     static int Identify(GDALOpenInfo *);
      52             :     static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
      53             :                                int nBands, GDALDataType eType,
      54             :                                char **papszOptions);
      55             : };
      56             : 
      57             : /************************************************************************/
      58             : /* ==================================================================== */
      59             : /*                              CTable2Dataset                          */
      60             : /* ==================================================================== */
      61             : /************************************************************************/
      62             : 
      63             : /************************************************************************/
      64             : /*                             CTable2Dataset()                          */
      65             : /************************************************************************/
      66             : 
      67           6 : CTable2Dataset::CTable2Dataset() : fpImage(nullptr)
      68             : {
      69           6 :     m_oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
      70           6 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
      71             : 
      72           6 :     memset(adfGeoTransform, 0, sizeof(adfGeoTransform));
      73           6 : }
      74             : 
      75             : /************************************************************************/
      76             : /*                            ~CTable2Dataset()                         */
      77             : /************************************************************************/
      78             : 
      79          12 : CTable2Dataset::~CTable2Dataset()
      80             : 
      81             : {
      82           6 :     CTable2Dataset::Close();
      83          12 : }
      84             : 
      85             : /************************************************************************/
      86             : /*                              Close()                                 */
      87             : /************************************************************************/
      88             : 
      89          11 : CPLErr CTable2Dataset::Close()
      90             : {
      91          11 :     CPLErr eErr = CE_None;
      92          11 :     if (nOpenFlags != OPEN_FLAGS_CLOSED)
      93             :     {
      94           6 :         if (CTable2Dataset::FlushCache(true) != CE_None)
      95           0 :             eErr = CE_Failure;
      96             : 
      97           6 :         if (fpImage != nullptr)
      98             :         {
      99           6 :             if (VSIFCloseL(fpImage) != 0)
     100             :             {
     101           0 :                 CPLError(CE_Failure, CPLE_FileIO, "I/O error");
     102           0 :                 eErr = CE_Failure;
     103             :             }
     104             :         }
     105             : 
     106           6 :         if (GDALPamDataset::Close() != CE_None)
     107           0 :             eErr = CE_Failure;
     108             :     }
     109          11 :     return eErr;
     110             : }
     111             : 
     112             : /************************************************************************/
     113             : /*                              Identify()                              */
     114             : /************************************************************************/
     115             : 
     116       51445 : int CTable2Dataset::Identify(GDALOpenInfo *poOpenInfo)
     117             : 
     118             : {
     119       51445 :     if (poOpenInfo->nHeaderBytes < 64)
     120       48068 :         return FALSE;
     121             : 
     122        3377 :     if (!STARTS_WITH_CI(
     123             :             reinterpret_cast<const char *>(poOpenInfo->pabyHeader + 0),
     124             :             "CTABLE V2"))
     125        3365 :         return FALSE;
     126             : 
     127          12 :     return TRUE;
     128             : }
     129             : 
     130             : /************************************************************************/
     131             : /*                                Open()                                */
     132             : /************************************************************************/
     133             : 
     134           6 : GDALDataset *CTable2Dataset::Open(GDALOpenInfo *poOpenInfo)
     135             : 
     136             : {
     137           6 :     if (!Identify(poOpenInfo) || !poOpenInfo->fpL)
     138           0 :         return nullptr;
     139             : 
     140             :     /* -------------------------------------------------------------------- */
     141             :     /*      Create a corresponding GDALDataset.                             */
     142             :     /* -------------------------------------------------------------------- */
     143          12 :     auto poDS = std::make_unique<CTable2Dataset>();
     144           6 :     poDS->eAccess = poOpenInfo->eAccess;
     145           6 :     std::swap(poDS->fpImage, poOpenInfo->fpL);
     146             : 
     147             :     /* -------------------------------------------------------------------- */
     148             :     /*      Read the file header.                                           */
     149             :     /* -------------------------------------------------------------------- */
     150             : 
     151           6 :     CPL_IGNORE_RET_VAL(VSIFSeekL(poDS->fpImage, 0, SEEK_SET));
     152             : 
     153           6 :     char achHeader[160] = {'\0'};
     154           6 :     CPL_IGNORE_RET_VAL(VSIFReadL(achHeader, 1, 160, poDS->fpImage));
     155           6 :     achHeader[16 + 79] = '\0';
     156             : 
     157          12 :     CPLString osDescription = reinterpret_cast<const char *>(achHeader + 16);
     158           6 :     osDescription.Trim();
     159           6 :     poDS->SetMetadataItem("DESCRIPTION", osDescription);
     160             : 
     161             :     /* -------------------------------------------------------------------- */
     162             :     /*      Convert from LSB to local machine byte order.                   */
     163             :     /* -------------------------------------------------------------------- */
     164           6 :     CPL_LSBPTR64(achHeader + 96);
     165           6 :     CPL_LSBPTR64(achHeader + 104);
     166           6 :     CPL_LSBPTR64(achHeader + 112);
     167           6 :     CPL_LSBPTR64(achHeader + 120);
     168           6 :     CPL_LSBPTR32(achHeader + 128);
     169           6 :     CPL_LSBPTR32(achHeader + 132);
     170             : 
     171             :     /* -------------------------------------------------------------------- */
     172             :     /*      Extract size, and geotransform.                                 */
     173             :     /* -------------------------------------------------------------------- */
     174             :     int nRasterXSize, nRasterYSize;
     175           6 :     memcpy(&nRasterXSize, achHeader + 128, 4);
     176           6 :     memcpy(&nRasterYSize, achHeader + 132, 4);
     177          12 :     if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize) ||
     178             :         /* to avoid overflow in later -8 * nRasterXSize computation */
     179           6 :         nRasterXSize >= INT_MAX / 8)
     180             :     {
     181           0 :         return nullptr;
     182             :     }
     183             : 
     184           6 :     poDS->nRasterXSize = nRasterXSize;
     185           6 :     poDS->nRasterYSize = nRasterYSize;
     186             : 
     187             :     double adfValues[4];
     188           6 :     memcpy(adfValues, achHeader + 96, sizeof(double) * 4);
     189             : 
     190          30 :     for (int i = 0; i < 4; i++)
     191          24 :         adfValues[i] *= 180 / M_PI;  // Radians to degrees.
     192             : 
     193           6 :     poDS->adfGeoTransform[0] = adfValues[0] - adfValues[2] * 0.5;
     194           6 :     poDS->adfGeoTransform[1] = adfValues[2];
     195           6 :     poDS->adfGeoTransform[2] = 0.0;
     196          12 :     poDS->adfGeoTransform[3] =
     197           6 :         adfValues[1] + adfValues[3] * (nRasterYSize - 0.5);
     198           6 :     poDS->adfGeoTransform[4] = 0.0;
     199           6 :     poDS->adfGeoTransform[5] = -adfValues[3];
     200             : 
     201             :     /* -------------------------------------------------------------------- */
     202             :     /*      Setup the bands.                                                */
     203             :     /* -------------------------------------------------------------------- */
     204             :     auto poBand =
     205          12 :         RawRasterBand::Create(poDS.get(), 1, poDS->fpImage,
     206             :                               160 + 4 +
     207           6 :                                   static_cast<vsi_l_offset>(nRasterXSize) *
     208           6 :                                       (nRasterYSize - 1) * 2 * 4,
     209             :                               8, -8 * nRasterXSize, GDT_Float32,
     210             :                               RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
     211          12 :                               RawRasterBand::OwnFP::NO);
     212           6 :     if (!poBand)
     213           0 :         return nullptr;
     214           6 :     poBand->SetDescription("Latitude Offset (radians)");
     215           6 :     poDS->SetBand(1, std::move(poBand));
     216             : 
     217             :     poBand =
     218          12 :         RawRasterBand::Create(poDS.get(), 2, poDS->fpImage,
     219           6 :                               160 + static_cast<vsi_l_offset>(nRasterXSize) *
     220           6 :                                         (nRasterYSize - 1) * 2 * 4,
     221             :                               8, -8 * nRasterXSize, GDT_Float32,
     222             :                               RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
     223           6 :                               RawRasterBand::OwnFP::NO);
     224           6 :     if (!poBand)
     225           0 :         return nullptr;
     226           6 :     poBand->SetDescription("Longitude Offset (radians)");
     227           6 :     poBand->SetMetadataItem("positive_value", "west");
     228           6 :     poDS->SetBand(2, std::move(poBand));
     229             : 
     230             :     /* -------------------------------------------------------------------- */
     231             :     /*      Initialize any PAM information.                                 */
     232             :     /* -------------------------------------------------------------------- */
     233           6 :     poDS->SetDescription(poOpenInfo->pszFilename);
     234           6 :     poDS->TryLoadXML();
     235             : 
     236             :     /* -------------------------------------------------------------------- */
     237             :     /*      Check for overviews.                                            */
     238             :     /* -------------------------------------------------------------------- */
     239           6 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     240             : 
     241           6 :     return poDS.release();
     242             : }
     243             : 
     244             : /************************************************************************/
     245             : /*                          GetGeoTransform()                           */
     246             : /************************************************************************/
     247             : 
     248           2 : CPLErr CTable2Dataset::GetGeoTransform(double *padfTransform)
     249             : 
     250             : {
     251           2 :     memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
     252           2 :     return CE_None;
     253             : }
     254             : 
     255             : /************************************************************************/
     256             : /*                          SetGeoTransform()                           */
     257             : /************************************************************************/
     258             : 
     259           2 : CPLErr CTable2Dataset::SetGeoTransform(double *padfTransform)
     260             : 
     261             : {
     262           2 :     if (eAccess == GA_ReadOnly)
     263             :     {
     264           0 :         CPLError(CE_Failure, CPLE_NoWriteAccess,
     265             :                  "Unable to update geotransform on readonly file.");
     266           0 :         return CE_Failure;
     267             :     }
     268             : 
     269           2 :     if (padfTransform[2] != 0.0 || padfTransform[4] != 0.0)
     270             :     {
     271           0 :         CPLError(
     272             :             CE_Failure, CPLE_AppDefined,
     273             :             "Rotated and sheared geotransforms not supported for CTable2.");
     274           0 :         return CE_Failure;
     275             :     }
     276             : 
     277           2 :     memcpy(adfGeoTransform, padfTransform, sizeof(double) * 6);
     278             : 
     279             :     /* -------------------------------------------------------------------- */
     280             :     /*      Update grid header.                                             */
     281             :     /* -------------------------------------------------------------------- */
     282           2 :     const double dfDegToRad = M_PI / 180.0;
     283             : 
     284             :     // read grid header
     285           2 :     CPL_IGNORE_RET_VAL(VSIFSeekL(fpImage, 0, SEEK_SET));
     286             : 
     287           2 :     char achHeader[160] = {'\0'};
     288           2 :     CPL_IGNORE_RET_VAL(VSIFReadL(achHeader, 1, sizeof(achHeader), fpImage));
     289             : 
     290             :     // lower left origin (longitude, center of pixel, radians)
     291           2 :     double dfValue =
     292           2 :         (adfGeoTransform[0] + adfGeoTransform[1] * 0.5) * dfDegToRad;
     293           2 :     CPL_LSBPTR64(&dfValue);
     294           2 :     memcpy(achHeader + 96, &dfValue, 8);
     295             : 
     296             :     // lower left origin (latitude, center of pixel, radians)
     297           2 :     dfValue = (adfGeoTransform[3] + adfGeoTransform[5] * (nRasterYSize - 0.5)) *
     298             :               dfDegToRad;
     299           2 :     CPL_LSBPTR64(&dfValue);
     300           2 :     memcpy(achHeader + 104, &dfValue, 8);
     301             : 
     302             :     // pixel width (radians)
     303           2 :     dfValue = adfGeoTransform[1] * dfDegToRad;
     304           2 :     CPL_LSBPTR64(&dfValue);
     305           2 :     memcpy(achHeader + 112, &dfValue, 8);
     306             : 
     307             :     // pixel height (radians)
     308           2 :     dfValue = adfGeoTransform[5] * -1 * dfDegToRad;
     309           2 :     CPL_LSBPTR64(&dfValue);
     310           2 :     memcpy(achHeader + 120, &dfValue, 8);
     311             : 
     312             :     // write grid header.
     313           2 :     CPL_IGNORE_RET_VAL(VSIFSeekL(fpImage, 0, SEEK_SET));
     314           2 :     CPL_IGNORE_RET_VAL(VSIFWriteL(achHeader, 1, sizeof(achHeader), fpImage));
     315             : 
     316           2 :     return CE_None;
     317             : }
     318             : 
     319             : /************************************************************************/
     320             : /*                               Create()                               */
     321             : /************************************************************************/
     322             : 
     323          49 : GDALDataset *CTable2Dataset::Create(const char *pszFilename, int nXSize,
     324             :                                     int nYSize, int /* nBandsIn */,
     325             :                                     GDALDataType eType, char **papszOptions)
     326             : {
     327          49 :     if (eType != GDT_Float32)
     328             :     {
     329          46 :         CPLError(
     330             :             CE_Failure, CPLE_AppDefined,
     331             :             "Attempt to create CTable2 file with unsupported data type '%s'.",
     332             :             GDALGetDataTypeName(eType));
     333          46 :         return nullptr;
     334             :     }
     335             : 
     336             :     /* -------------------------------------------------------------------- */
     337             :     /*      Try to open or create file.                                     */
     338             :     /* -------------------------------------------------------------------- */
     339           3 :     VSILFILE *fp = VSIFOpenL(pszFilename, "wb");
     340             : 
     341           3 :     if (fp == nullptr)
     342             :     {
     343           0 :         CPLError(CE_Failure, CPLE_OpenFailed,
     344             :                  "Attempt to create file `%s' failed.\n", pszFilename);
     345           0 :         return nullptr;
     346             :     }
     347             : 
     348             :     /* -------------------------------------------------------------------- */
     349             :     /*      Create a file header, with a defaulted georeferencing.          */
     350             :     /* -------------------------------------------------------------------- */
     351           3 :     char achHeader[160] = {'\0'};
     352             : 
     353           3 :     memset(achHeader, 0, sizeof(achHeader));
     354             : 
     355           3 :     memcpy(achHeader + 0, "CTABLE V2.0     ", 16);
     356             : 
     357           3 :     if (CSLFetchNameValue(papszOptions, "DESCRIPTION") != nullptr)
     358           0 :         strncpy(achHeader + 16, CSLFetchNameValue(papszOptions, "DESCRIPTION"),
     359             :                 80);
     360             : 
     361             :     // lower left origin (longitude, center of pixel, radians)
     362           3 :     double dfValue = 0;
     363           3 :     CPL_LSBPTR64(&dfValue);
     364           3 :     memcpy(achHeader + 96, &dfValue, 8);
     365             : 
     366             :     // lower left origin (latitude, center of pixel, radians)
     367           3 :     dfValue = 0;
     368           3 :     CPL_LSBPTR64(&dfValue);
     369           3 :     memcpy(achHeader + 104, &dfValue, 8);
     370             : 
     371             :     // pixel width (radians)
     372           3 :     dfValue = 0.01 * M_PI / 180.0;
     373           3 :     CPL_LSBPTR64(&dfValue);
     374           3 :     memcpy(achHeader + 112, &dfValue, 8);
     375             : 
     376             :     // pixel height (radians)
     377           3 :     dfValue = 0.01 * M_PI / 180.0;
     378           3 :     CPL_LSBPTR64(&dfValue);
     379           3 :     memcpy(achHeader + 120, &dfValue, 8);
     380             : 
     381             :     // raster width in pixels
     382           3 :     int nValue32 = nXSize;
     383           3 :     CPL_LSBPTR32(&nValue32);
     384           3 :     memcpy(achHeader + 128, &nValue32, 4);
     385             : 
     386             :     // raster width in pixels
     387           3 :     nValue32 = nYSize;
     388           3 :     CPL_LSBPTR32(&nValue32);
     389           3 :     memcpy(achHeader + 132, &nValue32, 4);
     390             : 
     391           3 :     CPL_IGNORE_RET_VAL(VSIFWriteL(achHeader, 1, sizeof(achHeader), fp));
     392             : 
     393             :     /* -------------------------------------------------------------------- */
     394             :     /*      Write zeroed grid data.                                         */
     395             :     /* -------------------------------------------------------------------- */
     396           3 :     float *pafLine = static_cast<float *>(CPLCalloc(sizeof(float) * 2, nXSize));
     397             : 
     398         213 :     for (int i = 0; i < nYSize; i++)
     399             :     {
     400         210 :         if (static_cast<int>(
     401         210 :                 VSIFWriteL(pafLine, sizeof(float) * 2, nXSize, fp)) != nXSize)
     402             :         {
     403           0 :             CPLError(CE_Failure, CPLE_FileIO,
     404             :                      "Write failed at line %d, perhaps the disk is full?", i);
     405           0 :             return nullptr;
     406             :         }
     407             :     }
     408             : 
     409             :     /* -------------------------------------------------------------------- */
     410             :     /*      Cleanup and return.                                             */
     411             :     /* -------------------------------------------------------------------- */
     412           3 :     CPLFree(pafLine);
     413             : 
     414           3 :     if (VSIFCloseL(fp) != 0)
     415             :     {
     416           0 :         CPLError(CE_Failure, CPLE_FileIO, "I/O error");
     417           0 :         return nullptr;
     418             :     }
     419             : 
     420           3 :     return static_cast<GDALDataset *>(GDALOpen(pszFilename, GA_Update));
     421             : }
     422             : 
     423             : /************************************************************************/
     424             : /*                         GDALRegister_CTable2()                       */
     425             : /************************************************************************/
     426             : 
     427        1595 : void GDALRegister_CTable2()
     428             : 
     429             : {
     430        1595 :     if (GDALGetDriverByName("CTable2") != nullptr)
     431         302 :         return;
     432             : 
     433        1293 :     GDALDriver *poDriver = new GDALDriver();
     434             : 
     435        1293 :     poDriver->SetDescription("CTable2");
     436        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     437        1293 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "CTable2 Datum Grid Shift");
     438        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     439             : 
     440        1293 :     poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Float32");
     441             : 
     442        1293 :     poDriver->pfnOpen = CTable2Dataset::Open;
     443        1293 :     poDriver->pfnIdentify = CTable2Dataset::Identify;
     444        1293 :     poDriver->pfnCreate = CTable2Dataset::Create;
     445             : 
     446        1293 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     447             : }