ctable2dataset.cpp

/******************************************************************************
 *
 * Project:  Horizontal Datum Formats
 * Purpose:  Implementation of the CTable2 format, a PROJ.4 specific format
 *           that is more compact (due to a lack of unused error band) than NTv2
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 2012, Frank Warmerdam
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ****************************************************************************/

#include "cpl_string.h"
#include "gdal_frmts.h"
#include "ogr_srs_api.h"
#include "rawdataset.h"

#include <algorithm>

/************************************************************************/
/* ==================================================================== */
/*                              CTable2Dataset                          */
/* ==================================================================== */
/************************************************************************/

class CTable2Dataset final : public RawDataset
{
    VSILFILE *fpImage;  // image data file.

    double adfGeoTransform[6];
    OGRSpatialReference m_oSRS{};

    CPL_DISALLOW_COPY_ASSIGN(CTable2Dataset)

    CPLErr Close() override;

  public:
    CTable2Dataset();
    ~CTable2Dataset() override;

    CPLErr SetGeoTransform(double *padfTransform) override;
    CPLErr GetGeoTransform(double *padfTransform) override;

    const OGRSpatialReference *GetSpatialRef() const override
    {
        return &m_oSRS;
    }

    static GDALDataset *Open(GDALOpenInfo *);
    static int Identify(GDALOpenInfo *);
    static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
                               int nBands, GDALDataType eType,
                               char **papszOptions);
};

/************************************************************************/
/* ==================================================================== */
/*                              CTable2Dataset                          */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                             CTable2Dataset()                          */
/************************************************************************/

CTable2Dataset::CTable2Dataset() : fpImage(nullptr)
{
    m_oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
    m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);

    memset(adfGeoTransform, 0, sizeof(adfGeoTransform));
}

/************************************************************************/
/*                            ~CTable2Dataset()                         */
/************************************************************************/

CTable2Dataset::~CTable2Dataset()

{
    CTable2Dataset::Close();
}

/************************************************************************/
/*                              Close()                                 */
/************************************************************************/

CPLErr CTable2Dataset::Close()
{
    CPLErr eErr = CE_None;
    if (nOpenFlags != OPEN_FLAGS_CLOSED)
    {
        if (CTable2Dataset::FlushCache(true) != CE_None)
            eErr = CE_Failure;

        if (fpImage != nullptr)
        {
            if (VSIFCloseL(fpImage) != 0)
            {
                CPLError(CE_Failure, CPLE_FileIO, "I/O error");
                eErr = CE_Failure;
            }
        }

        if (GDALPamDataset::Close() != CE_None)
            eErr = CE_Failure;
    }
    return eErr;
}

/************************************************************************/
/*                              Identify()                              */
/************************************************************************/

int CTable2Dataset::Identify(GDALOpenInfo *poOpenInfo)

{
    if (poOpenInfo->nHeaderBytes < 64)
        return FALSE;

    if (!STARTS_WITH_CI(
            reinterpret_cast<const char *>(poOpenInfo->pabyHeader + 0),
            "CTABLE V2"))
        return FALSE;

    return TRUE;
}

/************************************************************************/
/*                                Open()                                */
/************************************************************************/

GDALDataset *CTable2Dataset::Open(GDALOpenInfo *poOpenInfo)

