/******************************************************************************

 *

 * Project:  FlatGeobuf driver

 * Purpose:  Implements OGRFlatGeobufLayer class.

 * Author:   Björn Harrtell <bjorn at wololo dot org>

 *

 ******************************************************************************

 * Copyright (c) 2018-2020, Björn Harrtell <bjorn at wololo dot org>

 *

 * SPDX-License-Identifier: MIT

 ****************************************************************************/

#include "ogrsf_frmts.h"

#include "cpl_vsi_virtual.h"

#include "cpl_conv.h"

#include "cpl_json.h"

#include "cpl_http.h"

#include "cpl_time.h"

#include "ogr_p.h"

#include "ograrrowarrayhelper.h"

#include "ogrlayerarrow.h"

#include "ogr_recordbatch.h"

#include "ogr_flatgeobuf.h"

#include "cplerrors.h"

#include "geometryreader.h"

#include "geometrywriter.h"

#include <algorithm>

#include <cmath>

#include <limits>

#include <new>

#include <stdexcept>

using namespace flatbuffers;

using namespace FlatGeobuf;

using namespace ogr_flatgeobuf;

static OGRErr CPLErrorMemoryAllocation(const char *message)

{

    CPLError(CE_Failure, CPLE_AppDefined, "Could not allocate memory: %s",

             message);

    return OGRERR_NOT_ENOUGH_MEMORY;

}

static OGRErr CPLErrorIO(const char *message)

{

    CPLError(CE_Failure, CPLE_AppDefined, "Unexpected I/O failure: %s",

             message);

    return OGRERR_FAILURE;

}

OGRFlatGeobufLayer::OGRFlatGeobufLayer(const Header *poHeader, GByte *headerBuf,

                                       const char *pszFilename, VSILFILE *poFp,

                                       uint64_t offset)

{

    m_poHeader = poHeader;

    CPLAssert(poHeader);

    m_headerBuf = headerBuf;

    CPLAssert(pszFilename);

    if (pszFilename)

        m_osFilename = pszFilename;

    m_poFp = poFp;

    m_offsetFeatures = offset;

    m_offset = offset;

    m_create = false;

    m_featuresCount = m_poHeader->features_count();

    m_geometryType = m_poHeader->geometry_type();

    m_indexNodeSize = m_poHeader->index_node_size();

    m_hasZ = m_poHeader->has_z();

    m_hasM = m_poHeader->has_m();

    m_hasT = m_poHeader->has_t();

    const auto envelope = m_poHeader->envelope();

    if (envelope && envelope->size() == 4 && std::isfinite((*envelope)[0]) &&

        std::isfinite((*envelope)[1]) && std::isfinite((*envelope)[2]) &&

        std::isfinite((*envelope)[3]))

    {

        m_sExtent.MinX = (*envelope)[0];

        m_sExtent.MinY = (*envelope)[1];

        m_sExtent.MaxX = (*envelope)[2];

        m_sExtent.MaxY = (*envelope)[3];

    }

    CPLDebugOnly("FlatGeobuf", "geometryType: %d, hasZ: %d, hasM: %d, hasT: %d",

                 (int)m_geometryType, m_hasZ, m_hasM, m_hasT);

    const auto crs = m_poHeader->crs();

    if (crs != nullptr)

    {

        m_poSRS = new OGRSpatialReference();

        m_poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);

        const auto org = crs->org();

        const auto code = crs->code();

        const auto crs_wkt = crs->wkt();

        CPLString wkt = crs_wkt ? crs_wkt->c_str() : "";

        double dfCoordEpoch = 0;

        if (STARTS_WITH_CI(wkt.c_str(), "COORDINATEMETADATA["))

        {

            size_t nPos = std::string::npos;

            // We don't want to match FRAMEEPOCH[

            for (const char *pszEpoch :

                 {",EPOCH[", " EPOCH[", "\tEPOCH[", "\nEPOCH[", "\rEPOCH["})

            {

                nPos = wkt.ifind(pszEpoch);

                if (nPos != std::string::npos)

                    break;

            }

            if (nPos != std::string::npos)

            {

                dfCoordEpoch = CPLAtof(wkt.c_str() + nPos + strlen(",EPOCH["));

                wkt.resize(nPos);

                wkt = wkt.substr(strlen("COORDINATEMETADATA["));

            }

        }

        if ((org == nullptr || EQUAL(org->c_str(), "EPSG")) && code != 0)

        {

            m_poSRS->importFromEPSG(code);

        }

        else if (org && code != 0)

        {

            CPLString osCode;

            osCode.Printf("%s:%d", org->c_str(), code);

            if (m_poSRS->SetFromUserInput(

                    osCode.c_str(),

                    OGRSpatialReference::

                        SET_FROM_USER_INPUT_LIMITATIONS_get()) != OGRERR_NONE &&

                !wkt.empty())

            {

                m_poSRS->importFromWkt(wkt.c_str());

            }

        }

        else if (!wkt.empty())

        {

            m_poSRS->importFromWkt(wkt.c_str());

        }

        if (dfCoordEpoch > 0)

            m_poSRS->SetCoordinateEpoch(dfCoordEpoch);

    }

    m_eGType = getOGRwkbGeometryType();

    if (const auto title = poHeader->title())

        SetMetadataItem("TITLE", title->c_str());

    if (const auto description = poHeader->description())

        SetMetadataItem("DESCRIPTION", description->c_str());

    if (const auto metadata = poHeader->metadata())

    {

        CPLJSONDocument oDoc;

        CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);

        if (oDoc.LoadMemory(metadata->c_str()) &&

            oDoc.GetRoot().GetType() == CPLJSONObject::Type::Object)

        {

            for (const auto &oItem : oDoc.GetRoot().GetChildren())

            {

                if (oItem.GetType() == CPLJSONObject::Type::String)

                {

                    SetMetadataItem(oItem.GetName().c_str(),

                                    oItem.ToString().c_str());

                }

            }

        }

    }

    const char *pszName =

        m_poHeader->name() ? m_poHeader->name()->c_str() : "unknown";

    m_poFeatureDefn = new OGRFeatureDefn(pszName);

    SetDescription(m_poFeatureDefn->GetName());

    m_poFeatureDefn->SetGeomType(wkbNone);

    auto poGeomFieldDefn =

        std::make_unique<OGRGeomFieldDefn>(nullptr, m_eGType);

    if (m_poSRS != nullptr)

        poGeomFieldDefn->SetSpatialRef(m_poSRS);

    m_poFeatureDefn->AddGeomFieldDefn(std::move(poGeomFieldDefn));

    readColumns();

    m_poFeatureDefn->Reference();

}

