/src/gdal/ogr/ogrsf_frmts/gpkg/ogrgeopackageutility.cpp

Source
/******************************************************************************
 *
 * Project:  GeoPackage Translator
 * Purpose:  Utility functions for OGR GeoPackage driver.
 * Author:   Paul Ramsey, pramsey@boundlessgeo.com
 *
 ******************************************************************************
 * Copyright (c) 2013, Paul Ramsey <pramsey@boundlessgeo.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "ogrgeopackageutility.h"
#include "ogr_p.h"
#include "ogr_wkb.h"
#include "sqlite/ogrsqlitebase.h"
#include <limits>

/* Requirement 20: A GeoPackage SHALL store feature table geometries */
/* with the basic simple feature geometry types (Geometry, Point, */
/* LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon, */
/* GeomCollection) */
/* http://opengis.github.io/geopackage/#geometry_types */
OGRwkbGeometryType GPkgGeometryTypeToWKB(const char *pszGpkgType, bool bHasZ,
                                         bool bHasM)
{
    OGRwkbGeometryType oType;

    if (EQUAL("Geometry", pszGpkgType))
        oType = wkbUnknown;
    /* The 1.0 spec is not completely clear on what should be used... */
    else if (EQUAL("GeomCollection", pszGpkgType) ||
             EQUAL("GeometryCollection", pszGpkgType))
        oType = wkbGeometryCollection;
    else
    {
        oType = OGRFromOGCGeomType(pszGpkgType);
        if (oType == wkbUnknown)
            oType = wkbNone;
    }

    if ((oType != wkbNone) && bHasZ)
    {
        oType = wkbSetZ(oType);
    }
    if ((oType != wkbNone) && bHasM)
    {
        oType = wkbSetM(oType);
    }

    return oType;
}

/* Requirement 5: The columns of tables in a GeoPackage SHALL only be */
/* declared using one of the data types specified in table GeoPackage */
/* Data Types. */
/* http://opengis.github.io/geopackage/#table_column_data_types */
// return a OGRFieldType value or OFTMaxType + 1
int GPkgFieldToOGR(const char *pszGpkgType, OGRFieldSubType &eSubType,
                   int &nMaxWidth)
{
    eSubType = OFSTNone;
    nMaxWidth = 0;

    /* Integer types */
    if (STRNCASECMP("INT", pszGpkgType, 3) == 0)
    {
        if (!EQUAL("INT", pszGpkgType) && !EQUAL("INTEGER", pszGpkgType))
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                     "Field format '%s' not supported. "
                     "Interpreted as INT",
                     pszGpkgType);
        }
        return OFTInteger64;
    }
    else if (EQUAL("MEDIUMINT", pszGpkgType))
        return OFTInteger;
    else if (EQUAL("SMALLINT", pszGpkgType))
    {
        eSubType = OFSTInt16;
        return OFTInteger;
    }
    else if (EQUAL("TINYINT", pszGpkgType))
        return OFTInteger;  // [-128, 127]
    else if (EQUAL("BOOLEAN", pszGpkgType))
    {
        eSubType = OFSTBoolean;
        return OFTInteger;
    }

    /* Real types */
    else if (EQUAL("FLOAT", pszGpkgType))
    {
        eSubType = OFSTFloat32;
        return OFTReal;
    }
    else if (EQUAL("DOUBLE", pszGpkgType))
        return OFTReal;
    else if (EQUAL("REAL", pszGpkgType))
        return OFTReal;

    // Only used normally in gpkg_data_column_constraints table, and we
    // need this only is reading it through ExecuteSQL()
    else if (EQUAL("NUMERIC", pszGpkgType))
        return OFTReal;

    /* String/binary types */
    else if (STRNCASECMP("TEXT", pszGpkgType, 4) == 0)
    {
        if (pszGpkgType[4] == '(')
            nMaxWidth = atoi(pszGpkgType + 5);
        else if (pszGpkgType[4] != '\0')
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                     "Field format '%s' not supported. "
                     "Interpreted as TEXT",
                     pszGpkgType);
        }
        return OFTString;
    }

    else if (STRNCASECMP("BLOB", pszGpkgType, 4) == 0)
    {
        if (pszGpkgType[4] != '(' && pszGpkgType[4] != '\0')
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                     "Field format '%s' not supported. "
                     "Interpreted as BLOB",
                     pszGpkgType);
        }
        return OFTBinary;
    }

    /* Date types */
    else if (EQUAL("DATE", pszGpkgType))
        return OFTDate;
    else if (EQUAL("DATETIME", pszGpkgType))
        return OFTDateTime;

    /* Illegal! */
    else
    {
        if (GPkgGeometryTypeToWKB(pszGpkgType, false, false) == wkbNone)
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                     "Field format '%s' not supported", pszGpkgType);
        }
        return OFTMaxType + 1;
    }
}

/* Requirement 5: The columns of tables in a GeoPackage SHALL only be */
/* declared using one of the data types specified in table GeoPackage */
/* Data Types. */
/* http://opengis.github.io/geopackage/#table_column_data_types */
const char *GPkgFieldFromOGR(OGRFieldType eType, OGRFieldSubType eSubType,
                             int nMaxWidth)
{
    switch (eType)
    {
        case OFTInteger:
        {
            if (eSubType == OFSTBoolean)
                return "BOOLEAN";
            else if (eSubType == OFSTInt16)
                return "SMALLINT";
            else
                return "MEDIUMINT";
        }
        case OFTInteger64:
            return "INTEGER";
        case OFTReal:
        {
            if (eSubType == OFSTFloat32)
                return "FLOAT";
            else
                return "REAL";
        }
        case OFTString:
        {
            if (nMaxWidth > 0)
                return CPLSPrintf("TEXT(%d)", nMaxWidth);
            else
                return "TEXT";
        }
        case OFTBinary:
            return "BLOB";
        case OFTDate:
            return "DATE";
        case OFTDateTime:
            return "DATETIME";
        default:
            return "TEXT";
    }
}

/* Requirement 19: A GeoPackage SHALL store feature table geometries
 *  with or without optional elevation (Z) and/or measure (M) values in SQL
 *  BLOBs using the Standard GeoPackageBinary format specified in table
 * GeoPackage SQL Geometry Binary Format and clause Geometry Encoding.
 *
 *  http://opengis.github.io/geopackage/#gpb_format
 *
 *   GeoPackageBinaryHeader {
 *     byte[2] magic = 0x4750;
 *     byte version;
 *     byte flags;
 *     int32 srs_id;
 *     double[] envelope;
 *    }
 *
 *   StandardGeoPackageBinary {
 *     GeoPackageBinaryHeader header;
 *     WKBGeometry geometry;
 *   }
 *
 *  Flags byte contents:
 *  Bit 7: Reserved for future
 *  Bit 6: Reserved for future
 *  Bit 5: Using Extended GPKG Binary?
 *  Bit 4: Geometry is Empty?
 *  Bit 3,2,1: Envelope contents (0 none, 1=X/Y, 2=X/Y/Z, 3=X/Y/M, 4=X/Y/Z/M)
 *  Bit 0: Byte order of header (0=big/XDR, 1=little/NDR)
 *
 */