{
    if (!Identify(poOpenInfo) || !poOpenInfo->fpL)
        return nullptr;

    /* -------------------------------------------------------------------- */
    /*      Create a corresponding GDALDataset.                             */
    /* -------------------------------------------------------------------- */
    auto poDS = std::make_unique<CTable2Dataset>();
    poDS->eAccess = poOpenInfo->eAccess;
    std::swap(poDS->fpImage, poOpenInfo->fpL);

    /* -------------------------------------------------------------------- */
    /*      Read the file header.                                           */
    /* -------------------------------------------------------------------- */

    CPL_IGNORE_RET_VAL(VSIFSeekL(poDS->fpImage, 0, SEEK_SET));

    char achHeader[160] = {'\0'};
    CPL_IGNORE_RET_VAL(VSIFReadL(achHeader, 1, 160, poDS->fpImage));
    achHeader[16 + 79] = '\0';

    CPLString osDescription = reinterpret_cast<const char *>(achHeader + 16);
    osDescription.Trim();
    poDS->SetMetadataItem("DESCRIPTION", osDescription);

    /* -------------------------------------------------------------------- */
    /*      Convert from LSB to local machine byte order.                   */
    /* -------------------------------------------------------------------- */
    CPL_LSBPTR64(achHeader + 96);
    CPL_LSBPTR64(achHeader + 104);
    CPL_LSBPTR64(achHeader + 112);
    CPL_LSBPTR64(achHeader + 120);
    CPL_LSBPTR32(achHeader + 128);
    CPL_LSBPTR32(achHeader + 132);

    /* -------------------------------------------------------------------- */
    /*      Extract size, and geotransform.                                 */
    /* -------------------------------------------------------------------- */
    int nRasterXSize, nRasterYSize;
    memcpy(&nRasterXSize, achHeader + 128, 4);
    memcpy(&nRasterYSize, achHeader + 132, 4);
    if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize) ||
        /* to avoid overflow in later -8 * nRasterXSize computation */
        nRasterXSize >= INT_MAX / 8)
    {
        return nullptr;
    }

    poDS->nRasterXSize = nRasterXSize;
    poDS->nRasterYSize = nRasterYSize;

    double adfValues[4];
    memcpy(adfValues, achHeader + 96, sizeof(double) * 4);

    for (int i = 0; i < 4; i++)
        adfValues[i] *= 180 / M_PI;  // Radians to degrees.

    poDS->adfGeoTransform[0] = adfValues[0] - adfValues[2] * 0.5;
    poDS->adfGeoTransform[1] = adfValues[2];
    poDS->adfGeoTransform[2] = 0.0;
    poDS->adfGeoTransform[3] =
        adfValues[1] + adfValues[3] * (nRasterYSize - 0.5);
    poDS->adfGeoTransform[4] = 0.0;
    poDS->adfGeoTransform[5] = -adfValues[3];

    /* -------------------------------------------------------------------- */
    /*      Setup the bands.                                                */
    /* -------------------------------------------------------------------- */
    auto poBand =
        RawRasterBand::Create(poDS.get(), 1, poDS->fpImage,
                              160 + 4 +
                                  static_cast<vsi_l_offset>(nRasterXSize) *
                                      (nRasterYSize - 1) * 2 * 4,
                              8, -8 * nRasterXSize, GDT_Float32,
                              RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
                              RawRasterBand::OwnFP::NO);
    if (!poBand)
        return nullptr;
    poBand->SetDescription("Latitude Offset (radians)");
    poDS->SetBand(1, std::move(poBand));

    poBand =
        RawRasterBand::Create(poDS.get(), 2, poDS->fpImage,
                              160 + static_cast<vsi_l_offset>(nRasterXSize) *
                                        (nRasterYSize - 1) * 2 * 4,
                              8, -8 * nRasterXSize, GDT_Float32,
                              RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
                              RawRasterBand::OwnFP::NO);
    if (!poBand)
        return nullptr;
    poBand->SetDescription("Longitude Offset (radians)");
    poBand->SetMetadataItem("positive_value", "west");
    poDS->SetBand(2, std::move(poBand));

    /* -------------------------------------------------------------------- */
    /*      Initialize any PAM information.                                 */
    /* -------------------------------------------------------------------- */
    poDS->SetDescription(poOpenInfo->pszFilename);
    poDS->TryLoadXML();

    /* -------------------------------------------------------------------- */
    /*      Check for overviews.                                            */
    /* -------------------------------------------------------------------- */
    poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);

    return poDS.release();
}

/************************************************************************/
/*                          GetGeoTransform()                           */
/************************************************************************/

CPLErr CTable2Dataset::GetGeoTransform(double *padfTransform)

{
    memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
    return CE_None;
}

/************************************************************************/
/*                          SetGeoTransform()                           */
/************************************************************************/

CPLErr CTable2Dataset::SetGeoTransform(double *padfTransform)

{
    if (eAccess == GA_ReadOnly)
    {
        CPLError(CE_Failure, CPLE_NoWriteAccess,
                 "Unable to update geotransform on readonly file.");
        return CE_Failure;
    }

    if (padfTransform[2] != 0.0 || padfTransform[4] != 0.0)
    {
        CPLError(
            CE_Failure, CPLE_AppDefined,
            "Rotated and sheared geotransforms not supported for CTable2.");
        return CE_Failure;
    }

    memcpy(adfGeoTransform, padfTransform, sizeof(double) * 6);

    /* -------------------------------------------------------------------- */
    /*      Update grid header.                                             */
    /* -------------------------------------------------------------------- */
    const double dfDegToRad = M_PI / 180.0;

    // read grid header
    CPL_IGNORE_RET_VAL(VSIFSeekL(fpImage, 0, SEEK_SET));

    char achHeader[160] = {'\0'};
    CPL_IGNORE_RET_VAL(VSIFReadL(achHeader, 1, sizeof(achHeader), fpImage));

    // lower left origin (longitude, center of pixel, radians)
    double dfValue =
        (adfGeoTransform[0] + adfGeoTransform[1] * 0.5) * dfDegToRad;
    CPL_LSBPTR64(&dfValue);
    memcpy(achHeader + 96, &dfValue, 8);

    // lower left origin (latitude, center of pixel, radians)
    dfValue = (adfGeoTransform[3] + adfGeoTransform[5] * (nRasterYSize - 0.5)) *
              dfDegToRad;
    CPL_LSBPTR64(&dfValue);
    memcpy(achHeader + 104, &dfValue, 8);

    // pixel width (radians)
    dfValue = adfGeoTransform[1] * dfDegToRad;
    CPL_LSBPTR64(&dfValue);
    memcpy(achHeader + 112, &dfValue, 8);

    // pixel height (radians)
    dfValue = adfGeoTransform[5] * -1 * dfDegToRad;
    CPL_LSBPTR64(&dfValue);
    memcpy(achHeader + 120, &dfValue, 8);

    // write grid header.
    CPL_IGNORE_RET_VAL(VSIFSeekL(fpImage, 0, SEEK_SET));
    CPL_IGNORE_RET_VAL(VSIFWriteL(achHeader, 1, sizeof(achHeader), fpImage));

    return CE_None;
}