OGRFlatGeobufLayer::OGRFlatGeobufLayer(

    GDALDataset *poDS, const char *pszLayerName, const char *pszFilename,

    const OGRSpatialReference *poSpatialRef, OGRwkbGeometryType eGType,

    bool bCreateSpatialIndexAtClose, VSILFILE *poFpWrite,

    std::string &osTempFile, CSLConstList papszOptions)

    : m_eGType(eGType), m_poDS(poDS), m_create(true),

      m_bCreateSpatialIndexAtClose(bCreateSpatialIndexAtClose),

      m_poFpWrite(poFpWrite), m_aosCreationOption(papszOptions),

      m_osTempFile(osTempFile)

{

    if (pszLayerName)

        m_osLayerName = pszLayerName;

    if (pszFilename)

        m_osFilename = pszFilename;

    m_geometryType = GeometryWriter::translateOGRwkbGeometryType(eGType);

    if wkbHasZ (eGType)

        m_hasZ = true;

    if wkbHasM (eGType)

        m_hasM = true;

    if (poSpatialRef)

        m_poSRS = poSpatialRef->Clone();

    CPLDebugOnly("FlatGeobuf", "geometryType: %d, hasZ: %d, hasM: %d, hasT: %d",

                 (int)m_geometryType, m_hasZ, m_hasM, m_hasT);

    SetMetadataItem(OLMD_FID64, "YES");

    m_poFeatureDefn = new OGRFeatureDefn(pszLayerName);

    SetDescription(m_poFeatureDefn->GetName());

    m_poFeatureDefn->SetGeomType(eGType);

    m_poFeatureDefn->Reference();

}

OGRwkbGeometryType OGRFlatGeobufLayer::getOGRwkbGeometryType()

{

    OGRwkbGeometryType ogrType = OGRwkbGeometryType::wkbUnknown;

    if (static_cast<int>(m_geometryType) <= 17)

        ogrType = (OGRwkbGeometryType)m_geometryType;

    if (m_hasZ)

        ogrType = wkbSetZ(ogrType);

    if (m_hasM)

        ogrType = wkbSetM(ogrType);

    return ogrType;

}

