/src/gdal/frmts/raw/byndataset.cpp

Source (jump to first uncovered line)
/******************************************************************************
 *
 * Project:  Natural Resources Canada's Geoid BYN file format
 * Purpose:  Implementation of BYN format
 * Author:   Ivan Lucena, ivan.lucena@outlook.com
 *
 ******************************************************************************
 * Copyright (c) 2018, Ivan Lucena
 * Copyright (c) 2018, Even Rouault
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_port.h"
#include "byndataset.h"
#include "rawdataset.h"

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

#include <algorithm>
#include <cstdlib>
#include <limits>

// Specification at
// https://www.nrcan.gc.ca/sites/www.nrcan.gc.ca/files/earthsciences/pdf/gpshgrid_e.pdf

const static BYNEllipsoids EllipsoidTable[] = {
    {"GRS80", 6378137.0, 298.257222101},
    {"WGS84", 6378137.0, 298.257223564},
    {"ALT1", 6378136.3, 298.256415099},
    {"GRS67", 6378160.0, 298.247167427},
    {"ELLIP1", 6378136.46, 298.256415099},
    {"ALT2", 6378136.3, 298.257},
    {"ELLIP2", 6378136.0, 298.257},
    {"CLARKE 1866", 6378206.4, 294.9786982}};

/************************************************************************/
/*                            BYNRasterBand()                           */
/************************************************************************/

BYNRasterBand::BYNRasterBand(GDALDataset *poDSIn, int nBandIn,
                             VSILFILE *fpRawIn, vsi_l_offset nImgOffsetIn,
                             int nPixelOffsetIn, int nLineOffsetIn,
                             GDALDataType eDataTypeIn, int bNativeOrderIn)
    : RawRasterBand(poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
                    nLineOffsetIn, eDataTypeIn, bNativeOrderIn,
                    RawRasterBand::OwnFP::NO)
{
}

/************************************************************************/
/*                           ~BYNRasterBand()                           */
/************************************************************************/

BYNRasterBand::~BYNRasterBand()
{
}

/************************************************************************/
/*                           GetNoDataValue()                           */
/************************************************************************/

double BYNRasterBand::GetNoDataValue(int *pbSuccess)
{
    if (pbSuccess)
        *pbSuccess = TRUE;
    int bSuccess = FALSE;
    double dfNoData = GDALPamRasterBand::GetNoDataValue(&bSuccess);
    if (bSuccess)
    {
        return dfNoData;
    }
    const double dfFactor =
        reinterpret_cast<BYNDataset *>(poDS)->hHeader.dfFactor;
    return eDataType == GDT_Int16 ? 32767.0 : 9999.0 * dfFactor;
}

/************************************************************************/
/*                              GetScale()                              */
/************************************************************************/

double BYNRasterBand::GetScale(int *pbSuccess)
{
    if (pbSuccess != nullptr)
        *pbSuccess = TRUE;
    const double dfFactor =
        reinterpret_cast<BYNDataset *>(poDS)->hHeader.dfFactor;
    return (dfFactor != 0.0) ? 1.0 / dfFactor : 0.0;
}

/************************************************************************/
/* ==================================================================== */
/*                              BYNDataset                              */
/* ==================================================================== */
/************************************************************************/

BYNDataset::BYNDataset()
    : fpImage(nullptr), hHeader{0, 0, 0, 0, 0, 0,   0,   0, 0.0, 0,   0, 0,
                                0, 0, 0, 0, 0, 0.0, 0.0, 0, 0,   0.0, 0}
{
    m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    adfGeoTransform[0] = 0.0;
    adfGeoTransform[1] = 1.0;
    adfGeoTransform[2] = 0.0;
    adfGeoTransform[3] = 0.0;
    adfGeoTransform[4] = 0.0;
    adfGeoTransform[5] = 1.0;
}

/************************************************************************/
/*                            ~BYNDataset()                             */
/************************************************************************/

BYNDataset::~BYNDataset()

{
    BYNDataset::Close();
}

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

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

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

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

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

int BYNDataset::Identify(GDALOpenInfo *poOpenInfo)

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

/* -------------------------------------------------------------------- */
/*      Check file extension (.byn/.err)                                */
/* -------------------------------------------------------------------- */
#ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
    const char *pszFileExtension = poOpenInfo->osExtension.c_str();

    if (!EQUAL(pszFileExtension, "byn") && !EQUAL(pszFileExtension, "err"))
    {
        return FALSE;
    }
#endif

    /* -------------------------------------------------------------------- */
    /*      Check some value's ranges on header                             */
    /* -------------------------------------------------------------------- */

    BYNHeader hHeader = {0, 0, 0, 0, 0, 0,   0,   0, 0.0, 0,   0, 0,
                         0, 0, 0, 0, 0, 0.0, 0.0, 0, 0,   0.0, 0};

    buffer2header(poOpenInfo->pabyHeader, &hHeader);

    if (hHeader.nGlobal < 0 || hHeader.nGlobal > 1 || hHeader.nType < 0 ||
        hHeader.nType > 9 || (hHeader.nSizeOf != 2 && hHeader.nSizeOf != 4) ||
        hHeader.nVDatum < 0 || hHeader.nVDatum > 3 || hHeader.nDescrip < 0 ||
        hHeader.nDescrip > 3 || hHeader.nSubType < 0 || hHeader.nSubType > 9 ||
        hHeader.nDatum < 0 || hHeader.nDatum > 1 || hHeader.nEllipsoid < 0 ||
        hHeader.nEllipsoid > 7 || hHeader.nByteOrder < 0 ||
        hHeader.nByteOrder > 1 || hHeader.nScale < 0 || hHeader.nScale > 1)
        return FALSE;