GByte *GPkgGeometryFromOGR(const OGRGeometry *poGeometry, int iSrsId,
                           const OGRGeomCoordinateBinaryPrecision *psPrecision,
                           size_t *pnWkbLen)
{
    CPLAssert(poGeometry != nullptr);

    GByte byFlags = 0;
    GByte byEnv = 1;
    OGRwkbExportOptions wkbExportOptions;
    if (psPrecision)
        wkbExportOptions.sPrecision = *psPrecision;
    wkbExportOptions.eByteOrder = static_cast<OGRwkbByteOrder>(CPL_IS_LSB);
    OGRErr err;
    bool bPoint = (wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
    bool bEmpty = poGeometry->IsEmpty();
    /* We voluntarily use getCoordinateDimension() so as to get only 2 for
     * XY/XYM */
    /* and 3 for XYZ/XYZM as we currently don't write envelopes with M extent.
     */
    int iDims = poGeometry->getCoordinateDimension();

    /* Header has 8 bytes for sure, and optional extra space for bounds */
    size_t nHeaderLen = 2 + 1 + 1 + 4;
    if (!bPoint && !bEmpty)
    {
        nHeaderLen += 8 * 2 * iDims;
    }

    /* Total BLOB size is header + WKB size */
    size_t nWkbLen = nHeaderLen + poGeometry->WkbSize();
    if (nWkbLen > static_cast<size_t>(std::numeric_limits<int>::max()))
    {
        CPLError(CE_Failure, CPLE_NotSupported, "too big geometry blob");
        return nullptr;
    }
    GByte *pabyWkb = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nWkbLen));
    if (!pabyWkb)
        return nullptr;
    if (pnWkbLen)
        *pnWkbLen = nWkbLen;

    /* Header Magic */
    pabyWkb[0] = 0x47;
    pabyWkb[1] = 0x50;

    /* GPKG BLOB Version */
    pabyWkb[2] = 0;

    /* Extended? No. */

    /* Envelope dimensionality? */

    /* Don't write envelope for point type */
    if (bPoint)
        byEnv = 0;
    else
        /* 3D envelope for 3D data */
        if (iDims == 3)
            byEnv = 2;
        /* 2D envelope otherwise */
        else
            byEnv = 1;

    /* Empty? No envelope then. */
    if (bEmpty)
    {
        byEnv = 0;
        /* Set empty flag */
        byFlags |= (1 << 4);
    }

    /* Set envelope flags */
    byFlags |= (byEnv << 1);

    /* Byte order of header? */
    /* Use native endianness */
    byFlags |= wkbExportOptions.eByteOrder;

    /* Write flags byte */
    pabyWkb[3] = byFlags;

    /* Write srs_id */
    memcpy(pabyWkb + 4, &iSrsId, 4);

    /* Write envelope */
    if (!bEmpty && !bPoint)
    {
        double *padPtr = reinterpret_cast<double *>(pabyWkb + 8);
        if (iDims == 3)
        {
            OGREnvelope3D oEnv3d;
            poGeometry->getEnvelope(&oEnv3d);
            padPtr[0] = oEnv3d.MinX;
            padPtr[1] = oEnv3d.MaxX;
            padPtr[2] = oEnv3d.MinY;
            padPtr[3] = oEnv3d.MaxY;
            padPtr[4] = oEnv3d.MinZ;
            padPtr[5] = oEnv3d.MaxZ;
        }
        else
        {
            OGREnvelope oEnv;
            poGeometry->getEnvelope(&oEnv);
            padPtr[0] = oEnv.MinX;
            padPtr[1] = oEnv.MaxX;
            padPtr[2] = oEnv.MinY;
            padPtr[3] = oEnv.MaxY;
        }
    }

    GByte *pabyPtr = pabyWkb + nHeaderLen;

    /* Use the wkbVariantIso for ISO SQL/MM output (differs for 3d geometry) */
    wkbExportOptions.eWkbVariant = wkbVariantIso;
    err = poGeometry->exportToWkb(pabyPtr, &wkbExportOptions);
    if (err != OGRERR_NONE)
    {
        CPLFree(pabyWkb);
        return nullptr;
    }

    return pabyWkb;
}