static ColumnType toColumnType(const char *pszFieldName, OGRFieldType type,

                               OGRFieldSubType subType)

{

    switch (type)

    {

        case OGRFieldType::OFTInteger:

            return subType == OFSTBoolean ? ColumnType::Bool

                   : subType == OFSTInt16 ? ColumnType::Short

                                          : ColumnType::Int;

        case OGRFieldType::OFTInteger64:

            return ColumnType::Long;

        case OGRFieldType::OFTReal:

            return subType == OFSTFloat32 ? ColumnType::Float

                                          : ColumnType::Double;

        case OGRFieldType::OFTString:

            return ColumnType::String;

        case OGRFieldType::OFTDate:

            return ColumnType::DateTime;

        case OGRFieldType::OFTTime:

            return ColumnType::DateTime;

        case OGRFieldType::OFTDateTime:

            return ColumnType::DateTime;

        case OGRFieldType::OFTBinary:

            return ColumnType::Binary;

        default:

            CPLError(CE_Warning, CPLE_AppDefined,

                     "toColumnType: %s field is of type %s, which is not "

                     "handled natively. Falling back to String.",

                     pszFieldName, OGRFieldDefn::GetFieldTypeName(type));

    }

    return ColumnType::String;

}

static OGRFieldType toOGRFieldType(ColumnType type, OGRFieldSubType &eSubType)

{

    eSubType = OFSTNone;

    switch (type)

    {

        case ColumnType::Byte:

            return OGRFieldType::OFTInteger;

        case ColumnType::UByte:

            return OGRFieldType::OFTInteger;

        case ColumnType::Bool:

            eSubType = OFSTBoolean;

            return OGRFieldType::OFTInteger;

        case ColumnType::Short:

            eSubType = OFSTInt16;

            return OGRFieldType::OFTInteger;

        case ColumnType::UShort:

            return OGRFieldType::OFTInteger;

        case ColumnType::Int:

            return OGRFieldType::OFTInteger;

        case ColumnType::UInt:

            return OGRFieldType::OFTInteger64;

        case ColumnType::Long:

            return OGRFieldType::OFTInteger64;

        case ColumnType::ULong:

            return OGRFieldType::OFTReal;

        case ColumnType::Float:

            eSubType = OFSTFloat32;

            return OGRFieldType::OFTReal;

        case ColumnType::Double:

            return OGRFieldType::OFTReal;

        case ColumnType::String:

            return OGRFieldType::OFTString;

        case ColumnType::Json:

            return OGRFieldType::OFTString;

        case ColumnType::DateTime:

            return OGRFieldType::OFTDateTime;

        case ColumnType::Binary:

            return OGRFieldType::OFTBinary;

    }

    return OGRFieldType::OFTString;

}

const std::vector<Offset<Column>>

OGRFlatGeobufLayer::writeColumns(FlatBufferBuilder &fbb)

{

    std::vector<Offset<Column>> columns;

    for (int i = 0; i < m_poFeatureDefn->GetFieldCount(); i++)

    {

        const auto field = m_poFeatureDefn->GetFieldDefn(i);

        const auto name = field->GetNameRef();

        const auto columnType =

            toColumnType(name, field->GetType(), field->GetSubType());

        auto title = field->GetAlternativeNameRef();

        if (EQUAL(title, ""))

            title = nullptr;

        const std::string &osComment = field->GetComment();

        const char *description =

            !osComment.empty() ? osComment.c_str() : nullptr;

        auto width = -1;

        auto precision = -1;

        auto scale = field->GetPrecision();

        if (scale == 0)

            scale = -1;

        if (columnType == ColumnType::Float || columnType == ColumnType::Double)

            precision = field->GetWidth();

        else

            width = field->GetWidth();

        auto nullable = CPL_TO_BOOL(field->IsNullable());

        auto unique = CPL_TO_BOOL(field->IsUnique());

        auto primaryKey = false;

        // CPLDebugOnly("FlatGeobuf", "Create column %s (index %d)", name, i);

        const auto column =

            CreateColumnDirect(fbb, name, columnType, title, description, width,

                               precision, scale, nullable, unique, primaryKey);

        columns.push_back(column);

        // CPLDebugOnly("FlatGeobuf", "DEBUG writeColumns: Created column %s

        // added as index %d", name, i);

    }

    CPLDebugOnly("FlatGeobuf", "Created %lu columns for writing",

                 static_cast<long unsigned int>(columns.size()));

    return columns;

}