#if 0
    // We have disabled those checks as invalid values are often found in some
    // datasets, such as http://s3.microsurvey.com/os/fieldgenius/geoids/Lithuania.zip
    // We don't use those fields, so we may just ignore them.
    if((hHeader.nTideSys   < 0 || hHeader.nTideSys   > 2 ||
        hHeader.nPtType    < 0 || hHeader.nPtType    > 1 ))
    {
        // Some datasets use 0xCC as a marker for invalidity for
        // records starting from Geopotential Wo
        for( int i = 52; i < 78; i++ )
        {
            if( poOpenInfo->pabyHeader[i] != 0xCC )
                return FALSE;
        }
    }
#endif

    if (hHeader.nScale == 0)
    {
        if ((std::abs(static_cast<GIntBig>(hHeader.nSouth) -
                      (hHeader.nDLat / 2)) > BYN_MAX_LAT) ||
            (std::abs(static_cast<GIntBig>(hHeader.nNorth) +
                      (hHeader.nDLat / 2)) > BYN_MAX_LAT) ||
            (std::abs(static_cast<GIntBig>(hHeader.nWest) -
                      (hHeader.nDLon / 2)) > BYN_MAX_LON) ||
            (std::abs(static_cast<GIntBig>(hHeader.nEast) +
                      (hHeader.nDLon / 2)) > BYN_MAX_LON))
            return FALSE;
    }
    else
    {
        if ((std::abs(static_cast<GIntBig>(hHeader.nSouth) -
                      (hHeader.nDLat / 2)) > BYN_MAX_LAT_SCL) ||
            (std::abs(static_cast<GIntBig>(hHeader.nNorth) +
                      (hHeader.nDLat / 2)) > BYN_MAX_LAT_SCL) ||
            (std::abs(static_cast<GIntBig>(hHeader.nWest) -
                      (hHeader.nDLon / 2)) > BYN_MAX_LON_SCL) ||
            (std::abs(static_cast<GIntBig>(hHeader.nEast) +
                      (hHeader.nDLon / 2)) > BYN_MAX_LON_SCL))
            return FALSE;
    }

    return TRUE;
}

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

