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