void OGRFlatGeobufLayer::readColumns()

{

    const auto columns = m_poHeader->columns();

    if (columns == nullptr)

        return;

    for (uint32_t i = 0; i < columns->size(); i++)

    {

        const auto column = columns->Get(i);

        const auto type = column->type();

        const auto name = column->name()->c_str();

        const auto title =

            column->title() != nullptr ? column->title()->c_str() : nullptr;

        const auto width = column->width();

        const auto precision = column->precision();

        const auto scale = column->scale();

        const auto nullable = column->nullable();

        const auto unique = column->unique();

        OGRFieldSubType eSubType = OFSTNone;

        const auto ogrType = toOGRFieldType(column->type(), eSubType);

        OGRFieldDefn field(name, ogrType);

        field.SetSubType(eSubType);

        field.SetAlternativeName(title);

        if (column->description())

            field.SetComment(column->description()->str());

        if (width != -1 && type != ColumnType::Float &&

            type != ColumnType::Double)

            field.SetWidth(width);

        if (precision != -1)

            field.SetWidth(precision);

        field.SetPrecision(scale != -1 ? scale : 0);

        field.SetNullable(nullable);

        field.SetUnique(unique);

        m_poFeatureDefn->AddFieldDefn(&field);

        // CPLDebugOnly("FlatGeobuf", "DEBUG readColumns: Read column %s added

        // as index %d", name, i);

    }

    CPLDebugOnly("FlatGeobuf",

                 "Read %lu columns and added to feature definition",

                 static_cast<long unsigned int>(columns->size()));

}

void OGRFlatGeobufLayer::writeHeader(VSILFILE *poFp, uint64_t featuresCount,

                                     std::vector<double> *extentVector)