GDALDataset *BYNDataset::Open(GDALOpenInfo *poOpenInfo)

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

    /* -------------------------------------------------------------------- */
    /*      Create a corresponding GDALDataset.                             */
    /* -------------------------------------------------------------------- */

    auto poDS = std::make_unique<BYNDataset>();

    poDS->eAccess = poOpenInfo->eAccess;
    std::swap(poDS->fpImage, poOpenInfo->fpL);

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

    buffer2header(poOpenInfo->pabyHeader, &poDS->hHeader);

    /********************************/
    /* Scale boundaries and spacing */
    /********************************/

    double dfSouth = poDS->hHeader.nSouth;
    double dfNorth = poDS->hHeader.nNorth;
    double dfWest = poDS->hHeader.nWest;
    double dfEast = poDS->hHeader.nEast;
    double dfDLat = poDS->hHeader.nDLat;
    double dfDLon = poDS->hHeader.nDLon;

    if (poDS->hHeader.nScale == 1)
    {
        dfSouth *= BYN_SCALE;
        dfNorth *= BYN_SCALE;
        dfWest *= BYN_SCALE;
        dfEast *= BYN_SCALE;
        dfDLat *= BYN_SCALE;
        dfDLon *= BYN_SCALE;
    }

    /******************************/
    /* Calculate rows and columns */
    /******************************/

    double dfXSize = -1;
    double dfYSize = -1;

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

    if (dfDLat != 0.0 && dfDLon != 0.0)
    {
        dfXSize = ((dfEast - dfWest + 1.0) / dfDLon) + 1.0;
        dfYSize = ((dfNorth - dfSouth + 1.0) / dfDLat) + 1.0;
    }

    if (dfXSize > 0.0 && dfXSize < std::numeric_limits<double>::max() &&
        dfYSize > 0.0 && dfYSize < std::numeric_limits<double>::max())
    {
        poDS->nRasterXSize = static_cast<GInt32>(dfXSize);
        poDS->nRasterYSize = static_cast<GInt32>(dfYSize);
    }

    if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
    {
        return nullptr;
    }

    /*****************************/
    /* Build GeoTransform matrix */
    /*****************************/

    poDS->adfGeoTransform[0] = (dfWest - (dfDLon / 2.0)) / 3600.0;
    poDS->adfGeoTransform[1] = dfDLon / 3600.0;
    poDS->adfGeoTransform[2] = 0.0;
    poDS->adfGeoTransform[3] = (dfNorth + (dfDLat / 2.0)) / 3600.0;
    poDS->adfGeoTransform[4] = 0.0;
    poDS->adfGeoTransform[5] = -1 * dfDLat / 3600.0;

    /*********************/
    /* Set data type     */
    /*********************/

    GDALDataType eDT = GDT_Unknown;

    if (poDS->hHeader.nSizeOf == 2)
        eDT = GDT_Int16;
    else if (poDS->hHeader.nSizeOf == 4)
        eDT = GDT_Int32;
    else
    {
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*      Create band information object.                                 */
    /* -------------------------------------------------------------------- */

    const int nDTSize = GDALGetDataTypeSizeBytes(eDT);

    int bIsLSB = poDS->hHeader.nByteOrder == 1 ? 1 : 0;

    auto poBand = std::make_unique<BYNRasterBand>(
        poDS.get(), 1, poDS->fpImage, BYN_HDR_SZ, nDTSize,
        poDS->nRasterXSize * nDTSize, eDT, CPL_IS_LSB == bIsLSB);
    if (!poBand->IsValid())
        return nullptr;
    poDS->SetBand(1, 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 BYNDataset::GetGeoTransform(double *padfTransform)

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

/************************************************************************/
/*                          GetSpatialRef()                             */
/************************************************************************/

const OGRSpatialReference *BYNDataset::GetSpatialRef() const

{
    if (!m_oSRS.IsEmpty())
        return &m_oSRS;

    /* Try to use a prefefined EPSG compound CS */

    if (hHeader.nDatum == 1 && hHeader.nVDatum == 2)
    {
        m_oSRS.importFromEPSG(BYN_DATUM_1_VDATUM_2);
        return &m_oSRS;
    }

    /* Build the GEOGCS based on Datum ( or Ellipsoid )*/

    bool bNoGeogCS = false;

    if (hHeader.nDatum == 0)
        m_oSRS.importFromEPSG(BYN_DATUM_0);
    else if (hHeader.nDatum == 1)
        m_oSRS.importFromEPSG(BYN_DATUM_1);
    else
    {
        /* Build GEOGCS based on Ellipsoid (Table 3) */

        if (hHeader.nEllipsoid > -1 &&
            hHeader.nEllipsoid <
                static_cast<GInt16>(CPL_ARRAYSIZE(EllipsoidTable)))
            m_oSRS.SetGeogCS(
                CPLSPrintf("BYN Ellipsoid(%d)", hHeader.nEllipsoid),
                "Unspecified", EllipsoidTable[hHeader.nEllipsoid].pszName,
                EllipsoidTable[hHeader.nEllipsoid].dfSemiMajor,
                EllipsoidTable[hHeader.nEllipsoid].dfInvFlattening);
        else
            bNoGeogCS = true;
    }

    /* Build the VERT_CS based on VDatum */

    OGRSpatialReference oSRSComp;
    OGRSpatialReference oSRSVert;

    int nVertCS = 0;

    if (hHeader.nVDatum == 1)
        nVertCS = BYN_VDATUM_1;
    else if (hHeader.nVDatum == 2)
        nVertCS = BYN_VDATUM_2;
    else if (hHeader.nVDatum == 3)
        nVertCS = BYN_VDATUM_3;
    else
    {
        /* Return GEOGCS ( .err files ) */

        if (bNoGeogCS)
            return nullptr;

        return &m_oSRS;
    }

    oSRSVert.importFromEPSG(nVertCS);

    /* Create CPMPD_CS with GEOGCS and VERT_CS */

    if (oSRSComp.SetCompoundCS(CPLSPrintf("BYN Datum(%d) & VDatum(%d)",
                                          hHeader.nDatum, hHeader.nDatum),
                               &m_oSRS, &oSRSVert) == CE_None)
    {
        /* Return COMPD_CS with GEOGCS and VERT_CS */

        m_oSRS = std::move(oSRSComp);
        m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
        return &m_oSRS;
    }

    return nullptr;
}

/*----------------------------------------------------------------------*/
/*                           buffer2header()                            */
/*----------------------------------------------------------------------*/

void BYNDataset::buffer2header(const GByte *pabyBuf, BYNHeader *pohHeader)

{
    memcpy(&pohHeader->nSouth, pabyBuf, 4);
    memcpy(&pohHeader->nNorth, pabyBuf + 4, 4);
    memcpy(&pohHeader->nWest, pabyBuf + 8, 4);
    memcpy(&pohHeader->nEast, pabyBuf + 12, 4);
    memcpy(&pohHeader->nDLat, pabyBuf + 16, 2);
    memcpy(&pohHeader->nDLon, pabyBuf + 18, 2);
    memcpy(&pohHeader->nGlobal, pabyBuf + 20, 2);
    memcpy(&pohHeader->nType, pabyBuf + 22, 2);
    memcpy(&pohHeader->dfFactor, pabyBuf + 24, 8);
    memcpy(&pohHeader->nSizeOf, pabyBuf + 32, 2);
    memcpy(&pohHeader->nVDatum, pabyBuf + 34, 2);
    memcpy(&pohHeader->nDescrip, pabyBuf + 40, 2);
    memcpy(&pohHeader->nSubType, pabyBuf + 42, 2);
    memcpy(&pohHeader->nDatum, pabyBuf + 44, 2);
    memcpy(&pohHeader->nEllipsoid, pabyBuf + 46, 2);
    memcpy(&pohHeader->nByteOrder, pabyBuf + 48, 2);
    memcpy(&pohHeader->nScale, pabyBuf + 50, 2);
    memcpy(&pohHeader->dfWo, pabyBuf + 52, 8);
    memcpy(&pohHeader->dfGM, pabyBuf + 60, 8);
    memcpy(&pohHeader->nTideSys, pabyBuf + 68, 2);
    memcpy(&pohHeader->nRealiz, pabyBuf + 70, 2);
    memcpy(&pohHeader->dEpoch, pabyBuf + 72, 4);
    memcpy(&pohHeader->nPtType, pabyBuf + 76, 2);

#if defined(CPL_MSB)
    CPL_LSBPTR32(&pohHeader->nSouth);
    CPL_LSBPTR32(&pohHeader->nNorth);
    CPL_LSBPTR32(&pohHeader->nWest);
    CPL_LSBPTR32(&pohHeader->nEast);
    CPL_LSBPTR16(&pohHeader->nDLat);
    CPL_LSBPTR16(&pohHeader->nDLon);
    CPL_LSBPTR16(&pohHeader->nGlobal);
    CPL_LSBPTR16(&pohHeader->nType);
    CPL_LSBPTR64(&pohHeader->dfFactor);
    CPL_LSBPTR16(&pohHeader->nSizeOf);
    CPL_LSBPTR16(&pohHeader->nVDatum);
    CPL_LSBPTR16(&pohHeader->nDescrip);
    CPL_LSBPTR16(&pohHeader->nSubType);
    CPL_LSBPTR16(&pohHeader->nDatum);
    CPL_LSBPTR16(&pohHeader->nEllipsoid);
    CPL_LSBPTR16(&pohHeader->nByteOrder);
    CPL_LSBPTR16(&pohHeader->nScale);
    CPL_LSBPTR64(&pohHeader->dfWo);
    CPL_LSBPTR64(&pohHeader->dfGM);
    CPL_LSBPTR16(&pohHeader->nTideSys);
    CPL_LSBPTR16(&pohHeader->nRealiz);
    CPL_LSBPTR32(&pohHeader->dEpoch);
    CPL_LSBPTR16(&pohHeader->nPtType);
#endif

#if DEBUG
    CPLDebug("BYN", "South         = %d", pohHeader->nSouth);
    CPLDebug("BYN", "North         = %d", pohHeader->nNorth);
    CPLDebug("BYN", "West          = %d", pohHeader->nWest);
    CPLDebug("BYN", "East          = %d", pohHeader->nEast);
    CPLDebug("BYN", "DLat          = %d", pohHeader->nDLat);
    CPLDebug("BYN", "DLon          = %d", pohHeader->nDLon);
    CPLDebug("BYN", "DGlobal       = %d", pohHeader->nGlobal);
    CPLDebug("BYN", "DType         = %d", pohHeader->nType);
    CPLDebug("BYN", "Factor        = %f", pohHeader->dfFactor);
    CPLDebug("BYN", "SizeOf        = %d", pohHeader->nSizeOf);
    CPLDebug("BYN", "VDatum        = %d", pohHeader->nVDatum);
    CPLDebug("BYN", "Data          = %d", pohHeader->nDescrip);
    CPLDebug("BYN", "SubType       = %d", pohHeader->nSubType);
    CPLDebug("BYN", "Datum         = %d", pohHeader->nDatum);
    CPLDebug("BYN", "Ellipsoid     = %d", pohHeader->nEllipsoid);
    CPLDebug("BYN", "ByteOrder     = %d", pohHeader->nByteOrder);
    CPLDebug("BYN", "Scale         = %d", pohHeader->nScale);
    CPLDebug("BYN", "Wo            = %f", pohHeader->dfWo);
    CPLDebug("BYN", "GM            = %f", pohHeader->dfGM);
    CPLDebug("BYN", "TideSystem    = %d", pohHeader->nTideSys);
    CPLDebug("BYN", "RefRealzation = %d", pohHeader->nRealiz);
    CPLDebug("BYN", "Epoch         = %f", pohHeader->dEpoch);
    CPLDebug("BYN", "PtType        = %d", pohHeader->nPtType);
#endif
}

/************************************************************************/
/*                          GDALRegister_BYN()                          */
/************************************************************************/

void GDALRegister_BYN()

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

    GDALDriver *poDriver = new GDALDriver();

    poDriver->SetDescription("BYN");
    poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
                              "Natural Resources Canada's Geoid");
    poDriver->SetMetadataItem(GDAL_DMD_EXTENSIONS, "byn err");
    poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/byn.html");

    poDriver->pfnOpen = BYNDataset::Open;
    poDriver->pfnIdentify = BYNDataset::Identify;

    GetGDALDriverManager()->RegisterDriver(poDriver);
}