OGRErr GPkgHeaderFromWKB(const GByte *pabyGpkg, size_t nGpkgLen,
                         GPkgHeader *poHeader)
{
    CPLAssert(pabyGpkg != nullptr);
    CPLAssert(poHeader != nullptr);

    /* Magic (match required) */
    if (nGpkgLen < 8 || pabyGpkg[0] != 0x47 || pabyGpkg[1] != 0x50 ||
        pabyGpkg[2] != 0) /* Version (only 0 supported at this time)*/
    {
        memset(poHeader, 0, sizeof(*poHeader));
        return OGRERR_FAILURE;
    }

    /* Flags */
    GByte byFlags = pabyGpkg[3];
    poHeader->bEmpty = (byFlags & (0x01 << 4)) >> 4;
    poHeader->bExtended = (byFlags & (0x01 << 5)) >> 5;
    poHeader->eByteOrder = static_cast<OGRwkbByteOrder>(byFlags & 0x01);
    poHeader->bExtentHasXY = false;
    poHeader->bExtentHasZ = false;
#ifdef notdef
    poHeader->bExtentHasM = false;
#endif
    bool bSwap = OGR_SWAP(poHeader->eByteOrder);

    /* Envelope */
    int iEnvelope = (byFlags & (0x07 << 1)) >> 1;
    int nEnvelopeDim = 0;
    if (iEnvelope)
    {
        poHeader->bExtentHasXY = true;
        if (iEnvelope == 1)
        {
            nEnvelopeDim = 2; /* 2D envelope */
        }
        else if (iEnvelope == 2)
        {
            poHeader->bExtentHasZ = true;
            nEnvelopeDim = 3; /* 2D+Z envelope */
        }
        else if (iEnvelope == 3)
        {
#ifdef notdef
            poHeader->bExtentHasM = true;
#endif
            nEnvelopeDim = 3; /* 2D+M envelope */
        }
        else if (iEnvelope == 4)
        {
            poHeader->bExtentHasZ = true;
#ifdef notdef
            poHeader->bExtentHasM = true;
#endif
            nEnvelopeDim = 4; /* 2D+ZM envelope */
        }
        else
        {
            return OGRERR_FAILURE;
        }
    }

    /* SrsId */
    int iSrsId = 0;
    memcpy(&iSrsId, pabyGpkg + 4, 4);
    if (bSwap)
    {
        iSrsId = CPL_SWAP32(iSrsId);
    }
    poHeader->iSrsId = iSrsId;

    if (nGpkgLen < static_cast<size_t>(8 + 8 * 2 * nEnvelopeDim))
    {
        // Not enough bytes
        return OGRERR_FAILURE;
    }

    /* Envelope */
    const double *padPtr = reinterpret_cast<const double *>(pabyGpkg + 8);
    if (poHeader->bExtentHasXY)
    {
        poHeader->MinX = padPtr[0];
        poHeader->MaxX = padPtr[1];
        poHeader->MinY = padPtr[2];
        poHeader->MaxY = padPtr[3];
        if (bSwap)
        {
            CPL_SWAPDOUBLE(&(poHeader->MinX));
            CPL_SWAPDOUBLE(&(poHeader->MaxX));
            CPL_SWAPDOUBLE(&(poHeader->MinY));
            CPL_SWAPDOUBLE(&(poHeader->MaxY));
        }
    }
    if (poHeader->bExtentHasZ)
    {
        poHeader->MinZ = padPtr[4];
        poHeader->MaxZ = padPtr[5];
        if (bSwap)
        {
            CPL_SWAPDOUBLE(&(poHeader->MinZ));
            CPL_SWAPDOUBLE(&(poHeader->MaxZ));
        }
    }
#ifdef notdef
    if (poHeader->bExtentHasM)
    {
        poHeader->MinM = padPtr[(poHeader->bExtentHasZ) ? 6 : 4];
        poHeader->MaxM = padPtr[(poHeader->bExtentHasZ) ? 7 : 5];
        if (bSwap)
        {
            CPL_SWAPDOUBLE(&(poHeader->MinM));
            CPL_SWAPDOUBLE(&(poHeader->MaxM));
        }
    }
#endif

    /* Header size in byte stream */
    poHeader->nHeaderLen = 8 + 8 * 2 * nEnvelopeDim;

#ifdef DEBUG_VERBOSE
    std::string s;
    for (size_t i = poHeader->nHeaderLen; i < nGpkgLen; ++i)
    {
        s += CPLSPrintf("%02X ", pabyGpkg[i]);
    }
    CPLDebug("GPKG", "Bytes after GPKG header: %s", s.c_str());
#endif

    // Workaround for a Spatialite bug. The CastToXYZ() function, when
    // called on an empty geometry, and EnableGpkgMode() is enabled,
    // returns an empty geometry, not tagged as such, but with an invalid
    // bounding box.
    // Cf https://github.com/OSGeo/gdal/issues/13557
    if (!poHeader->bEmpty && poHeader->bExtentHasXY &&
        poHeader->nHeaderLen == 40 &&
        poHeader->MinX == std::numeric_limits<double>::max() &&
        poHeader->MaxX == -std::numeric_limits<double>::max() &&
        poHeader->MinY == std::numeric_limits<double>::max() &&
        poHeader->MaxY == -std::numeric_limits<double>::max() &&
        // CastToXYZ(POLYGON EMPTY) returns just the GPKG header for some reason
        (nGpkgLen == 40 ||
         // CastToXYZ(LINESTRING EMPTY) returns a LINESTRING Z EMPTY
         (nGpkgLen == 49 &&
          memcmp(pabyGpkg + 40, "\x01\xEA\x03\x00\x00\x00\x00\x00\x00", 9) ==
              0) ||
         // CastToXYZ(POINT EMPTY) returns a Point Z (NaN NaN 0)
         // CastToXYZ(POINT Z EMPTY) returns a Point Z (NaN NaN NaN)
         (nGpkgLen == 69 && memcmp(pabyGpkg + 40,
                                   "\x01"
                                   "\xE9\x03\x00\x00"
                                   "\x00\x00\x00\x00\x00\x00\xF8\x7F"
                                   "\x00\x00\x00\x00\x00\x00\xF8\x7F",
                                   21) == 0)))
    {
        CPLDebugOnce("GPKG",
                     "Work arounding a Spatialite bug with empty geometries");
        poHeader->bEmpty = true;
        poHeader->bExtentHasXY = false;
    }

    return OGRERR_NONE;
}

bool GPkgUpdateHeader(GByte *pabyGpkg, size_t nGpkgLen, int nSrsId, double MinX,
                      double MaxX, double MinY, double MaxY, double MinZ,
                      double MaxZ)
{
    CPLAssert(nGpkgLen >= 8);

    /* Flags */
    const GByte byFlags = pabyGpkg[3];
    const auto eByteOrder = static_cast<OGRwkbByteOrder>(byFlags & 0x01);
    const bool bSwap = OGR_SWAP(eByteOrder);

    /* SrsId */
    if (bSwap)
    {
        nSrsId = CPL_SWAP32(nSrsId);
    }
    memcpy(pabyGpkg + 4, &nSrsId, 4);

    /* Envelope */
    const int iEnvelope = (byFlags & (0x07 << 1)) >> 1;
    int nEnvelopeDim = 0;
    if (iEnvelope)
    {
        if (iEnvelope == 1)
        {
            nEnvelopeDim = 2; /* 2D envelope */
        }
        else if (iEnvelope == 2)
        {
            nEnvelopeDim = 3; /* 2D+Z envelope */
        }
        else if (iEnvelope == 3)
        {
            nEnvelopeDim = 3; /* 2D+M envelope */
        }
        else if (iEnvelope == 4)
        {
            nEnvelopeDim = 4; /* 2D+ZM envelope */
        }
        else
        {
            return false;
        }
    }
    else
    {
        return true;
    }

    if (nGpkgLen < static_cast<size_t>(8 + 8 * 2 * nEnvelopeDim))
    {
        // Not enough bytes
        return false;
    }

    /* Envelope */
    if (bSwap)
    {
        CPL_SWAPDOUBLE(&(MinX));
        CPL_SWAPDOUBLE(&(MaxX));
        CPL_SWAPDOUBLE(&(MinY));
        CPL_SWAPDOUBLE(&(MaxY));
        CPL_SWAPDOUBLE(&(MinZ));
        CPL_SWAPDOUBLE(&(MaxZ));
    }

    double *padPtr = reinterpret_cast<double *>(pabyGpkg + 8);
    memcpy(&padPtr[0], &MinX, sizeof(double));
    memcpy(&padPtr[1], &MaxX, sizeof(double));
    memcpy(&padPtr[2], &MinY, sizeof(double));
    memcpy(&padPtr[3], &MaxY, sizeof(double));

    if (iEnvelope == 2 || iEnvelope == 4)
    {
        memcpy(&padPtr[4], &MinZ, sizeof(double));
        memcpy(&padPtr[5], &MaxZ, sizeof(double));
    }

    return true;
}

OGRGeometry *GPkgGeometryToOGR(const GByte *pabyGpkg, size_t nGpkgLen,
                               OGRSpatialReference *poSrs)
{
    CPLAssert(pabyGpkg != nullptr);

    GPkgHeader oHeader;

    /* Read header */
    OGRErr err = GPkgHeaderFromWKB(pabyGpkg, nGpkgLen, &oHeader);
    if (err != OGRERR_NONE)
        return nullptr;

    /* WKB pointer */
    const GByte *pabyWkb = pabyGpkg + oHeader.nHeaderLen;
    size_t nWkbLen = nGpkgLen - oHeader.nHeaderLen;

    /* Parse WKB */
    OGRGeometry *poGeom = nullptr;
    err = OGRGeometryFactory::createFromWkb(pabyWkb, poSrs, &poGeom,
                                            static_cast<int>(nWkbLen));
    if (err != OGRERR_NONE)
        return nullptr;

    return poGeom;
}