{

    size_t c;

    c = VSIFWriteL(&magicbytes, sizeof(magicbytes), 1, poFp);

    CPLDebugOnly("FlatGeobuf", "Wrote magicbytes (%lu bytes)",

                 static_cast<long unsigned int>(c * sizeof(magicbytes)));

    m_writeOffset += sizeof(magicbytes);

    FlatBufferBuilder fbb;

    fbb.TrackMinAlign(8);

    auto columns = writeColumns(fbb);

    flatbuffers::Offset<Crs> crs = 0;

    if (m_poSRS)

    {

        int nAuthorityCode = 0;

        const char *pszAuthorityName = m_poSRS->GetAuthorityName();

        if (pszAuthorityName == nullptr || strlen(pszAuthorityName) == 0)

        {

            // Try to force identify an EPSG code.

            m_poSRS->AutoIdentifyEPSG();

            pszAuthorityName = m_poSRS->GetAuthorityName();

            if (pszAuthorityName != nullptr && EQUAL(pszAuthorityName, "EPSG"))

            {

                const char *pszAuthorityCode = m_poSRS->GetAuthorityCode();

                if (pszAuthorityCode != nullptr && strlen(pszAuthorityCode) > 0)

                {

                    /* Import 'clean' SRS */

                    m_poSRS->importFromEPSG(atoi(pszAuthorityCode));

                    pszAuthorityName = m_poSRS->GetAuthorityName();

                }

            }

        }

        if (pszAuthorityName != nullptr && strlen(pszAuthorityName) > 0)

        {

            // For the root authority name 'EPSG', the authority code

            // should always be integral

            nAuthorityCode = atoi(m_poSRS->GetAuthorityCode());

        }

        // Translate SRS to WKT.

        char *pszWKT = nullptr;

        const char *const apszOptionsWkt[] = {"FORMAT=WKT2_2019", nullptr};

        m_poSRS->exportToWkt(&pszWKT, apszOptionsWkt);

        if (pszWKT && pszWKT[0] == '\0')

        {

            CPLFree(pszWKT);

            pszWKT = nullptr;

        }

        if (pszWKT && m_poSRS->GetCoordinateEpoch() > 0)

        {

            std::string osCoordinateEpoch =

                CPLSPrintf("%f", m_poSRS->GetCoordinateEpoch());

            if (osCoordinateEpoch.find('.') != std::string::npos)

            {

                while (osCoordinateEpoch.back() == '0')

                    osCoordinateEpoch.pop_back();

            }

            std::string osWKT("COORDINATEMETADATA[");

            osWKT += pszWKT;

            osWKT += ",EPOCH[";

            osWKT += osCoordinateEpoch;

            osWKT += "]]";

            CPLFree(pszWKT);

            pszWKT = CPLStrdup(osWKT.c_str());

        }

        if (pszWKT && !CPLIsUTF8(pszWKT, -1))

        {

            char *pszWKTtmp = CPLForceToASCII(pszWKT, -1, '?');

            CPLFree(pszWKT);

            pszWKT = pszWKTtmp;

        }

        crs = CreateCrsDirect(fbb, pszAuthorityName, nAuthorityCode,

                              m_poSRS->GetName(), nullptr, pszWKT);

        CPLFree(pszWKT);

    }

    std::string osTitle(m_aosCreationOption.FetchNameValueDef("TITLE", ""));

    std::string osDescription(

        m_aosCreationOption.FetchNameValueDef("DESCRIPTION", ""));

    std::string osMetadata;

    CPLJSONObject oMetadataJSONObj;

    bool bEmptyMetadata = true;

    for (GDALMajorObject *poContainer :

         {static_cast<GDALMajorObject *>(this),

          static_cast<GDALMajorObject *>(

              m_poDS && m_poDS->GetLayerCount() == 1 ? m_poDS : nullptr)})

    {

        if (poContainer)

        {

            if (CSLConstList papszMD = poContainer->GetMetadata())

            {

                for (CSLConstList papszIter = papszMD; *papszIter; ++papszIter)

                {

                    char *pszKey = nullptr;

                    const char *pszValue =

                        CPLParseNameValue(*papszIter, &pszKey);

                    if (pszKey && pszValue && !EQUAL(pszKey, OLMD_FID64))

                    {

                        if (EQUAL(pszKey, "TITLE"))

                        {

                            if (osTitle.empty())

                                osTitle = pszValue;

                        }

                        else if (EQUAL(pszKey, "DESCRIPTION"))

                        {

                            if (osDescription.empty())

                                osDescription = pszValue;

                        }

                        else

                        {

                            bEmptyMetadata = false;

                            oMetadataJSONObj.Add(pszKey, pszValue);

                        }

                    }

                    CPLFree(pszKey);

                }

            }

        }

    }

    if (!bEmptyMetadata)

    {

        osMetadata =

            oMetadataJSONObj.Format(CPLJSONObject::PrettyFormat::Plain);

    }

    const auto header = CreateHeaderDirect(

        fbb, m_osLayerName.c_str(), extentVector, m_geometryType, m_hasZ,

        m_hasM, m_hasT, m_hasTM, &columns, featuresCount, m_indexNodeSize, crs,

        osTitle.empty() ? nullptr : osTitle.c_str(),

        osDescription.empty() ? nullptr : osDescription.c_str(),

        osMetadata.empty() ? nullptr : osMetadata.c_str());

    fbb.FinishSizePrefixed(header);

    c = VSIFWriteL(fbb.GetBufferPointer(), 1, fbb.GetSize(), poFp);

    CPLDebugOnly("FlatGeobuf", "Wrote header (%lu bytes)",

                 static_cast<long unsigned int>(c));

    m_writeOffset += c;

}

static bool SupportsSeekWhileWriting(const std::string &osFilename)

{

    return (!STARTS_WITH(osFilename.c_str(), "/vsi")) ||

           STARTS_WITH(osFilename.c_str(), "/vsimem/");

}

bool OGRFlatGeobufLayer::CreateFinalFile()