Coverage Report

Created: 2025-06-09 08:44

Line	Count	Source (jump to first uncovered line)
1		/******************************************************************************
2		*
3		* Project: Natural Resources Canada's Geoid BYN file format
4		* Purpose: Implementation of BYN format
5		* Author: Ivan Lucena, ivan.lucena@outlook.com
6		*
7		******************************************************************************
8		* Copyright (c) 2018, Ivan Lucena
9		* Copyright (c) 2018, Even Rouault
10		*
11		* SPDX-License-Identifier: MIT
12		****************************************************************************/
13
14		#include "cpl_port.h"
15		#include "byndataset.h"
16		#include "rawdataset.h"
17
18		#include "cpl_string.h"
19		#include "gdal_frmts.h"
20		#include "ogr_spatialref.h"
21		#include "ogr_srs_api.h"
22
23		#include <algorithm>
24		#include <cstdlib>
25		#include <limits>
26
27		// Specification at
28		// https://www.nrcan.gc.ca/sites/www.nrcan.gc.ca/files/earthsciences/pdf/gpshgrid_e.pdf
29
30		const static BYNEllipsoids EllipsoidTable[] = {
31		{"GRS80", 6378137.0, 298.257222101},
32		{"WGS84", 6378137.0, 298.257223564},
33		{"ALT1", 6378136.3, 298.256415099},
34		{"GRS67", 6378160.0, 298.247167427},
35		{"ELLIP1", 6378136.46, 298.256415099},
36		{"ALT2", 6378136.3, 298.257},
37		{"ELLIP2", 6378136.0, 298.257},
38		{"CLARKE 1866", 6378206.4, 294.9786982}};
39
40		/************************************************************************/
41		/* BYNRasterBand() */
42		/************************************************************************/
43
44		BYNRasterBand::BYNRasterBand(GDALDataset *poDSIn, int nBandIn,
45		VSILFILE *fpRawIn, vsi_l_offset nImgOffsetIn,
46		int nPixelOffsetIn, int nLineOffsetIn,
47		GDALDataType eDataTypeIn, int bNativeOrderIn)
48	650	: RawRasterBand(poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
49	650	nLineOffsetIn, eDataTypeIn, bNativeOrderIn,
50	650	RawRasterBand::OwnFP::NO)
51	650	{
52	650	}
53
54		/************************************************************************/
55		/* ~BYNRasterBand() */
56		/************************************************************************/
57
58		BYNRasterBand::~BYNRasterBand()
59	650	{
60	650	}
61
62		/************************************************************************/
63		/* GetNoDataValue() */
64		/************************************************************************/
65
66		double BYNRasterBand::GetNoDataValue(int *pbSuccess)
67	182k	{
68	182k	if (pbSuccess)
69	182k	*pbSuccess = TRUE;
70	182k	int bSuccess = FALSE;
71	182k	double dfNoData = GDALPamRasterBand::GetNoDataValue(&bSuccess);
72	182k	if (bSuccess)
73	0	{
74	0	return dfNoData;
75	0	}
76	182k	const double dfFactor =
77	182k	reinterpret_cast<BYNDataset *>(poDS)->hHeader.dfFactor;
78	182k	return eDataType == GDT_Int16 ? 32767.0 : 9999.0 * dfFactor;
79	182k	}
80
81		/************************************************************************/
82		/* GetScale() */
83		/************************************************************************/
84
85		double BYNRasterBand::GetScale(int *pbSuccess)
86	16.2k	{
87	16.2k	if (pbSuccess != nullptr)
88	75	*pbSuccess = TRUE;
89	16.2k	const double dfFactor =
90	16.2k	reinterpret_cast<BYNDataset *>(poDS)->hHeader.dfFactor;
91	16.2k	return (dfFactor != 0.0) ? 1.0 / dfFactor : 0.0;
92	16.2k	}
93
94		/************************************************************************/
95		/* ==================================================================== */
96		/* BYNDataset */
97		/* ==================================================================== */
98		/************************************************************************/
99
100		BYNDataset::BYNDataset()
101	658	: fpImage(nullptr), hHeader{0, 0, 0, 0, 0, 0, 0, 0, 0.0, 0, 0, 0,
102	658	0, 0, 0, 0, 0, 0.0, 0.0, 0, 0, 0.0, 0}
103	658	{
104	658	m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
105	658	adfGeoTransform[0] = 0.0;
106	658	adfGeoTransform[1] = 1.0;
107	658	adfGeoTransform[2] = 0.0;
108	658	adfGeoTransform[3] = 0.0;
109	658	adfGeoTransform[4] = 0.0;
110	658	adfGeoTransform[5] = 1.0;
111	658	}
112
113		/************************************************************************/
114		/* ~BYNDataset() */
115		/************************************************************************/
116
117		BYNDataset::~BYNDataset()
118
119	658	{
120	658	BYNDataset::Close();
121	658	}
122
123		/************************************************************************/
124		/* Close() */
125		/************************************************************************/
126
127		CPLErr BYNDataset::Close()
128	1.30k	{
129	1.30k	CPLErr eErr = CE_None;
130	1.30k	if (nOpenFlags != OPEN_FLAGS_CLOSED)
131	658	{
132	658	if (BYNDataset::FlushCache(true) != CE_None)
133	13	eErr = CE_Failure;
134
135	658	if (fpImage != nullptr)
136	658	{
137	658	if (VSIFCloseL(fpImage) != 0)
138	0	{
139	0	eErr = CE_Failure;
140	0	CPLError(CE_Failure, CPLE_FileIO, "I/O error");
141	0	}
142	658	}
143
144	658	if (GDALPamDataset::Close() != CE_None)
145	0	eErr = CE_Failure;
146	658	}
147	1.30k	return eErr;
148	1.30k	}
149
150		/************************************************************************/
151		/* Identify() */
152		/************************************************************************/
153
154		int BYNDataset::Identify(GDALOpenInfo *poOpenInfo)
155
156	860k	{
157	860k	if (poOpenInfo->nHeaderBytes < BYN_HDR_SZ)
158	787k	return FALSE;
159
160		/* -------------------------------------------------------------------- */
161		/* Check file extension (.byn/.err) */
162		/* -------------------------------------------------------------------- */
163		#ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
164		const char *pszFileExtension = poOpenInfo->osExtension.c_str();
165
166		if (!EQUAL(pszFileExtension, "byn") && !EQUAL(pszFileExtension, "err"))
167		{
168		return FALSE;
169		}
170		#endif
171
172		/* -------------------------------------------------------------------- */
173		/* Check some value's ranges on header */
174		/* -------------------------------------------------------------------- */
175
176	73.1k	BYNHeader hHeader = {0, 0, 0, 0, 0, 0, 0, 0, 0.0, 0, 0, 0,
177	73.1k	0, 0, 0, 0, 0, 0.0, 0.0, 0, 0, 0.0, 0};
178
179	73.1k	buffer2header(poOpenInfo->pabyHeader, &hHeader);
180
181	73.1k	if (hHeader.nGlobal < 0 \|\| hHeader.nGlobal > 1 \|\| hHeader.nType < 0 \|\|
182	73.1k	hHeader.nType > 9 \|\| (hHeader.nSizeOf != 2 && hHeader.nSizeOf != 4) \|\|
183	73.1k	hHeader.nVDatum < 0 \|\| hHeader.nVDatum > 3 \|\| hHeader.nDescrip < 0 \|\|
184	73.1k	hHeader.nDescrip > 3 \|\| hHeader.nSubType < 0 \|\| hHeader.nSubType > 9 \|\|
185	73.1k	hHeader.nDatum < 0 \|\| hHeader.nDatum > 1 \|\| hHeader.nEllipsoid < 0 \|\|
186	73.1k	hHeader.nEllipsoid > 7 \|\| hHeader.nByteOrder < 0 \|\|
187	73.1k	hHeader.nByteOrder > 1 \|\| hHeader.nScale < 0 \|\| hHeader.nScale > 1)
188	71.7k	return FALSE;
189
190		#if 0
191		// We have disabled those checks as invalid values are often found in some
192		// datasets, such as http://s3.microsurvey.com/os/fieldgenius/geoids/Lithuania.zip
193		// We don't use those fields, so we may just ignore them.
194		if((hHeader.nTideSys < 0 \|\| hHeader.nTideSys > 2 \|\|
195		hHeader.nPtType < 0 \|\| hHeader.nPtType > 1 ))
196		{
197		// Some datasets use 0xCC as a marker for invalidity for
198		// records starting from Geopotential Wo
199		for( int i = 52; i < 78; i++ )
200		{
201		if( poOpenInfo->pabyHeader[i] != 0xCC )
202		return FALSE;
203		}
204		}
205		#endif
206
207	1.34k	if (hHeader.nScale == 0)
208	1.34k	{
209	1.34k	if ((std::abs(static_cast<GIntBig>(hHeader.nSouth) -
210	1.34k	(hHeader.nDLat / 2)) > BYN_MAX_LAT) \|\|
211	1.34k	(std::abs(static_cast<GIntBig>(hHeader.nNorth) +
212	1.32k	(hHeader.nDLat / 2)) > BYN_MAX_LAT) \|\|
213	1.34k	(std::abs(static_cast<GIntBig>(hHeader.nWest) -
214	1.32k	(hHeader.nDLon / 2)) > BYN_MAX_LON) \|\|
215	1.34k	(std::abs(static_cast<GIntBig>(hHeader.nEast) +
216	1.31k	(hHeader.nDLon / 2)) > BYN_MAX_LON))
217	24	return FALSE;
218	1.34k	}
219	2	else
220	2	{
221	2	if ((std::abs(static_cast<GIntBig>(hHeader.nSouth) -
222	2	(hHeader.nDLat / 2)) > BYN_MAX_LAT_SCL) \|\|
223	2	(std::abs(static_cast<GIntBig>(hHeader.nNorth) +
224	1	(hHeader.nDLat / 2)) > BYN_MAX_LAT_SCL) \|\|
225	2	(std::abs(static_cast<GIntBig>(hHeader.nWest) -
226	0	(hHeader.nDLon / 2)) > BYN_MAX_LON_SCL) \|\|
227	2	(std::abs(static_cast<GIntBig>(hHeader.nEast) +
228	0	(hHeader.nDLon / 2)) > BYN_MAX_LON_SCL))
229	2	return FALSE;
230	2	}
231
232	1.31k	return TRUE;
233	1.34k	}
234
235		/************************************************************************/
236		/* Open() */
237		/************************************************************************/
238
239		GDALDataset BYNDataset::Open(GDALOpenInfo poOpenInfo)
240
241	658	{
242	658	if (!Identify(poOpenInfo) \|\| poOpenInfo->fpL == nullptr \|\|
243	658	poOpenInfo->eAccess == GA_Update)
244	0	return nullptr;
245
246		/* -------------------------------------------------------------------- */
247		/* Create a corresponding GDALDataset. */
248		/* -------------------------------------------------------------------- */
249
250	658	auto poDS = std::make_unique<BYNDataset>();
251
252	658	poDS->eAccess = poOpenInfo->eAccess;
253	658	std::swap(poDS->fpImage, poOpenInfo->fpL);
254
255		/* -------------------------------------------------------------------- */
256		/* Read the header. */
257		/* -------------------------------------------------------------------- */
258
259	658	buffer2header(poOpenInfo->pabyHeader, &poDS->hHeader);
260
261		/********************************/
262		/* Scale boundaries and spacing */
263		/********************************/
264
265	658	double dfSouth = poDS->hHeader.nSouth;
266	658	double dfNorth = poDS->hHeader.nNorth;
267	658	double dfWest = poDS->hHeader.nWest;
268	658	double dfEast = poDS->hHeader.nEast;
269	658	double dfDLat = poDS->hHeader.nDLat;
270	658	double dfDLon = poDS->hHeader.nDLon;
271
272	658	if (poDS->hHeader.nScale == 1)
273	0	{
274	0	dfSouth *= BYN_SCALE;
275	0	dfNorth *= BYN_SCALE;
276	0	dfWest *= BYN_SCALE;
277	0	dfEast *= BYN_SCALE;
278	0	dfDLat *= BYN_SCALE;
279	0	dfDLon *= BYN_SCALE;
280	0	}
281
282		/******************************/
283		/* Calculate rows and columns */
284		/******************************/
285
286	658	double dfXSize = -1;
287	658	double dfYSize = -1;
288
289	658	poDS->nRasterXSize = -1;
290	658	poDS->nRasterYSize = -1;
291
292	658	if (dfDLat != 0.0 && dfDLon != 0.0)
293	654	{
294	654	dfXSize = ((dfEast - dfWest + 1.0) / dfDLon) + 1.0;
295	654	dfYSize = ((dfNorth - dfSouth + 1.0) / dfDLat) + 1.0;
296	654	}
297
298	658	if (dfXSize > 0.0 && dfXSize < std::numeric_limits<double>::max() &&
299	658	dfYSize > 0.0 && dfYSize < std::numeric_limits<double>::max())
300	650	{
301	650	poDS->nRasterXSize = static_cast<GInt32>(dfXSize);
302	650	poDS->nRasterYSize = static_cast<GInt32>(dfYSize);
303	650	}
304
305	658	if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
306	8	{
307	8	return nullptr;
308	8	}
309
310		/*****************************/
311		/* Build GeoTransform matrix */
312		/*****************************/
313
314	650	poDS->adfGeoTransform[0] = (dfWest - (dfDLon / 2.0)) / 3600.0;
315	650	poDS->adfGeoTransform[1] = dfDLon / 3600.0;
316	650	poDS->adfGeoTransform[2] = 0.0;
317	650	poDS->adfGeoTransform[3] = (dfNorth + (dfDLat / 2.0)) / 3600.0;
318	650	poDS->adfGeoTransform[4] = 0.0;
319	650	poDS->adfGeoTransform[5] = -1 * dfDLat / 3600.0;
320
321		/*********************/
322		/* Set data type */
323		/*********************/
324
325	650	GDALDataType eDT = GDT_Unknown;
326
327	650	if (poDS->hHeader.nSizeOf == 2)
328	200	eDT = GDT_Int16;
329	450	else if (poDS->hHeader.nSizeOf == 4)
330	450	eDT = GDT_Int32;
331	0	else
332	0	{
333	0	return nullptr;
334	0	}
335
336		/* -------------------------------------------------------------------- */
337		/* Create band information object. */
338		/* -------------------------------------------------------------------- */
339
340	650	const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
341
342	650	int bIsLSB = poDS->hHeader.nByteOrder == 1 ? 1 : 0;
343
344	650	auto poBand = std::make_unique<BYNRasterBand>(
345	650	poDS.get(), 1, poDS->fpImage, BYN_HDR_SZ, nDTSize,
346	650	poDS->nRasterXSize * nDTSize, eDT, CPL_IS_LSB == bIsLSB);
347	650	if (!poBand->IsValid())
348	0	return nullptr;
349	650	poDS->SetBand(1, std::move(poBand));
350
351		/* -------------------------------------------------------------------- */
352		/* Initialize any PAM information. */
353		/* -------------------------------------------------------------------- */
354
355	650	poDS->SetDescription(poOpenInfo->pszFilename);
356	650	poDS->TryLoadXML();
357
358		/* -------------------------------------------------------------------- */
359		/* Check for overviews. */
360		/* -------------------------------------------------------------------- */
361
362	650	poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
363
364	650	return poDS.release();
365	650	}
366
367		/************************************************************************/
368		/* GetGeoTransform() */
369		/************************************************************************/
370
371		CPLErr BYNDataset::GetGeoTransform(double *padfTransform)
372
373	786	{
374	786	memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
375	786	return CE_None;
376	786	}
377
378		/************************************************************************/
379		/* GetSpatialRef() */
380		/************************************************************************/
381
382		const OGRSpatialReference *BYNDataset::GetSpatialRef() const
383
384	693	{
385	693	if (!m_oSRS.IsEmpty())
386	74	return &m_oSRS;
387
388		/* Try to use a prefefined EPSG compound CS */
389
390	619	if (hHeader.nDatum == 1 && hHeader.nVDatum == 2)
391	18	{
392	18	m_oSRS.importFromEPSG(BYN_DATUM_1_VDATUM_2);
393	18	return &m_oSRS;
394	18	}
395
396		/* Build the GEOGCS based on Datum ( or Ellipsoid )*/
397
398	601	bool bNoGeogCS = false;
399
400	601	if (hHeader.nDatum == 0)
401	270	m_oSRS.importFromEPSG(BYN_DATUM_0);
402	331	else if (hHeader.nDatum == 1)
403	331	m_oSRS.importFromEPSG(BYN_DATUM_1);
404	0	else
405	0	{
406		/* Build GEOGCS based on Ellipsoid (Table 3) */
407
408	0	if (hHeader.nEllipsoid > -1 &&
409	0	hHeader.nEllipsoid <
410	0	static_cast<GInt16>(CPL_ARRAYSIZE(EllipsoidTable)))
411	0	m_oSRS.SetGeogCS(
412	0	CPLSPrintf("BYN Ellipsoid(%d)", hHeader.nEllipsoid),
413	0	"Unspecified", EllipsoidTable[hHeader.nEllipsoid].pszName,
414	0	EllipsoidTable[hHeader.nEllipsoid].dfSemiMajor,
415	0	EllipsoidTable[hHeader.nEllipsoid].dfInvFlattening);
416	0	else
417	0	bNoGeogCS = true;
418	0	}
419
420		/* Build the VERT_CS based on VDatum */
421
422	601	OGRSpatialReference oSRSComp;
423	601	OGRSpatialReference oSRSVert;
424
425	601	int nVertCS = 0;
426
427	601	if (hHeader.nVDatum == 1)
428	436	nVertCS = BYN_VDATUM_1;
429	165	else if (hHeader.nVDatum == 2)
430	152	nVertCS = BYN_VDATUM_2;
431	13	else if (hHeader.nVDatum == 3)
432	5	nVertCS = BYN_VDATUM_3;
433	8	else
434	8	{
435		/* Return GEOGCS ( .err files ) */
436
437	8	if (bNoGeogCS)
438	0	return nullptr;
439
440	8	return &m_oSRS;
441	8	}
442
443	593	oSRSVert.importFromEPSG(nVertCS);
444
445		/* Create CPMPD_CS with GEOGCS and VERT_CS */
446
447	593	if (oSRSComp.SetCompoundCS(CPLSPrintf("BYN Datum(%d) & VDatum(%d)",
448	593	hHeader.nDatum, hHeader.nDatum),
449	593	&m_oSRS, &oSRSVert) == CE_None)
450	593	{
451		/* Return COMPD_CS with GEOGCS and VERT_CS */
452
453	593	m_oSRS = std::move(oSRSComp);
454	593	m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
455	593	return &m_oSRS;
456	593	}
457
458	0	return nullptr;
459	593	}
460
461		/----------------------------------------------------------------------/
462		/* buffer2header() */
463		/----------------------------------------------------------------------/
464
465		void BYNDataset::buffer2header(const GByte pabyBuf, BYNHeader pohHeader)
466
467	73.7k	{
468	73.7k	memcpy(&pohHeader->nSouth, pabyBuf, 4);
469	73.7k	memcpy(&pohHeader->nNorth, pabyBuf + 4, 4);
470	73.7k	memcpy(&pohHeader->nWest, pabyBuf + 8, 4);
471	73.7k	memcpy(&pohHeader->nEast, pabyBuf + 12, 4);
472	73.7k	memcpy(&pohHeader->nDLat, pabyBuf + 16, 2);
473	73.7k	memcpy(&pohHeader->nDLon, pabyBuf + 18, 2);
474	73.7k	memcpy(&pohHeader->nGlobal, pabyBuf + 20, 2);
475	73.7k	memcpy(&pohHeader->nType, pabyBuf + 22, 2);
476	73.7k	memcpy(&pohHeader->dfFactor, pabyBuf + 24, 8);
477	73.7k	memcpy(&pohHeader->nSizeOf, pabyBuf + 32, 2);
478	73.7k	memcpy(&pohHeader->nVDatum, pabyBuf + 34, 2);
479	73.7k	memcpy(&pohHeader->nDescrip, pabyBuf + 40, 2);
480	73.7k	memcpy(&pohHeader->nSubType, pabyBuf + 42, 2);
481	73.7k	memcpy(&pohHeader->nDatum, pabyBuf + 44, 2);
482	73.7k	memcpy(&pohHeader->nEllipsoid, pabyBuf + 46, 2);
483	73.7k	memcpy(&pohHeader->nByteOrder, pabyBuf + 48, 2);
484	73.7k	memcpy(&pohHeader->nScale, pabyBuf + 50, 2);
485	73.7k	memcpy(&pohHeader->dfWo, pabyBuf + 52, 8);
486	73.7k	memcpy(&pohHeader->dfGM, pabyBuf + 60, 8);
487	73.7k	memcpy(&pohHeader->nTideSys, pabyBuf + 68, 2);
488	73.7k	memcpy(&pohHeader->nRealiz, pabyBuf + 70, 2);
489	73.7k	memcpy(&pohHeader->dEpoch, pabyBuf + 72, 4);
490	73.7k	memcpy(&pohHeader->nPtType, pabyBuf + 76, 2);
491
492		#if defined(CPL_MSB)
493		CPL_LSBPTR32(&pohHeader->nSouth);
494		CPL_LSBPTR32(&pohHeader->nNorth);
495		CPL_LSBPTR32(&pohHeader->nWest);
496		CPL_LSBPTR32(&pohHeader->nEast);
497		CPL_LSBPTR16(&pohHeader->nDLat);
498		CPL_LSBPTR16(&pohHeader->nDLon);
499		CPL_LSBPTR16(&pohHeader->nGlobal);
500		CPL_LSBPTR16(&pohHeader->nType);
501		CPL_LSBPTR64(&pohHeader->dfFactor);
502		CPL_LSBPTR16(&pohHeader->nSizeOf);
503		CPL_LSBPTR16(&pohHeader->nVDatum);
504		CPL_LSBPTR16(&pohHeader->nDescrip);
505		CPL_LSBPTR16(&pohHeader->nSubType);
506		CPL_LSBPTR16(&pohHeader->nDatum);
507		CPL_LSBPTR16(&pohHeader->nEllipsoid);
508		CPL_LSBPTR16(&pohHeader->nByteOrder);
509		CPL_LSBPTR16(&pohHeader->nScale);
510		CPL_LSBPTR64(&pohHeader->dfWo);
511		CPL_LSBPTR64(&pohHeader->dfGM);
512		CPL_LSBPTR16(&pohHeader->nTideSys);
513		CPL_LSBPTR16(&pohHeader->nRealiz);
514		CPL_LSBPTR32(&pohHeader->dEpoch);
515		CPL_LSBPTR16(&pohHeader->nPtType);
516		#endif
517
518		#if DEBUG
519		CPLDebug("BYN", "South = %d", pohHeader->nSouth);
520		CPLDebug("BYN", "North = %d", pohHeader->nNorth);
521		CPLDebug("BYN", "West = %d", pohHeader->nWest);
522		CPLDebug("BYN", "East = %d", pohHeader->nEast);
523		CPLDebug("BYN", "DLat = %d", pohHeader->nDLat);
524		CPLDebug("BYN", "DLon = %d", pohHeader->nDLon);
525		CPLDebug("BYN", "DGlobal = %d", pohHeader->nGlobal);
526		CPLDebug("BYN", "DType = %d", pohHeader->nType);
527		CPLDebug("BYN", "Factor = %f", pohHeader->dfFactor);
528		CPLDebug("BYN", "SizeOf = %d", pohHeader->nSizeOf);
529		CPLDebug("BYN", "VDatum = %d", pohHeader->nVDatum);
530		CPLDebug("BYN", "Data = %d", pohHeader->nDescrip);
531		CPLDebug("BYN", "SubType = %d", pohHeader->nSubType);
532		CPLDebug("BYN", "Datum = %d", pohHeader->nDatum);
533		CPLDebug("BYN", "Ellipsoid = %d", pohHeader->nEllipsoid);
534		CPLDebug("BYN", "ByteOrder = %d", pohHeader->nByteOrder);
535		CPLDebug("BYN", "Scale = %d", pohHeader->nScale);
536		CPLDebug("BYN", "Wo = %f", pohHeader->dfWo);
537		CPLDebug("BYN", "GM = %f", pohHeader->dfGM);
538		CPLDebug("BYN", "TideSystem = %d", pohHeader->nTideSys);
539		CPLDebug("BYN", "RefRealzation = %d", pohHeader->nRealiz);
540		CPLDebug("BYN", "Epoch = %f", pohHeader->dEpoch);
541		CPLDebug("BYN", "PtType = %d", pohHeader->nPtType);
542		#endif
543	73.7k	}
544
545		/************************************************************************/
546		/* GDALRegister_BYN() */
547		/************************************************************************/
548
549		void GDALRegister_BYN()
550
551	24	{
552	24	if (GDALGetDriverByName("BYN") != nullptr)
553	0	return;
554
555	24	GDALDriver *poDriver = new GDALDriver();
556
557	24	poDriver->SetDescription("BYN");
558	24	poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
559	24	poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
560	24	"Natural Resources Canada's Geoid");
561	24	poDriver->SetMetadataItem(GDAL_DMD_EXTENSIONS, "byn err");
562	24	poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
563	24	poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/byn.html");
564
565	24	poDriver->pfnOpen = BYNDataset::Open;
566	24	poDriver->pfnIdentify = BYNDataset::Identify;
567
568	24	GetGDALDriverManager()->RegisterDriver(poDriver);
569	24	}