/************************************************************************/
/*                       OGRGeoPackageGetHeader()                       */
/************************************************************************/

bool OGRGeoPackageGetHeader(sqlite3_context * /*pContext*/, int /*argc*/,
                            sqlite3_value **argv, GPkgHeader *psHeader,
                            bool bNeedExtent, bool bNeedExtent3D, int iGeomIdx)
{

    // Extent3D implies extent
    const bool bNeedAnyExtent{bNeedExtent || bNeedExtent3D};

    if (sqlite3_value_type(argv[iGeomIdx]) != SQLITE_BLOB)
    {
        memset(psHeader, 0, sizeof(*psHeader));
        return false;
    }
    const int nBLOBLen = sqlite3_value_bytes(argv[iGeomIdx]);
    const GByte *pabyBLOB =
        reinterpret_cast<const GByte *>(sqlite3_value_blob(argv[iGeomIdx]));

    if (nBLOBLen < 8)
    {
        memset(psHeader, 0, sizeof(*psHeader));
        return false;
    }
    else if (GPkgHeaderFromWKB(pabyBLOB, nBLOBLen, psHeader) != OGRERR_NONE)
    {
        bool bEmpty = false;
        memset(psHeader, 0, sizeof(*psHeader));
        if (OGRSQLiteGetSpatialiteGeometryHeader(
                pabyBLOB, nBLOBLen, &(psHeader->iSrsId), nullptr, &bEmpty,
                &(psHeader->MinX), &(psHeader->MinY), &(psHeader->MaxX),
                &(psHeader->MaxY)) == OGRERR_NONE)
        {
            psHeader->bEmpty = bEmpty;
            psHeader->bExtentHasXY = !bEmpty;
            if (!bNeedExtent3D && !(bEmpty && bNeedAnyExtent))
                return true;
        }

        return false;
    }

    if (psHeader->bEmpty && bNeedAnyExtent)
    {
        return false;
    }
    else if (!psHeader->bExtentHasXY && bNeedExtent && !bNeedExtent3D)
    {
        OGREnvelope sEnvelope;
        if (OGRWKBGetBoundingBox(pabyBLOB + psHeader->nHeaderLen,
                                 static_cast<size_t>(nBLOBLen) -
                                     psHeader->nHeaderLen,
                                 sEnvelope))
        {
            psHeader->MinX = sEnvelope.MinX;
            psHeader->MaxX = sEnvelope.MaxX;
            psHeader->MinY = sEnvelope.MinY;
            psHeader->MaxY = sEnvelope.MaxY;
            return true;
        }
        return false;
    }
    else if (!psHeader->bExtentHasZ && bNeedExtent3D)
    {
        OGREnvelope3D sEnvelope3D;
        if (OGRWKBGetBoundingBox(pabyBLOB + psHeader->nHeaderLen,
                                 static_cast<size_t>(nBLOBLen) -
                                     psHeader->nHeaderLen,
                                 sEnvelope3D))
        {
            psHeader->MinX = sEnvelope3D.MinX;
            psHeader->MaxX = sEnvelope3D.MaxX;
            psHeader->MinY = sEnvelope3D.MinY;
            psHeader->MaxY = sEnvelope3D.MaxY;
            psHeader->MinZ = sEnvelope3D.MinZ;
            psHeader->MaxZ = sEnvelope3D.MaxZ;
            return true;
        }
        return false;
    }
    return true;
}

Coverage Report

Created: 2026-05-16 08:20