{

    // no spatial index requested, we are (almost) done

    if (!m_bCreateSpatialIndexAtClose)

    {

        if (m_poFpWrite == nullptr || !SupportsSeekWhileWriting(m_osFilename))

        {

            return true;

        }

        // Rewrite header

        VSIFSeekL(m_poFpWrite, 0, SEEK_SET);

        m_writeOffset = 0;

        std::vector<double> extentVector;

        if (!m_sExtent.IsInit())

        {

            extentVector.resize(4, std::numeric_limits<double>::quiet_NaN());

        }

        else

        {

            extentVector.push_back(m_sExtent.MinX);

            extentVector.push_back(m_sExtent.MinY);

            extentVector.push_back(m_sExtent.MaxX);

            extentVector.push_back(m_sExtent.MaxY);

        }

        writeHeader(m_poFpWrite, m_featuresCount, &extentVector);

        // Sanity check to verify that the dummy header and the real header

        // have the same size.

        if (m_featuresCount)

        {

            CPLAssert(m_writeOffset == m_offsetAfterHeader);

        }

        CPL_IGNORE_RET_VAL(m_writeOffset);  // otherwise checkers might tell the

                                            // member is not used

        return true;

    }

    m_poFp = VSIFOpenL(m_osFilename.c_str(), "wb");

    if (m_poFp == nullptr)

    {

        CPLError(CE_Failure, CPLE_OpenFailed, "Failed to create %s:\n%s",

                 m_osFilename.c_str(), VSIStrerror(errno));

        return false;

    }

    // check if something has been written, if not write empty layer and bail

    if (m_writeOffset == 0 || m_featuresCount == 0)

    {

        CPLDebugOnly("FlatGeobuf", "Writing empty layer");

        writeHeader(m_poFp, 0, nullptr);

        return true;

    }

    CPLDebugOnly("FlatGeobuf", "Writing second pass sorted by spatial index");

    const uint64_t nTempFileSize = m_writeOffset;

    m_writeOffset = 0;

    m_indexNodeSize = 16;

    size_t c;

    if (m_featuresCount >= std::numeric_limits<size_t>::max() / 8)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "Too many features for this architecture");

        return false;

    }

    NodeItem extent = calcExtent(m_featureItems);

    auto extentVector = extent.toVector();

    writeHeader(m_poFp, m_featuresCount, &extentVector);

    CPLDebugOnly("FlatGeobuf", "Sorting items for Packed R-tree");

    hilbertSort(m_featureItems);

    CPLDebugOnly("FlatGeobuf", "Calc new feature offsets");

    uint64_t featureOffset = 0;

    for (auto &item : m_featureItems)

    {

        item.nodeItem.offset = featureOffset;

        featureOffset += item.size;

    }

    CPLDebugOnly("FlatGeobuf", "Creating Packed R-tree");

    c = 0;

    try

    {

        const auto fillNodeItems = [this](NodeItem *dest)

        {

            size_t i = 0;

            for (const auto &featureItem : m_featureItems)

            {

                dest[i] = featureItem.nodeItem;

                ++i;

            }

        };

        PackedRTree tree(fillNodeItems, m_featureItems.size(), extent);

        CPLDebugOnly("FlatGeobuf", "PackedRTree extent %f, %f, %f, %f",

                     extentVector[0], extentVector[1], extentVector[2],

                     extentVector[3]);

        tree.streamWrite([this, &c](uint8_t *data, size_t size)

                         { c += VSIFWriteL(data, 1, size, m_poFp); });

    }

    catch (const std::exception &e)

    {

        CPLError(CE_Failure, CPLE_AppDefined, "Create: %s", e.what());

        return false;

    }

    CPLDebugOnly("FlatGeobuf", "Wrote tree (%lu bytes)",

                 static_cast<long unsigned int>(c));

    m_writeOffset += c;

    CPLDebugOnly("FlatGeobuf", "Writing feature buffers at offset %lu",

                 static_cast<long unsigned int>(m_writeOffset));

    c = 0;

    // For temporary files not in memory, we use a batch strategy to write the

    // final file. That is to say we try to separate reads in the source

    // temporary file and writes in the target file as much as possible, and by

    // reading source features in increasing offset within a batch.

    const bool bUseBatchStrategy =

        !STARTS_WITH(m_osTempFile.c_str(), "/vsimem/");

    if (bUseBatchStrategy)

    {

        const uint32_t nMaxBufferSize = std::max(

            m_maxFeatureSize,

            static_cast<uint32_t>(std::min(

                static_cast<uint64_t>(100 * 1024 * 1024), nTempFileSize)));

        if (ensureFeatureBuf(nMaxBufferSize) != OGRERR_NONE)

            return false;

        uint32_t offsetInBuffer = 0;

        struct BatchItem

        {

            size_t featureIdx;  // index of m_featureItems[]

            uint32_t offsetInBuffer;

        };

        std::vector<BatchItem> batch;

        const auto flushBatch = [this, &batch, &offsetInBuffer]()

        {

            // Sort by increasing source offset

            std::sort(batch.begin(), batch.end(),

                      [this](const BatchItem &a, const BatchItem &b)

                      {

                          return m_featureItems[a.featureIdx].offset <

                                 m_featureItems[b.featureIdx].offset;

                      });

            // Read source features

            for (const auto &batchItem : batch)

            {

                const auto &item = m_featureItems[batchItem.featureIdx];

                if (VSIFSeekL(m_poFpWrite, item.offset, SEEK_SET) == -1)

                {

                    CPLErrorIO("seeking to temp feature location");

                    return false;

                }

                if (VSIFReadL(m_featureBuf + batchItem.offsetInBuffer, 1,

                              item.size, m_poFpWrite) != item.size)

                {

                    CPLErrorIO("reading temp feature");

                    return false;

                }

            }

            // Write target features

            if (offsetInBuffer > 0 &&

                VSIFWriteL(m_featureBuf, 1, offsetInBuffer, m_poFp) !=

                    offsetInBuffer)

            {

                CPLErrorIO("writing feature");

                return false;

            }

            batch.clear();

            offsetInBuffer = 0;

            return true;

        };

        for (size_t i = 0; i < m_featuresCount; i++)

        {

            const auto &featureItem = m_featureItems[i];

            const auto featureSize = featureItem.size;

            if (offsetInBuffer + featureSize > m_featureBufSize)

            {

                if (!flushBatch())

                {

                    return false;

                }

            }

            BatchItem bachItem;

            bachItem.offsetInBuffer = offsetInBuffer;

            bachItem.featureIdx = i;

            batch.emplace_back(bachItem);

            offsetInBuffer += featureSize;

            c += featureSize;

        }

        if (!flushBatch())

        {

            return false;

        }

    }

    else

    {

        const auto err = ensureFeatureBuf(m_maxFeatureSize);

        if (err != OGRERR_NONE)

            return false;

        for (const auto &featureItem : m_featureItems)

        {

            const auto featureSize = featureItem.size;

            // CPLDebugOnly("FlatGeobuf", "featureItem.offset: %lu",

            // static_cast<long unsigned int>(featureItem.offset));

            // CPLDebugOnly("FlatGeobuf", "featureSize: %d", featureSize);

            if (VSIFSeekL(m_poFpWrite, featureItem.offset, SEEK_SET) == -1)

            {

                CPLErrorIO("seeking to temp feature location");

                return false;

            }

            if (VSIFReadL(m_featureBuf, 1, featureSize, m_poFpWrite) !=

                featureSize)

            {

                CPLErrorIO("reading temp feature");

                return false;

            }

            if (VSIFWriteL(m_featureBuf, 1, featureSize, m_poFp) != featureSize)

            {

                CPLErrorIO("writing feature");

                return false;

            }

            c += featureSize;

        }

    }

    CPLDebugOnly("FlatGeobuf", "Wrote feature buffers (%lu bytes)",

                 static_cast<long unsigned int>(c));

    m_writeOffset += c;

    CPLDebugOnly("FlatGeobuf", "Now at offset %lu",

                 static_cast<long unsigned int>(m_writeOffset));

    return true;

}