/************************************************************************/
/*                               Create()                               */
/************************************************************************/

GDALDataset *CTable2Dataset::Create(const char *pszFilename, int nXSize,
                                    int nYSize, int /* nBandsIn */,
                                    GDALDataType eType, char **papszOptions)
{
    if (eType != GDT_Float32)
    {
        CPLError(
            CE_Failure, CPLE_AppDefined,
            "Attempt to create CTable2 file with unsupported data type '%s'.",
            GDALGetDataTypeName(eType));
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*      Try to open or create file.                                     */
    /* -------------------------------------------------------------------- */
    VSILFILE *fp = VSIFOpenL(pszFilename, "wb");

    if (fp == nullptr)
    {
        CPLError(CE_Failure, CPLE_OpenFailed,
                 "Attempt to create file `%s' failed.\n", pszFilename);
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*      Create a file header, with a defaulted georeferencing.          */
    /* -------------------------------------------------------------------- */
    char achHeader[160] = {'\0'};

    memset(achHeader, 0, sizeof(achHeader));

    memcpy(achHeader + 0, "CTABLE V2.0     ", 16);

    if (CSLFetchNameValue(papszOptions, "DESCRIPTION") != nullptr)
        strncpy(achHeader + 16, CSLFetchNameValue(papszOptions, "DESCRIPTION"),
                80);

    // lower left origin (longitude, center of pixel, radians)
    double dfValue = 0;
    CPL_LSBPTR64(&dfValue);
    memcpy(achHeader + 96, &dfValue, 8);

    // lower left origin (latitude, center of pixel, radians)
    dfValue = 0;
    CPL_LSBPTR64(&dfValue);
    memcpy(achHeader + 104, &dfValue, 8);

    // pixel width (radians)
    dfValue = 0.01 * M_PI / 180.0;
    CPL_LSBPTR64(&dfValue);
    memcpy(achHeader + 112, &dfValue, 8);

    // pixel height (radians)
    dfValue = 0.01 * M_PI / 180.0;
    CPL_LSBPTR64(&dfValue);
    memcpy(achHeader + 120, &dfValue, 8);

    // raster width in pixels
    int nValue32 = nXSize;
    CPL_LSBPTR32(&nValue32);
    memcpy(achHeader + 128, &nValue32, 4);

    // raster width in pixels
    nValue32 = nYSize;
    CPL_LSBPTR32(&nValue32);
    memcpy(achHeader + 132, &nValue32, 4);

    CPL_IGNORE_RET_VAL(VSIFWriteL(achHeader, 1, sizeof(achHeader), fp));

    /* -------------------------------------------------------------------- */
    /*      Write zeroed grid data.                                         */
    /* -------------------------------------------------------------------- */
    float *pafLine = static_cast<float *>(CPLCalloc(sizeof(float) * 2, nXSize));

    for (int i = 0; i < nYSize; i++)
    {
        if (static_cast<int>(
                VSIFWriteL(pafLine, sizeof(float) * 2, nXSize, fp)) != nXSize)
        {
            CPLError(CE_Failure, CPLE_FileIO,
                     "Write failed at line %d, perhaps the disk is full?", i);
            return nullptr;
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Cleanup and return.                                             */
    /* -------------------------------------------------------------------- */
    CPLFree(pafLine);

    if (VSIFCloseL(fp) != 0)
    {
        CPLError(CE_Failure, CPLE_FileIO, "I/O error");
        return nullptr;
    }

    return static_cast<GDALDataset *>(GDALOpen(pszFilename, GA_Update));
}

/************************************************************************/
/*                         GDALRegister_CTable2()                       */
/************************************************************************/

void GDALRegister_CTable2()

{
    if (GDALGetDriverByName("CTable2") != nullptr)
        return;

    GDALDriver *poDriver = new GDALDriver();

    poDriver->SetDescription("CTable2");
    poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "CTable2 Datum Grid Shift");
    poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");

    poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Float32");

    poDriver->pfnOpen = CTable2Dataset::Open;
    poDriver->pfnIdentify = CTable2Dataset::Identify;
    poDriver->pfnCreate = CTable2Dataset::Create;

    GetGDALDriverManager()->RegisterDriver(poDriver);
}