Line	Count	Source
1		/******************************************************************************
2		*
3		* Project: GeoPackage Translator
4		* Purpose: Utility functions for OGR GeoPackage driver.
5		* Author: Paul Ramsey, pramsey@boundlessgeo.com
6		*
7		******************************************************************************
8		* Copyright (c) 2013, Paul Ramsey <pramsey@boundlessgeo.com>
9		*
10		* SPDX-License-Identifier: MIT
11		****************************************************************************/
12
13		#include "ogrgeopackageutility.h"
14		#include "ogr_p.h"
15		#include "ogr_wkb.h"
16		#include "sqlite/ogrsqlitebase.h"
17		#include <limits>
18
19		/* Requirement 20: A GeoPackage SHALL store feature table geometries */
20		/* with the basic simple feature geometry types (Geometry, Point, */
21		/* LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon, */
22		/* GeomCollection) */
23		/* http://opengis.github.io/geopackage/#geometry_types */
24		OGRwkbGeometryType GPkgGeometryTypeToWKB(const char *pszGpkgType, bool bHasZ,
25		bool bHasM)
26	5.58k	{
27	5.58k	OGRwkbGeometryType oType;
28
29	5.58k	if (EQUAL("Geometry", pszGpkgType))
30	366	oType = wkbUnknown;
31		/* The 1.0 spec is not completely clear on what should be used... */
32	5.22k	else if (EQUAL("GeomCollection", pszGpkgType) \|\|
33	5.22k	EQUAL("GeometryCollection", pszGpkgType))
34	0	oType = wkbGeometryCollection;
35	5.22k	else
36	5.22k	{
37	5.22k	oType = OGRFromOGCGeomType(pszGpkgType);
38	5.22k	if (oType == wkbUnknown)
39	3.47k	oType = wkbNone;
40	5.22k	}
41
42	5.58k	if ((oType != wkbNone) && bHasZ)
43	313	{
44	313	oType = wkbSetZ(oType);
45	313	}
46	5.58k	if ((oType != wkbNone) && bHasM)
47	528	{
48	528	oType = wkbSetM(oType);
49	528	}
50
51	5.58k	return oType;
52	5.58k	}
53
54		/* Requirement 5: The columns of tables in a GeoPackage SHALL only be */
55		/* declared using one of the data types specified in table GeoPackage */
56		/* Data Types. */
57		/* http://opengis.github.io/geopackage/#table_column_data_types */
58		// return a OGRFieldType value or OFTMaxType + 1
59		int GPkgFieldToOGR(const char *pszGpkgType, OGRFieldSubType &eSubType,
60		int &nMaxWidth)
61	5.04k	{
62	5.04k	eSubType = OFSTNone;
63	5.04k	nMaxWidth = 0;
64
65		/* Integer types */
66	5.04k	if (STRNCASECMP("INT", pszGpkgType, 3) == 0)
67	939	{
68	939	if (!EQUAL("INT", pszGpkgType) && !EQUAL("INTEGER", pszGpkgType))
69	461	{
70	461	CPLError(CE_Warning, CPLE_AppDefined,
71	461	"Field format '%s' not supported. "
72	461	"Interpreted as INT",
73	461	pszGpkgType);
74	461	}
75	939	return OFTInteger64;
76	939	}
77	4.11k	else if (EQUAL("MEDIUMINT", pszGpkgType))
78	8	return OFTInteger;
79	4.10k	else if (EQUAL("SMALLINT", pszGpkgType))
80	0	{
81	0	eSubType = OFSTInt16;
82	0	return OFTInteger;
83	0	}
84	4.10k	else if (EQUAL("TINYINT", pszGpkgType))
85	0	return OFTInteger; // [-128, 127]
86	4.10k	else if (EQUAL("BOOLEAN", pszGpkgType))
87	0	{
88	0	eSubType = OFSTBoolean;
89	0	return OFTInteger;
90	0	}
91
92		/* Real types */
93	4.10k	else if (EQUAL("FLOAT", pszGpkgType))
94	24	{
95	24	eSubType = OFSTFloat32;
96	24	return OFTReal;
97	24	}
98	4.07k	else if (EQUAL("DOUBLE", pszGpkgType))
99	0	return OFTReal;
100	4.07k	else if (EQUAL("REAL", pszGpkgType))
101	316	return OFTReal;
102
103		// Only used normally in gpkg_data_column_constraints table, and we
104		// need this only is reading it through ExecuteSQL()
105	3.76k	else if (EQUAL("NUMERIC", pszGpkgType))
106	0	return OFTReal;
107
108		/* String/binary types */
109	3.76k	else if (STRNCASECMP("TEXT", pszGpkgType, 4) == 0)
110	117	{
111	117	if (pszGpkgType[4] == '(')
112	1	nMaxWidth = atoi(pszGpkgType + 5);
113	116	else if (pszGpkgType[4] != '\0')
114	0	{
115	0	CPLError(CE_Warning, CPLE_AppDefined,
116	0	"Field format '%s' not supported. "
117	0	"Interpreted as TEXT",
118	0	pszGpkgType);
119	0	}
120	117	return OFTString;
121	117	}
122
123	3.64k	else if (STRNCASECMP("BLOB", pszGpkgType, 4) == 0)
124	69	{
125	69	if (pszGpkgType[4] != '(' && pszGpkgType[4] != '\0')
126	11	{
127	11	CPLError(CE_Warning, CPLE_AppDefined,
128	11	"Field format '%s' not supported. "
129	11	"Interpreted as BLOB",
130	11	pszGpkgType);
131	11	}
132	69	return OFTBinary;
133	69	}
134
135		/* Date types */
136	3.57k	else if (EQUAL("DATE", pszGpkgType))
137	0	return OFTDate;
138	3.57k	else if (EQUAL("DATETIME", pszGpkgType))
139	5	return OFTDateTime;
140
141		/* Illegal! */
142	3.57k	else
143	3.57k	{
144	3.57k	if (GPkgGeometryTypeToWKB(pszGpkgType, false, false) == wkbNone)
145	2.98k	{
146	2.98k	CPLError(CE_Warning, CPLE_AppDefined,
147	2.98k	"Field format '%s' not supported", pszGpkgType);
148	2.98k	}
149	3.57k	return OFTMaxType + 1;
150	3.57k	}
151	5.04k	}
152
153		/* Requirement 5: The columns of tables in a GeoPackage SHALL only be */
154		/* declared using one of the data types specified in table GeoPackage */
155		/* Data Types. */
156		/* http://opengis.github.io/geopackage/#table_column_data_types */
157		const char *GPkgFieldFromOGR(OGRFieldType eType, OGRFieldSubType eSubType,
158		int nMaxWidth)
159	28.0k	{
160	28.0k	switch (eType)
161	28.0k	{
162	160	case OFTInteger:
163	160	{
164	160	if (eSubType == OFSTBoolean)
165	37	return "BOOLEAN";
166	123	else if (eSubType == OFSTInt16)
167	0	return "SMALLINT";
168	123	else
169	123	return "MEDIUMINT";
170	160	}
171	32	case OFTInteger64:
172	32	return "INTEGER";
173	82	case OFTReal:
174	82	{
175	82	if (eSubType == OFSTFloat32)
176	0	return "FLOAT";
177	82	else
178	82	return "REAL";
179	82	}
180	27.7k	case OFTString:
181	27.7k	{
182	27.7k	if (nMaxWidth > 0)
183	169	return CPLSPrintf("TEXT(%d)", nMaxWidth);
184	27.5k	else
185	27.5k	return "TEXT";
186	27.7k	}
187	0	case OFTBinary:
188	0	return "BLOB";
189	0	case OFTDate:
190	0	return "DATE";
191	0	case OFTDateTime:
192	0	return "DATETIME";
193	0	default:
194	0	return "TEXT";
195	28.0k	}
196	28.0k	}
197
198		/* Requirement 19: A GeoPackage SHALL store feature table geometries
199		* with or without optional elevation (Z) and/or measure (M) values in SQL
200		* BLOBs using the Standard GeoPackageBinary format specified in table
201		* GeoPackage SQL Geometry Binary Format and clause Geometry Encoding.
202		*
203		* http://opengis.github.io/geopackage/#gpb_format
204		*
205		* GeoPackageBinaryHeader {
206		* byte[2] magic = 0x4750;
207		* byte version;
208		* byte flags;
209		* int32 srs_id;
210		* double[] envelope;
211		* }
212		*
213		* StandardGeoPackageBinary {
214		* GeoPackageBinaryHeader header;
215		* WKBGeometry geometry;
216		* }
217		*
218		* Flags byte contents:
219		* Bit 7: Reserved for future
220		* Bit 6: Reserved for future
221		* Bit 5: Using Extended GPKG Binary?
222		* Bit 4: Geometry is Empty?
223		* Bit 3,2,1: Envelope contents (0 none, 1=X/Y, 2=X/Y/Z, 3=X/Y/M, 4=X/Y/Z/M)
224		* Bit 0: Byte order of header (0=big/XDR, 1=little/NDR)
225		*
226		*/
227
228		GByte GPkgGeometryFromOGR(const OGRGeometry poGeometry, int iSrsId,
229		const OGRGeomCoordinateBinaryPrecision *psPrecision,
230		size_t *pnWkbLen)
231	28.5k	{
232	28.5k	CPLAssert(poGeometry != nullptr);
233
234	28.5k	GByte byFlags = 0;
235	28.5k	GByte byEnv = 1;
236	28.5k	OGRwkbExportOptions wkbExportOptions;
237	28.5k	if (psPrecision)
238	28.5k	wkbExportOptions.sPrecision = *psPrecision;
239	28.5k	wkbExportOptions.eByteOrder = static_cast<OGRwkbByteOrder>(CPL_IS_LSB);
240	28.5k	OGRErr err;
241	28.5k	bool bPoint = (wkbFlatten(poGeometry->getGeometryType()) == wkbPoint);
242	28.5k	bool bEmpty = poGeometry->IsEmpty();
243		/* We voluntarily use getCoordinateDimension() so as to get only 2 for
244		* XY/XYM */
245		/* and 3 for XYZ/XYZM as we currently don't write envelopes with M extent.
246		*/
247	28.5k	int iDims = poGeometry->getCoordinateDimension();
248
249		/* Header has 8 bytes for sure, and optional extra space for bounds */
250	28.5k	size_t nHeaderLen = 2 + 1 + 1 + 4;
251	28.5k	if (!bPoint && !bEmpty)
252	3.99k	{
253	3.99k	nHeaderLen += 8 * 2 * iDims;
254	3.99k	}
255
256		/* Total BLOB size is header + WKB size */
257	28.5k	size_t nWkbLen = nHeaderLen + poGeometry->WkbSize();
258	28.5k	if (nWkbLen > static_cast<size_t>(std::numeric_limits<int>::max()))
259	0	{
260	0	CPLError(CE_Failure, CPLE_NotSupported, "too big geometry blob");
261	0	return nullptr;
262	0	}
263	28.5k	GByte pabyWkb = static_cast<GByte >(VSI_MALLOC_VERBOSE(nWkbLen));
264	28.5k	if (!pabyWkb)
265	0	return nullptr;
266	28.5k	if (pnWkbLen)
267	28.5k	*pnWkbLen = nWkbLen;
268
269		/* Header Magic */
270	28.5k	pabyWkb[0] = 0x47;
271	28.5k	pabyWkb[1] = 0x50;
272
273		/* GPKG BLOB Version */
274	28.5k	pabyWkb[2] = 0;
275
276		/* Extended? No. */
277
278		/* Envelope dimensionality? */
279
280		/* Don't write envelope for point type */
281	28.5k	if (bPoint)
282	11.6k	byEnv = 0;
283	16.9k	else
284		/* 3D envelope for 3D data */
285	16.9k	if (iDims == 3)
286	1.61k	byEnv = 2;
287		/* 2D envelope otherwise */
288	15.2k	else
289	15.2k	byEnv = 1;
290
291		/* Empty? No envelope then. */
292	28.5k	if (bEmpty)
293	13.1k	{
294	13.1k	byEnv = 0;
295		/* Set empty flag */
296	13.1k	byFlags \|= (1 << 4);
297	13.1k	}
298
299		/* Set envelope flags */
300	28.5k	byFlags \|= (byEnv << 1);
301
302		/* Byte order of header? */
303		/* Use native endianness */
304	28.5k	byFlags \|= wkbExportOptions.eByteOrder;
305
306		/* Write flags byte */
307	28.5k	pabyWkb[3] = byFlags;
308
309		/* Write srs_id */
310	28.5k	memcpy(pabyWkb + 4, &iSrsId, 4);
311
312		/* Write envelope */
313	28.5k	if (!bEmpty && !bPoint)
314	3.99k	{
315	3.99k	double padPtr = reinterpret_cast<double >(pabyWkb + 8);
316	3.99k	if (iDims == 3)
317	553	{
318	553	OGREnvelope3D oEnv3d;
319	553	poGeometry->getEnvelope(&oEnv3d);
320	553	padPtr[0] = oEnv3d.MinX;
321	553	padPtr[1] = oEnv3d.MaxX;
322	553	padPtr[2] = oEnv3d.MinY;
323	553	padPtr[3] = oEnv3d.MaxY;
324	553	padPtr[4] = oEnv3d.MinZ;
325	553	padPtr[5] = oEnv3d.MaxZ;
326	553	}
327	3.44k	else
328	3.44k	{
329	3.44k	OGREnvelope oEnv;
330	3.44k	poGeometry->getEnvelope(&oEnv);
331	3.44k	padPtr[0] = oEnv.MinX;
332	3.44k	padPtr[1] = oEnv.MaxX;
333	3.44k	padPtr[2] = oEnv.MinY;
334	3.44k	padPtr[3] = oEnv.MaxY;
335	3.44k	}
336	3.99k	}
337
338	28.5k	GByte *pabyPtr = pabyWkb + nHeaderLen;
339
340		/* Use the wkbVariantIso for ISO SQL/MM output (differs for 3d geometry) */
341	28.5k	wkbExportOptions.eWkbVariant = wkbVariantIso;
342	28.5k	err = poGeometry->exportToWkb(pabyPtr, &wkbExportOptions);
343	28.5k	if (err != OGRERR_NONE)
344	0	{
345	0	CPLFree(pabyWkb);
346	0	return nullptr;
347	0	}
348
349	28.5k	return pabyWkb;
350	28.5k	}
351
352		OGRErr GPkgHeaderFromWKB(const GByte *pabyGpkg, size_t nGpkgLen,
353		GPkgHeader *poHeader)
354	76.4k	{
355	76.4k	CPLAssert(pabyGpkg != nullptr);
356	76.4k	CPLAssert(poHeader != nullptr);
357
358		/* Magic (match required) */
359	76.4k	if (nGpkgLen < 8 \|\| pabyGpkg[0] != 0x47 \|\| pabyGpkg[1] != 0x50 \|\|
360	75.4k	pabyGpkg[2] != 0) /* Version (only 0 supported at this time)*/
361	2.15k	{
362	2.15k	memset(poHeader, 0, sizeof(*poHeader));
363	2.15k	return OGRERR_FAILURE;
364	2.15k	}
365
366		/* Flags */
367	74.3k	GByte byFlags = pabyGpkg[3];
368	74.3k	poHeader->bEmpty = (byFlags & (0x01 << 4)) >> 4;
369	74.3k	poHeader->bExtended = (byFlags & (0x01 << 5)) >> 5;
370	74.3k	poHeader->eByteOrder = static_cast<OGRwkbByteOrder>(byFlags & 0x01);
371	74.3k	poHeader->bExtentHasXY = false;
372	74.3k	poHeader->bExtentHasZ = false;
373		#ifdef notdef
374		poHeader->bExtentHasM = false;
375		#endif
376	74.3k	bool bSwap = OGR_SWAP(poHeader->eByteOrder);
377
378		/* Envelope */
379	74.3k	int iEnvelope = (byFlags & (0x07 << 1)) >> 1;
380	74.3k	int nEnvelopeDim = 0;
381	74.3k	if (iEnvelope)
382	21.4k	{
383	21.4k	poHeader->bExtentHasXY = true;
384	21.4k	if (iEnvelope == 1)
385	18.0k	{
386	18.0k	nEnvelopeDim = 2; /* 2D envelope */
387	18.0k	}
388	3.38k	else if (iEnvelope == 2)
389	2.77k	{
390	2.77k	poHeader->bExtentHasZ = true;
391	2.77k	nEnvelopeDim = 3; /* 2D+Z envelope */
392	2.77k	}
393	615	else if (iEnvelope == 3)
394	0	{
395		#ifdef notdef
396		poHeader->bExtentHasM = true;
397		#endif
398	0	nEnvelopeDim = 3; /* 2D+M envelope */
399	0	}
400	615	else if (iEnvelope == 4)
401	615	{
402	615	poHeader->bExtentHasZ = true;
403		#ifdef notdef
404		poHeader->bExtentHasM = true;
405		#endif
406	615	nEnvelopeDim = 4; /* 2D+ZM envelope */
407	615	}
408	0	else
409	0	{
410	0	return OGRERR_FAILURE;
411	0	}
412	21.4k	}
413
414		/* SrsId */
415	74.3k	int iSrsId = 0;
416	74.3k	memcpy(&iSrsId, pabyGpkg + 4, 4);
417	74.3k	if (bSwap)
418	530	{
419	530	iSrsId = CPL_SWAP32(iSrsId);
420	530	}
421	74.3k	poHeader->iSrsId = iSrsId;
422
423	74.3k	if (nGpkgLen < static_cast<size_t>(8 + 8 * 2 * nEnvelopeDim))
424	1.47k	{
425		// Not enough bytes
426	1.47k	return OGRERR_FAILURE;
427	1.47k	}
428
429		/* Envelope */
430	72.8k	const double padPtr = reinterpret_cast<const double >(pabyGpkg + 8);
431	72.8k	if (poHeader->bExtentHasXY)
432	19.9k	{
433	19.9k	poHeader->MinX = padPtr[0];
434	19.9k	poHeader->MaxX = padPtr[1];
435	19.9k	poHeader->MinY = padPtr[2];
436	19.9k	poHeader->MaxY = padPtr[3];
437	19.9k	if (bSwap)
438	0	{
439	0	CPL_SWAPDOUBLE(&(poHeader->MinX));
440	0	CPL_SWAPDOUBLE(&(poHeader->MaxX));
441	0	CPL_SWAPDOUBLE(&(poHeader->MinY));
442	0	CPL_SWAPDOUBLE(&(poHeader->MaxY));
443	0	}
444	19.9k	}
445	72.8k	if (poHeader->bExtentHasZ)
446	2.76k	{
447	2.76k	poHeader->MinZ = padPtr[4];
448	2.76k	poHeader->MaxZ = padPtr[5];
449	2.76k	if (bSwap)
450	0	{
451	0	CPL_SWAPDOUBLE(&(poHeader->MinZ));
452	0	CPL_SWAPDOUBLE(&(poHeader->MaxZ));
453	0	}
454	2.76k	}
455		#ifdef notdef
456		if (poHeader->bExtentHasM)
457		{
458		poHeader->MinM = padPtr[(poHeader->bExtentHasZ) ? 6 : 4];
459		poHeader->MaxM = padPtr[(poHeader->bExtentHasZ) ? 7 : 5];
460		if (bSwap)
461		{
462		CPL_SWAPDOUBLE(&(poHeader->MinM));
463		CPL_SWAPDOUBLE(&(poHeader->MaxM));
464		}
465		}
466		#endif
467
468		/* Header size in byte stream */
469	72.8k	poHeader->nHeaderLen = 8 + 8 * 2 * nEnvelopeDim;
470
471		#ifdef DEBUG_VERBOSE
472		std::string s;
473		for (size_t i = poHeader->nHeaderLen; i < nGpkgLen; ++i)
474		{
475		s += CPLSPrintf("%02X ", pabyGpkg[i]);
476		}
477		CPLDebug("GPKG", "Bytes after GPKG header: %s", s.c_str());
478		#endif
479
480		// Workaround for a Spatialite bug. The CastToXYZ() function, when
481		// called on an empty geometry, and EnableGpkgMode() is enabled,
482		// returns an empty geometry, not tagged as such, but with an invalid
483		// bounding box.
484		// Cf https://github.com/OSGeo/gdal/issues/13557
485	72.8k	if (!poHeader->bEmpty && poHeader->bExtentHasXY &&
486	19.9k	poHeader->nHeaderLen == 40 &&
487	17.2k	poHeader->MinX == std::numeric_limits<double>::max() &&
488	0	poHeader->MaxX == -std::numeric_limits<double>::max() &&
489	0	poHeader->MinY == std::numeric_limits<double>::max() &&
490	0	poHeader->MaxY == -std::numeric_limits<double>::max() &&
491		// CastToXYZ(POLYGON EMPTY) returns just the GPKG header for some reason
492	0	(nGpkgLen == 40 \|\|
493		// CastToXYZ(LINESTRING EMPTY) returns a LINESTRING Z EMPTY
494	0	(nGpkgLen == 49 &&
495	0	memcmp(pabyGpkg + 40, "\x01\xEA\x03\x00\x00\x00\x00\x00\x00", 9) ==
496	0	0) \|\|
497		// CastToXYZ(POINT EMPTY) returns a Point Z (NaN NaN 0)
498		// CastToXYZ(POINT Z EMPTY) returns a Point Z (NaN NaN NaN)
499	0	(nGpkgLen == 69 && memcmp(pabyGpkg + 40,
500	0	"\x01"
501	0	"\xE9\x03\x00\x00"
502	0	"\x00\x00\x00\x00\x00\x00\xF8\x7F"
503	0	"\x00\x00\x00\x00\x00\x00\xF8\x7F",
504	0	21) == 0)))
505	0	{
506	0	CPLDebugOnce("GPKG",
507	0	"Work arounding a Spatialite bug with empty geometries");
508	0	poHeader->bEmpty = true;
509	0	poHeader->bExtentHasXY = false;
510	0	}
511
512	72.8k	return OGRERR_NONE;
513	74.3k	}
514
515		bool GPkgUpdateHeader(GByte *pabyGpkg, size_t nGpkgLen, int nSrsId, double MinX,
516		double MaxX, double MinY, double MaxY, double MinZ,
517		double MaxZ)
518	0	{
519	0	CPLAssert(nGpkgLen >= 8);
520
521		/* Flags */
522	0	const GByte byFlags = pabyGpkg[3];
523	0	const auto eByteOrder = static_cast<OGRwkbByteOrder>(byFlags & 0x01);
524	0	const bool bSwap = OGR_SWAP(eByteOrder);
525
526		/* SrsId */
527	0	if (bSwap)
528	0	{
529	0	nSrsId = CPL_SWAP32(nSrsId);
530	0	}
531	0	memcpy(pabyGpkg + 4, &nSrsId, 4);
532
533		/* Envelope */
534	0	const int iEnvelope = (byFlags & (0x07 << 1)) >> 1;
535	0	int nEnvelopeDim = 0;
536	0	if (iEnvelope)
537	0	{
538	0	if (iEnvelope == 1)
539	0	{
540	0	nEnvelopeDim = 2; /* 2D envelope */
541	0	}
542	0	else if (iEnvelope == 2)
543	0	{
544	0	nEnvelopeDim = 3; /* 2D+Z envelope */
545	0	}
546	0	else if (iEnvelope == 3)
547	0	{
548	0	nEnvelopeDim = 3; /* 2D+M envelope */
549	0	}
550	0	else if (iEnvelope == 4)
551	0	{
552	0	nEnvelopeDim = 4; /* 2D+ZM envelope */
553	0	}
554	0	else
555	0	{
556	0	return false;
557	0	}
558	0	}
559	0	else
560	0	{
561	0	return true;
562	0	}
563
564	0	if (nGpkgLen < static_cast<size_t>(8 + 8 * 2 * nEnvelopeDim))
565	0	{
566		// Not enough bytes
567	0	return false;
568	0	}
569
570		/* Envelope */
571	0	if (bSwap)
572	0	{
573	0	CPL_SWAPDOUBLE(&(MinX));
574	0	CPL_SWAPDOUBLE(&(MaxX));
575	0	CPL_SWAPDOUBLE(&(MinY));
576	0	CPL_SWAPDOUBLE(&(MaxY));
577	0	CPL_SWAPDOUBLE(&(MinZ));
578	0	CPL_SWAPDOUBLE(&(MaxZ));
579	0	}
580
581	0	double padPtr = reinterpret_cast<double >(pabyGpkg + 8);
582	0	memcpy(&padPtr[0], &MinX, sizeof(double));
583	0	memcpy(&padPtr[1], &MaxX, sizeof(double));
584	0	memcpy(&padPtr[2], &MinY, sizeof(double));
585	0	memcpy(&padPtr[3], &MaxY, sizeof(double));
586
587	0	if (iEnvelope == 2 \|\| iEnvelope == 4)
588	0	{
589	0	memcpy(&padPtr[4], &MinZ, sizeof(double));
590	0	memcpy(&padPtr[5], &MaxZ, sizeof(double));
591	0	}
592
593	0	return true;
594	0	}
595
596		OGRGeometry GPkgGeometryToOGR(const GByte pabyGpkg, size_t nGpkgLen,
597		OGRSpatialReference *poSrs)
598	7.40k	{
599	7.40k	CPLAssert(pabyGpkg != nullptr);
600
601	7.40k	GPkgHeader oHeader;
602
603		/* Read header */
604	7.40k	OGRErr err = GPkgHeaderFromWKB(pabyGpkg, nGpkgLen, &oHeader);
605	7.40k	if (err != OGRERR_NONE)
606	2.51k	return nullptr;
607
608		/* WKB pointer */
609	4.89k	const GByte *pabyWkb = pabyGpkg + oHeader.nHeaderLen;
610	4.89k	size_t nWkbLen = nGpkgLen - oHeader.nHeaderLen;
611
612		/* Parse WKB */
613	4.89k	OGRGeometry *poGeom = nullptr;
614	4.89k	err = OGRGeometryFactory::createFromWkb(pabyWkb, poSrs, &poGeom,
615	4.89k	static_cast<int>(nWkbLen));
616	4.89k	if (err != OGRERR_NONE)
617	219	return nullptr;
618
619	4.67k	return poGeom;
620	4.89k	}
621
622		/************************************************************************/
623		/* OGRGeoPackageGetHeader() */
624		/************************************************************************/
625
626		bool OGRGeoPackageGetHeader(sqlite3_context * /pContext/, int /argc/,
627		sqlite3_value *argv, GPkgHeader psHeader,
628		bool bNeedExtent, bool bNeedExtent3D, int iGeomIdx)
629	69.0k	{
630
631		// Extent3D implies extent
632	69.0k	const bool bNeedAnyExtent{bNeedExtent \|\| bNeedExtent3D};
633
634	69.0k	if (sqlite3_value_type(argv[iGeomIdx]) != SQLITE_BLOB)
635	0	{
636	0	memset(psHeader, 0, sizeof(*psHeader));
637	0	return false;
638	0	}
639	69.0k	const int nBLOBLen = sqlite3_value_bytes(argv[iGeomIdx]);
640	69.0k	const GByte *pabyBLOB =
641	69.0k	reinterpret_cast<const GByte *>(sqlite3_value_blob(argv[iGeomIdx]));
642
643	69.0k	if (nBLOBLen < 8)
644	0	{
645	0	memset(psHeader, 0, sizeof(*psHeader));
646	0	return false;
647	0	}
648	69.0k	else if (GPkgHeaderFromWKB(pabyBLOB, nBLOBLen, psHeader) != OGRERR_NONE)
649	1.11k	{
650	1.11k	bool bEmpty = false;
651	1.11k	memset(psHeader, 0, sizeof(*psHeader));
652	1.11k	if (OGRSQLiteGetSpatialiteGeometryHeader(
653	1.11k	pabyBLOB, nBLOBLen, &(psHeader->iSrsId), nullptr, &bEmpty,
654	1.11k	&(psHeader->MinX), &(psHeader->MinY), &(psHeader->MaxX),
655	1.11k	&(psHeader->MaxY)) == OGRERR_NONE)
656	0	{
657	0	psHeader->bEmpty = bEmpty;
658	0	psHeader->bExtentHasXY = !bEmpty;
659	0	if (!bNeedExtent3D && !(bEmpty && bNeedAnyExtent))
660	0	return true;
661	0	}
662
663	1.11k	return false;
664	1.11k	}
665
666	67.9k	if (psHeader->bEmpty && bNeedAnyExtent)
667	0	{
668	0	return false;
669	0	}
670	67.9k	else if (!psHeader->bExtentHasXY && bNeedExtent && !bNeedExtent3D)
671	28.0k	{
672	28.0k	OGREnvelope sEnvelope;
673	28.0k	if (OGRWKBGetBoundingBox(pabyBLOB + psHeader->nHeaderLen,
674	28.0k	static_cast<size_t>(nBLOBLen) -
675	28.0k	psHeader->nHeaderLen,
676	28.0k	sEnvelope))
677	28.0k	{
678	28.0k	psHeader->MinX = sEnvelope.MinX;
679	28.0k	psHeader->MaxX = sEnvelope.MaxX;
680	28.0k	psHeader->MinY = sEnvelope.MinY;
681	28.0k	psHeader->MaxY = sEnvelope.MaxY;
682	28.0k	return true;
683	28.0k	}
684	0	return false;
685	28.0k	}
686	39.9k	else if (!psHeader->bExtentHasZ && bNeedExtent3D)
687	0	{
688	0	OGREnvelope3D sEnvelope3D;
689	0	if (OGRWKBGetBoundingBox(pabyBLOB + psHeader->nHeaderLen,
690	0	static_cast<size_t>(nBLOBLen) -
691	0	psHeader->nHeaderLen,
692	0	sEnvelope3D))
693	0	{
694	0	psHeader->MinX = sEnvelope3D.MinX;
695	0	psHeader->MaxX = sEnvelope3D.MaxX;
696	0	psHeader->MinY = sEnvelope3D.MinY;
697	0	psHeader->MaxY = sEnvelope3D.MaxY;
698	0	psHeader->MinZ = sEnvelope3D.MinZ;
699	0	psHeader->MaxZ = sEnvelope3D.MaxZ;
700	0	return true;
701	0	}
702	0	return false;
703	0	}
704	39.9k	return true;
705	67.9k	}