OGRFlatGeobufLayer::~OGRFlatGeobufLayer()

{

    OGRFlatGeobufLayer::Close();

    if (m_poFeatureDefn)

        m_poFeatureDefn->Release();

    if (m_poSRS)

        m_poSRS->Release();

    if (m_featureBuf)

        VSIFree(m_featureBuf);

    if (m_headerBuf)

        VSIFree(m_headerBuf);

}

CPLErr OGRFlatGeobufLayer::Close(GDALProgressFunc, void *)

{

    CPLErr eErr = CE_None;

    if (m_create)

    {

        if (!CreateFinalFile())

            eErr = CE_Failure;

        m_create = false;

    }

    if (m_poFp)

    {

        if (VSIFCloseL(m_poFp) != 0)

            eErr = CE_Failure;

        m_poFp = nullptr;

    }

    if (m_poFpWrite)

    {

        if (VSIFCloseL(m_poFpWrite) != 0)

            eErr = CE_Failure;

        m_poFpWrite = nullptr;

    }

    if (!m_osTempFile.empty())

    {

        VSIUnlink(m_osTempFile.c_str());

        m_osTempFile.clear();

    }

    return eErr;

}

OGRErr OGRFlatGeobufLayer::readFeatureOffset(uint64_t index,

                                             uint64_t &featureOffset)

{

    try

    {

        const auto treeSize =

            PackedRTree::size(m_featuresCount, m_indexNodeSize);

        const auto levelBounds =

            PackedRTree::generateLevelBounds(m_featuresCount, m_indexNodeSize);

        const auto bottomLevelOffset =

            m_offset - treeSize +

            (levelBounds.front().first * sizeof(NodeItem));

        const auto nodeItemOffset =

            bottomLevelOffset + (index * sizeof(NodeItem));

        const auto featureOffsetOffset = nodeItemOffset + (sizeof(double) * 4);

        if (VSIFSeekL(m_poFp, featureOffsetOffset, SEEK_SET) == -1)

            return CPLErrorIO("seeking feature offset");

        if (VSIFReadL(&featureOffset, sizeof(uint64_t), 1, m_poFp) != 1)

            return CPLErrorIO("reading feature offset");

        CPL_LSBPTR64(&featureOffset);

        return OGRERR_NONE;

    }

    catch (const std::exception &e)

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "Failed to calculate tree size: %s", e.what());

        return OGRERR_FAILURE;

    }

}

OGRFeature *OGRFlatGeobufLayer::GetFeature(GIntBig nFeatureId)

{

    if (m_indexNodeSize == 0)

    {

        return OGRLayer::GetFeature(nFeatureId);

    }

    else

    {

        if (nFeatureId < 0 ||

            static_cast<uint64_t>(nFeatureId) >= m_featuresCount)

        {

            return nullptr;

        }

        ResetReading();

        m_ignoreSpatialFilter = true;

        m_ignoreAttributeFilter = true;

        uint64_t featureOffset;

        const auto err = readFeatureOffset(nFeatureId, featureOffset);

        if (err != OGRERR_NONE)

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "Unexpected error reading feature offset from id");

            return nullptr;

        }

        m_offset = m_offsetFeatures + featureOffset;

        OGRFeature *poFeature = GetNextFeature();

        if (poFeature != nullptr)

            poFeature->SetFID(nFeatureId);

        ResetReading();

        return poFeature;

    }

}

OGRErr OGRFlatGeobufLayer::readIndex()

{

    if (m_queriedSpatialIndex || !m_poFilterGeom)

        return OGRERR_NONE;

    if (m_sFilterEnvelope.IsInit() && m_sExtent.IsInit() &&

        m_sFilterEnvelope.MinX <= m_sExtent.MinX &&

        m_sFilterEnvelope.MinY <= m_sExtent.MinY &&

        m_sFilterEnvelope.MaxX >= m_sExtent.MaxX &&

        m_sFilterEnvelope.MaxY >= m_sExtent.MaxY)

        return OGRERR_NONE;

    const auto indexNodeSize = m_poHeader->index_node_size();

    if (indexNodeSize == 0)

        return OGRERR_NONE;

    const auto featuresCount = m_poHeader->features_count();

    if (featuresCount == 0)

        return OGRERR_NONE;

    if (VSIFSeekL(m_poFp, static_cast<vsi_l_offset>(sizeof(magicbytes)),

                  SEEK_SET) == -1)  // skip magic bytes

        return CPLErrorIO("seeking past magic bytes");

    uoffset_t headerSize;

    if (VSIFReadL(&headerSize, sizeof(uoffset_t), 1, m_poFp) != 1)

        return CPLErrorIO("reading header size");

    CPL_LSBPTR32(&headerSize);

    try

    {

        const auto treeSize =

            indexNodeSize > 0 ? PackedRTree::size(featuresCount) : 0;

        if (treeSize > 0 && m_poFilterGeom && !m_ignoreSpatialFilter)

        {

            CPLDebugOnly("FlatGeobuf", "Attempting spatial index query");

            OGREnvelope env;

            m_poFilterGeom->getEnvelope(&env);

            NodeItem n{env.MinX, env.MinY, env.MaxX, env.MaxY, 0};

            CPLDebugOnly("FlatGeobuf", "Spatial index search on %f,%f,%f,%f",

                         env.MinX, env.MinY, env.MaxX, env.MaxY);

            const vsi_l_offset treeOffset =

                sizeof(magicbytes) + sizeof(uoffset_t) + headerSize;

            const auto readNode =

                [this, treeOffset](uint8_t *buf, size_t i, size_t s)

            {

                if (VSIFSeekL(m_poFp, treeOffset + i, SEEK_SET) == -1)

                    throw std::runtime_error("I/O seek failure");

                if (VSIFReadL(buf, 1, s, m_poFp) != s)

                    throw std::runtime_error("I/O read file");

            };

            m_foundItems = PackedRTree::streamSearch(

                featuresCount, indexNodeSize, n, readNode);

            m_featuresCount = m_foundItems.size();

            CPLDebugOnly("FlatGeobuf",

                         "%lu features found in spatial index search",