/src/gdal/frmts/netcdf/netcdfsg.cpp

Source
/******************************************************************************
 *
 * Project:  netCDF read/write Driver
 * Purpose:  GDAL bindings over netCDF library.
 * Author:   Winor Chen <wchen329 at wisc.edu>
 *
 ******************************************************************************
 * Copyright (c) 2019, Winor Chen <wchen329 at wisc.edu>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/
#include <cstdio>
#include <cstring>
#include <set>
#include <vector>
#include "netcdf.h"
#include "netcdfdataset.h"
#include "netcdfsg.h"

namespace nccfdriver
{

/* SGeometry_Reader
 * (implementations)
 *
 */
SGeometry_Reader::SGeometry_Reader(int ncId, int geoVarId)
    : gc_varId(geoVarId), touple_order(0)
{

    char container_name[NC_MAX_NAME + 1];
    memset(container_name, 0, NC_MAX_NAME + 1);

    // Get geometry container name
    if (nc_inq_varname(ncId, geoVarId, container_name) != NC_NOERR)
    {
        throw SG_Exception_Existential("new geometry container",
                                       "the variable of the given ID");
    }

    // Establish string version of container_name
    container_name_s = std::string(container_name);

    // Find geometry type
    this->type = nccfdriver::getGeometryType(ncId, geoVarId);

    if (this->type == NONE)
    {
        throw SG_Exception_Existential(
            static_cast<const char *>(container_name), CF_SG_GEOMETRY_TYPE);
    }

    // Get grid mapping variable, if it exists
    this->gm_varId = INVALID_VAR_ID;
    if (NCDFGetAttr(ncId, geoVarId, CF_GRD_MAPPING, gm_name_s) == CE_None &&
        !gm_name_s.empty())
    {
        const char *gm_name = gm_name_s.c_str();
        int gmVID;
        if (nc_inq_varid(ncId, gm_name, &gmVID) == NC_NOERR)
        {
            this->gm_varId = gmVID;
        }
    }

    // Find a list of node counts and part node count
    std::string nc_name_s;
    std::string pnc_name_s;
    std::string ir_name_s;
    int pnc_vid = INVALID_VAR_ID;
    int nc_vid = INVALID_VAR_ID;
    int ir_vid = INVALID_VAR_ID;
    int buf;
    size_t bound = 0;
    size_t total_node_count = 0;  // used in error checks later
    size_t total_part_node_count = 0;

    if (NCDFGetAttr(ncId, geoVarId, CF_SG_NODE_COUNT, nc_name_s) == CE_None &&
        !nc_name_s.empty())
    {
        const char *nc_name = nc_name_s.c_str();
        nc_inq_varid(ncId, nc_name, &nc_vid);
        while (nc_get_var1_int(ncId, nc_vid, &bound, &buf) == NC_NOERR)
        {
            if (buf < 0)
            {
                throw SG_Exception_Value_Violation(
                    static_cast<const char *>(container_name), CF_SG_NODE_COUNT,
                    "negative");
            }

            this->node_counts.push_back(buf);
            total_node_count += buf;
            bound++;
        }
    }

    if (NCDFGetAttr(ncId, geoVarId, CF_SG_PART_NODE_COUNT, pnc_name_s) ==
            CE_None &&
        !pnc_name_s.empty())
    {
        const char *pnc_name = pnc_name_s.c_str();
        bound = 0;
        nc_inq_varid(ncId, pnc_name, &pnc_vid);
        while (nc_get_var1_int(ncId, pnc_vid, &bound, &buf) == NC_NOERR)
        {
            if (buf < 0)
            {
                throw SG_Exception_Value_Violation(
                    static_cast<const char *>(container_name),
                    CF_SG_PART_NODE_COUNT, "negative");
            }

            this->pnode_counts.push_back(buf);
            total_part_node_count += buf;
            bound++;
        }
    }

    if (NCDFGetAttr(ncId, geoVarId, CF_SG_INTERIOR_RING, ir_name_s) ==
            CE_None &&
        !ir_name_s.empty())
    {
        const char *ir_name = ir_name_s.c_str();
        bound = 0;
        nc_inq_varid(ncId, ir_name, &ir_vid);
        while (nc_get_var1_int(ncId, ir_vid, &bound, &buf) == NC_NOERR)
        {
            bool store = buf == 0 ? false : true;

            if (buf != 0 && buf != 1)
            {
                throw SG_Exception_Value_Required(
                    static_cast<const char *>(container_name),
                    CF_SG_INTERIOR_RING, "0 or 1");
            }

            this->int_rings.push_back(store);
            bound++;
        }
    }

    /* Enforcement of well formed CF files
     * If these are not met then the dataset is malformed and will halt further
     * processing of simple geometries.
     */

    // part node count exists only when node count exists
    if (pnode_counts.size() > 0 && node_counts.size() == 0)
    {
        throw SG_Exception_Dep(static_cast<const char *>(container_name),
                               CF_SG_PART_NODE_COUNT, CF_SG_NODE_COUNT);
    }

    // part node counts and node counts don't match up in count
    if (pnc_vid != INVALID_VAR_ID && total_node_count != total_part_node_count)
    {
        throw SG_Exception_BadSum(static_cast<const char *>(container_name),
                                  CF_SG_PART_NODE_COUNT, CF_SG_PART_NODE_COUNT);
    }

    // interior rings only exist when part node counts exist
    if (int_rings.size() > 0 && pnode_counts.size() == 0)
    {
        throw SG_Exception_Dep(static_cast<const char *>(container_name),
                               CF_SG_INTERIOR_RING, CF_SG_PART_NODE_COUNT);
    }

    // cardinality of part_node_counts == cardinality of interior_ring (if
    // interior ring > 0)
    if (int_rings.size() > 0)
    {
        if (int_rings.size() != pnode_counts.size())
        {
            throw SG_Exception_Dim_MM(static_cast<const char *>(container_name),
                                      CF_SG_INTERIOR_RING,
                                      CF_SG_PART_NODE_COUNT);
        }
    }

    // lines and polygons require node_counts, multipolygons are checked with
    // part_node_count
    if (this->type == POLYGON || this->type == LINE)
    {
        if (node_counts.size() < 1)
        {
            throw SG_Exception_Existential(
                static_cast<const char *>(container_name), CF_SG_NODE_COUNT);
        }
    }

    /* Basic Safety checks End
     */

    // Create bound list
    size_t rc = 0;
    bound_list.push_back(0);  // start with 0

    if (node_counts.size() > 0)
    {
        for (size_t i = 0; i < node_counts.size() - 1; i++)
        {
            rc = rc + node_counts[i];
            bound_list.push_back(rc);
        }
    }

    std::string cart_s;
    // Node Coordinates
    if (NCDFGetAttr(ncId, geoVarId, CF_SG_NODE_COORDINATES, cart_s) !=
            CE_None ||
        cart_s.empty())
    {
        throw SG_Exception_Existential(container_name, CF_SG_NODE_COORDINATES);
    }

    // Create parts count list and an offset list for parts indexing
    if (this->node_counts.size() > 0)
    {
        int ind = 0;
        int parts = 0;
        int prog = 0;
        int c = 0;

        for (size_t pcnt = 0; pcnt < pnode_counts.size(); pcnt++)
        {
            if (prog == 0)
                pnc_bl.push_back(pcnt);

            if (int_rings.size() > 0 && !int_rings[pcnt])
                c++;

            prog = prog + pnode_counts[pcnt];
            parts++;

            if (prog == node_counts[ind])
            {
                ind++;
                this->parts_count.push_back(parts);
                if (int_rings.size() > 0)
                    this->poly_count.push_back(c);
                c = 0;
                prog = 0;
                parts = 0;
            }
            else if (prog > node_counts[ind])
            {
                throw SG_Exception_BadSum(container_name, CF_SG_PART_NODE_COUNT,
                                          CF_SG_NODE_COUNT);
            }
        }
    }

    // (1) the tuple order for a single point
    // (2) the variable ids with the relevant coordinate values
    int X = INVALID_VAR_ID;
    int Y = INVALID_VAR_ID;
    int Z = INVALID_VAR_ID;

    char cart[NC_MAX_NAME + 1];
    memset(cart, 0, NC_MAX_NAME + 1);
    strncpy(cart, cart_s.c_str(), NC_MAX_NAME);

    char *dim = strtok(cart, " ");
    int axis_id = 0;

    while (dim != nullptr)
    {
        if (nc_inq_varid(ncId, dim, &axis_id) == NC_NOERR)
        {

            // Check axis signature
            std::string a_sig;
            CPL_IGNORE_RET_VAL(NCDFGetAttr(ncId, axis_id, CF_AXIS, a_sig));

            // If valid signify axis correctly
            if (a_sig == "X")
            {
                X = axis_id;
            }
            else if (a_sig == "Y")
            {
                Y = axis_id;
            }
            else if (a_sig == "Z")
            {
                Z = axis_id;
            }
            else
            {
                throw SG_Exception_Dep(container_name,
                                       "A node_coordinates variable", CF_AXIS);
            }

            this->touple_order++;
        }
        else
        {
            throw SG_Exception_Existential(container_name, dim);
        }

        dim = strtok(nullptr, " ");
    }

    // Write axis in X, Y, Z order

    if (X != INVALID_VAR_ID)
        this->nodec_varIds.push_back(X);
    else
    {
        throw SG_Exception_Existential(container_name,
                                       "node_coordinates: X-axis");
    }
    if (Y != INVALID_VAR_ID)
        this->nodec_varIds.push_back(Y);
    else
    {
        throw SG_Exception_Existential(container_name,
                                       "node_coordinates: Y-axis");
    }
    if (Z != INVALID_VAR_ID)
        this->nodec_varIds.push_back(Z);

    /* Final Checks for node coordinates
     * (1) Each axis has one and only one dimension, dim length of each node
     * coordinates are all equal (2) total node count == dim length of each node
     * coordinates (if node_counts not empty) (3) there are at least two node
     * coordinate variable ids
     */

    int all_dim = INVALID_VAR_ID;
    bool dim_set = false;
    //(1) one dimension check, each node_coordinates have same dimension
    for (size_t nvitr = 0; nvitr < nodec_varIds.size(); nvitr++)
    {
        int dimC = 0;
        nc_inq_varndims(ncId, nodec_varIds[nvitr], &dimC);

        if (dimC != 1)
        {
            throw SG_Exception_Not1D();
        }

        // check that all have same dimension
        int inter_dim[1];
        if (nc_inq_vardimid(ncId, nodec_varIds[nvitr], inter_dim) != NC_NOERR)
        {
            throw SG_Exception_Existential(
                container_name, "one or more node_coordinate dimensions");
        }

        if (!dim_set)
        {
            all_dim = inter_dim[0];
            dim_set = true;
        }

        else
        {
            if (inter_dim[0] != all_dim)
                throw SG_Exception_Dim_MM(
                    container_name, "X, Y",
                    "in general all node coordinate axes");
        }
    }

    // (2) check equality one
    if (node_counts.size() > 0)
    {
        size_t diml = 0;
        nc_inq_dimlen(ncId, all_dim, &diml);

        if (diml != total_node_count)
            throw SG_Exception_BadSum(container_name, "node_count",
                                      "node coordinate dimension length");
    }

    // (3) check tuple order
    if (this->touple_order < 2)
    {
        throw SG_Exception_Existential(
            container_name,
            "insufficient node coordinates must have at least two axis");
    }

    /* Investigate for instance dimension
     * The procedure is as follows
     *
     * (1) if there's node_count, use the dimension used to index node count
     * (2) otherwise it's point (singleton) data, in this case use the node
     * coordinate dimension
     */
    size_t instance_dim_len = 0;

    if (node_counts.size() >= 1)
    {
        int nc_dims = 0;
        nc_inq_varndims(ncId, nc_vid, &nc_dims);

        if (nc_dims != 1)
            throw SG_Exception_Not1D();

        int nc_dim_id[1];

        if (nc_inq_vardimid(ncId, nc_vid, nc_dim_id) != NC_NOERR)
        {
            throw SG_Exception_Existential(container_name,
                                           "node_count dimension");
        }

        this->inst_dimId = nc_dim_id[0];
    }

    else
    {
        this->inst_dimId = all_dim;
    }

    nc_inq_dimlen(ncId, this->inst_dimId, &instance_dim_len);

    if (instance_dim_len == 0)
        throw SG_Exception_EmptyDim();

    // Set values accordingly
    this->inst_dimLen = instance_dim_len;
    this->pt_buffer = std::make_unique<Point>(this->touple_order);
    this->gc_varId = geoVarId;
    this->ncid = ncId;
}

Point &SGeometry_Reader::operator[](size_t index)
{
    for (int order = 0; order < touple_order; order++)
    {
        Point &pt = *(this->pt_buffer);
        double data;
        size_t real_ind = index;

        // Read a single coord
        int err =
            nc_get_var1_double(ncid, nodec_varIds[order], &real_ind, &data);

        if (err != NC_NOERR)
        {
            throw SG_Exception_BadPoint();
        }

        pt[order] = data;
    }

    return *(this->pt_buffer);
}

size_t SGeometry_Reader::get_geometry_count()
{
    if (type == POINT)
    {
        if (this->nodec_varIds.size() < 1)
            return 0;

        // If more than one dim, then error. Otherwise inquire its length and
        // return that
        int dims;
        if (nc_inq_varndims(this->ncid, nodec_varIds[0], &dims) != NC_NOERR)
            return 0;
        if (dims != 1)
            return 0;

        // Find which dimension is used for x
        int index;
        if (nc_inq_vardimid(this->ncid, nodec_varIds[0], &index) != NC_NOERR)
        {
            return 0;
        }

        // Finally find the length
        size_t len;
        if (nc_inq_dimlen(this->ncid, index, &len) != NC_NOERR)
        {
            return 0;
        }
        return len;
    }

    else
        return this->node_counts.size();
}

static void add_to_buffer(std::vector<unsigned char> &buffer, uint8_t v)
{
    buffer.push_back(v);
}

static void add_to_buffer(std::vector<unsigned char> &buffer, uint32_t v)
{
    const size_t old_size = buffer.size();
    buffer.resize(old_size + sizeof(v));
    memcpy(&buffer[old_size], &v, sizeof(v));
}

static void add_to_buffer(std::vector<unsigned char> &, int) = delete;

static void add_to_buffer(std::vector<unsigned char> &buffer, double v)
{
    const size_t old_size = buffer.size();
    buffer.resize(old_size + sizeof(v));
    memcpy(&buffer[old_size], &v, sizeof(v));
}

/* serializeToWKB
 * Takes the geometry in SGeometry at a given index and converts it into WKB
 * format. Converting SGeometry into WKB automatically allocates the required
 * buffer space and returns a buffer that MUST be free'd
 */
std::vector<unsigned char> SGeometry_Reader::serializeToWKB(size_t featureInd)
{
    std::vector<unsigned char> ret;
    int nc = 0;
    size_t sb = 0;

    // Points don't have node_count entry... only inspect and set node_counts if
    // not a point
    if (this->getGeometryType() != POINT)
    {
        nc = node_counts[featureInd];
        sb = bound_list[featureInd];
    }

    // Serialization occurs differently depending on geometry
    // The memory requirements also differ between geometries
    switch (this->getGeometryType())
    {
        case POINT:
            inPlaceSerialize_Point(this, featureInd, ret);
            break;

        case LINE:
            inPlaceSerialize_LineString(this, nc, sb, ret);
            break;

        case POLYGON:
            // A polygon has:
            // 1 byte header
            // 4 byte Type
            // 4 byte ring count (1 (exterior)) [assume only exterior rings,
            // otherwise multipolygon] For each ring: 4 byte point count, 8 byte
            // double x point count x # dimension (This is equivalent to
            // requesting enough for all points and a point count header for
            // each point) (Or 8 byte double x node count x # dimension + 4 byte
            // point count x part_count)

            // if interior ring, then assume that it will be a multipolygon
            // (maybe future work?)
            inPlaceSerialize_PolygonExtOnly(this, nc, sb, ret);
            break;

        case MULTIPOINT:
        {
            uint8_t header = PLATFORM_HEADER;
            uint32_t t = this->touple_order == 2   ? wkbMultiPoint
                         : this->touple_order == 3 ? wkbMultiPoint25D
                                                   : wkbNone;

            if (t == wkbNone)
                throw SG_Exception_BadFeature();

            // Add metadata
            add_to_buffer(ret, header);
            add_to_buffer(ret, t);
            add_to_buffer(ret, static_cast<uint32_t>(nc));

            // Add points
            for (int pts = 0; pts < nc; pts++)
            {
                inPlaceSerialize_Point(this, static_cast<size_t>(sb + pts),
                                       ret);
            }
        }

        break;

        case MULTILINE:
        {
            uint8_t header = PLATFORM_HEADER;
            uint32_t t = this->touple_order == 2   ? wkbMultiLineString
                         : this->touple_order == 3 ? wkbMultiLineString25D
                                                   : wkbNone;

            if (t == wkbNone)
                throw SG_Exception_BadFeature();
            int32_t lc = parts_count[featureInd];
            size_t seek_begin = sb;
            size_t pc_begin =
                pnc_bl[featureInd];  // initialize with first part count, list
                                     // of part counts is contiguous
            std::vector<int> pnc;

            // Build sub vector for part_node_counts
            // + Calculate wkbSize
            for (int itr = 0; itr < lc; itr++)
            {
                pnc.push_back(pnode_counts[pc_begin + itr]);
            }

            size_t cur_point = seek_begin;
            size_t pcount = pnc.size();

            // Begin Writing
            add_to_buffer(ret, header);
            add_to_buffer(ret, t);
            add_to_buffer(ret, static_cast<uint32_t>(pcount));

            for (size_t itr = 0; itr < pcount; itr++)
            {
                inPlaceSerialize_LineString(this, pnc[itr], cur_point, ret);
                cur_point = pnc[itr] + cur_point;
            }
        }

        break;

        case MULTIPOLYGON:
        {
            uint8_t header = PLATFORM_HEADER;
            uint32_t t = this->touple_order == 2   ? wkbMultiPolygon
                         : this->touple_order == 3 ? wkbMultiPolygon25D
                                                   : wkbNone;

            if (t == wkbNone)
                throw SG_Exception_BadFeature();
            bool noInteriors = this->int_rings.size() == 0 ? true : false;
            int32_t rc = parts_count[featureInd];
            size_t seek_begin = sb;
            size_t pc_begin =
                pnc_bl[featureInd];  // initialize with first part count, list
                                     // of part counts is contiguous
            std::vector<int> pnc;

            // Build sub vector for part_node_counts
            for (int itr = 0; itr < rc; itr++)
            {
                pnc.push_back(pnode_counts[pc_begin + itr]);
            }

            // Create Multipolygon headers
            add_to_buffer(ret, header);
            add_to_buffer(ret, t);

            if (noInteriors)
            {
                size_t cur_point = seek_begin;
                size_t pcount = pnc.size();
                add_to_buffer(ret, static_cast<uint32_t>(pcount));

                for (size_t itr = 0; itr < pcount; itr++)
                {
                    inPlaceSerialize_PolygonExtOnly(this, pnc[itr], cur_point,
                                                    ret);
                    cur_point = pnc[itr] + cur_point;
                }
            }

            else
            {
                int32_t polys = poly_count[featureInd];
                add_to_buffer(ret, static_cast<uint32_t>(polys));

                size_t base = pnc_bl[featureInd];  // beginning of parts_count
                                                   // for this multigeometry
                size_t seek = seek_begin;  // beginning of node range for this
                                           // multigeometry
                size_t ir_base = base + 1;

                // has interior rings,
                for (int32_t itr = 0; itr < polys; itr++)
                {
                    int rc_m = 1;

                    // count how many parts belong to each Polygon
                    while (ir_base < int_rings.size() && int_rings[ir_base])
                    {
                        rc_m++;
                        ir_base++;
                    }

                    if (rc_m == 1)
                        ir_base++;  // single ring case

                    std::vector<int> poly_parts;

                    // Make part node count sub vector
                    for (int itr_2 = 0; itr_2 < rc_m; itr_2++)
                    {
                        poly_parts.push_back(pnode_counts[base + itr_2]);
                    }

                    inPlaceSerialize_Polygon(this, poly_parts, rc_m, seek, ret);

                    // Update seek position
                    for (size_t itr_3 = 0; itr_3 < poly_parts.size(); itr_3++)
                    {
                        seek += poly_parts[itr_3];
                    }

                    base += poly_parts.size();
                    // cppcheck-suppress redundantAssignment
                    ir_base = base + 1;
                }
            }
        }
        break;

        default:

            throw SG_Exception_BadFeature();
    }

    return ret;
}

void SGeometry_PropertyScanner::open(int container_id)
{
    // First check for container_id, if variable doesn't exist error out
    if (nc_inq_var(this->nc, container_id, nullptr, nullptr, nullptr, nullptr,
                   nullptr) != NC_NOERR)
    {
        return;  // change to exception
    }

    // Now exists, see what variables refer to this one
    // First get name of this container
    char contname[NC_MAX_NAME + 1];
    memset(contname, 0, NC_MAX_NAME + 1);
    if (nc_inq_varname(this->nc, container_id, contname) != NC_NOERR)
    {
        return;
    }

    // Then scan throughout the netcdfDataset if those variables
    // geometry_container attribute matches the container
    int varCount = 0;
    if (nc_inq_nvars(this->nc, &varCount) != NC_NOERR)
    {
        return;
    }

    std::string buf;
    for (int curr = 0; curr < varCount; curr++)
    {
        if (NCDFGetAttr(this->nc, curr, CF_SG_GEOMETRY, buf) != CE_None ||
            buf.empty())
        {
            continue;
        }

        // If matches, then establish a reference by placing this variable's
        // data in both vectors
        if (!strcmp(contname, buf.c_str()))
        {
            std::string property_name;

            // look for special OGR original name field
            if (NCDFGetAttr(this->nc, curr, OGR_SG_ORIGINAL_LAYERNAME,
                            property_name) != CE_None)
            {
                // if it doesn't exist, then just use the variable name
                property_name.resize(NC_MAX_NAME + 1, 0);

                if (nc_inq_varname(this->nc, curr, property_name.data()) !=
                    NC_NOERR)
                {
                    throw SG_Exception_General_Malformed(contname);
                }
            }

            v_ids.push_back(curr);
            v_headers.push_back(std::move(property_name));
        }
    }
}

// Exception Class Implementations
SG_Exception_Dim_MM::SG_Exception_Dim_MM(const char *container_name,
                                         const char *field_1,
                                         const char *field_2)
{
    std::string cn_s(container_name);
    std::string field1_s(field_1);
    std::string field2_s(field_2);

    this->err_msg = "[" + cn_s + "] One or more dimensions of " + field1_s +
                    " and " + field2_s + " do not match but must match.";
}

SG_Exception_Existential::SG_Exception_Existential(const char *container_name,
                                                   const char *missing_name)
{
    std::string cn_s(container_name);
    std::string mn_s(missing_name);

    this->err_msg = "[" + cn_s +
                    "] The property or the variable associated with " + mn_s +
                    " is missing.";
}

SG_Exception_Dep::SG_Exception_Dep(const char *container_name, const char *arg1,
                                   const char *arg2)
{
    std::string cn_s(container_name);
    std::string arg1_s(arg1);
    std::string arg2_s(arg2);

    this->err_msg = "[" + cn_s + "] The attribute " + arg1_s +
                    " may not exist without the attribute " + arg2_s +
                    " existing.";
}

SG_Exception_BadSum::SG_Exception_BadSum(const char *container_name,
                                         const char *arg1, const char *arg2)
{
    std::string cn_s(container_name);
    std::string arg1_s(arg1);
    std::string arg2_s(arg2);

    this->err_msg = "[" + cn_s + "]" + " The sum of all values in " + arg1_s +
                    " and " + arg2_s + " do not match.";
}

SG_Exception_General_Malformed ::SG_Exception_General_Malformed(const char *arg)
{
    std::string arg1_s(arg);

    this->err_msg =
        "Corruption or malformed formatting has been detected in: " + arg1_s;
}

// to get past linker
SG_Exception::~SG_Exception()
{
}

// Helpers
// following is a short hand for a clean up and exit, since goto isn't allowed
double getCFVersion(int ncid)
{
    double ver = -1.0;
    std::string attrVal;

    // Fetch the CF attribute
    if (NCDFGetAttr(ncid, NC_GLOBAL, NCDF_CONVENTIONS, attrVal) != CE_None ||
        attrVal.empty())
    {
        return ver;
    }

    if (sscanf(attrVal.c_str(), "CF-%lf", &ver) != 1)
    {
        return -1.0;
    }

    return ver;
}

geom_t getGeometryType(int ncid, int varid)
{
    geom_t ret = UNSUPPORTED;
    std::string gt_name_s;

    if (NCDFGetAttr(ncid, varid, CF_SG_GEOMETRY_TYPE, gt_name_s) != CE_None ||
        gt_name_s.empty())
    {
        return NONE;
    }

    // Points
    if (!strcmp(gt_name_s.c_str(), CF_SG_TYPE_POINT))
    {
        // Node Count not present? Assume that it is a multipoint.
        if (nc_inq_att(ncid, varid, CF_SG_NODE_COUNT, nullptr, nullptr) ==
            NC_ENOTATT)
        {
            ret = POINT;
        }
        else
            ret = MULTIPOINT;
    }

    // Lines
    else if (!strcmp(gt_name_s.c_str(), CF_SG_TYPE_LINE))
    {
        // Part Node Count present? Assume multiline
        if (nc_inq_att(ncid, varid, CF_SG_PART_NODE_COUNT, nullptr, nullptr) ==
            NC_ENOTATT)
        {
            ret = LINE;
        }
        else
            ret = MULTILINE;
    }

    // Polygons
    else if (!strcmp(gt_name_s.c_str(), CF_SG_TYPE_POLY))
    {
        /* Polygons versus MultiPolygons, slightly ambiguous
         * Part Node Count & no Interior Ring - MultiPolygon
         * no Part Node Count & no Interior Ring - Polygon
         * Part Node Count & Interior Ring - assume that it is a MultiPolygon
         */
        int pnc_present =
            nc_inq_att(ncid, varid, CF_SG_PART_NODE_COUNT, nullptr, nullptr);
        int ir_present =
            nc_inq_att(ncid, varid, CF_SG_INTERIOR_RING, nullptr, nullptr);

        if (pnc_present == NC_ENOTATT && ir_present == NC_ENOTATT)
        {
            ret = POLYGON;
        }
        else
            ret = MULTIPOLYGON;
    }

    return ret;
}

void inPlaceSerialize_Point(SGeometry_Reader *ge, size_t seek_pos,
                            std::vector<unsigned char> &buffer)
{
    uint8_t order = PLATFORM_HEADER;
    uint32_t t = ge->get_axisCount() == 2   ? wkbPoint
                 : ge->get_axisCount() == 3 ? wkbPoint25D
                                            : wkbNone;

    if (t == wkbNone)
        throw SG_Exception_BadFeature();

    add_to_buffer(buffer, order);
    add_to_buffer(buffer, t);

    // Now get point data;
    Point &p = (*ge)[seek_pos];
    add_to_buffer(buffer, p[0]);
    add_to_buffer(buffer, p[1]);

    if (ge->get_axisCount() >= 3)
    {
        add_to_buffer(buffer, p[2]);
    }
}

void inPlaceSerialize_LineString(SGeometry_Reader *ge, int node_count,
                                 size_t seek_begin,
                                 std::vector<unsigned char> &buffer)
{
    uint8_t order = PLATFORM_HEADER;
    uint32_t t = ge->get_axisCount() == 2   ? wkbLineString
                 : ge->get_axisCount() == 3 ? wkbLineString25D
                                            : wkbNone;

    if (t == wkbNone)
        throw SG_Exception_BadFeature();
    uint32_t nc = (uint32_t)node_count;

    add_to_buffer(buffer, order);
    add_to_buffer(buffer, t);
    add_to_buffer(buffer, nc);

    // Now serialize points
    for (int ind = 0; ind < node_count; ind++)
    {
        Point &p = (*ge)[seek_begin + ind];
        add_to_buffer(buffer, p[0]);
        add_to_buffer(buffer, p[1]);

        if (ge->get_axisCount() >= 3)
        {
            add_to_buffer(buffer, p[2]);
        }
    }
}

void inPlaceSerialize_PolygonExtOnly(SGeometry_Reader *ge, int node_count,
                                     size_t seek_begin,
                                     std::vector<unsigned char> &buffer)
{
    uint8_t order = PLATFORM_HEADER;
    uint32_t t = ge->get_axisCount() == 2   ? wkbPolygon
                 : ge->get_axisCount() == 3 ? wkbPolygon25D
                                            : wkbNone;

    if (t == wkbNone)
        throw SG_Exception_BadFeature();
    uint32_t rc = 1;

    add_to_buffer(buffer, order);
    add_to_buffer(buffer, t);
    add_to_buffer(buffer, rc);
    add_to_buffer(buffer, static_cast<uint32_t>(node_count));
    for (int pind = 0; pind < node_count; pind++)
    {
        Point &pt = (*ge)[seek_begin + pind];
        add_to_buffer(buffer, pt[0]);
        add_to_buffer(buffer, pt[1]);

        if (ge->get_axisCount() >= 3)
        {
            add_to_buffer(buffer, pt[2]);
        }
    }
}

void inPlaceSerialize_Polygon(SGeometry_Reader *ge, std::vector<int> &pnc,
                              int ring_count, size_t seek_begin,
                              std::vector<unsigned char> &buffer)
{
    uint8_t order = PLATFORM_HEADER;
    uint32_t t = ge->get_axisCount() == 2   ? wkbPolygon
                 : ge->get_axisCount() == 3 ? wkbPolygon25D
                                            : wkbNone;

    if (t == wkbNone)
        throw SG_Exception_BadFeature();
    uint32_t rc = static_cast<uint32_t>(ring_count);

    add_to_buffer(buffer, order);
    add_to_buffer(buffer, t);
    add_to_buffer(buffer, rc);
    int cmoffset = 0;
    for (int ring_c = 0; ring_c < ring_count; ring_c++)
    {
        uint32_t node_count = pnc[ring_c];
        add_to_buffer(buffer, node_count);

        int pind = 0;
        for (pind = 0; pind < pnc[ring_c]; pind++)
        {
            Point &pt = (*ge)[seek_begin + cmoffset + pind];
            add_to_buffer(buffer, pt[0]);
            add_to_buffer(buffer, pt[1]);

            if (ge->get_axisCount() >= 3)
            {
                add_to_buffer(buffer, pt[2]);
            }
        }

        cmoffset += pind;
    }
}

int scanForGeometryContainers(int ncid, std::set<int> &r_ids)
{
    int nvars;
    if (nc_inq_nvars(ncid, &nvars) != NC_NOERR)
    {
        return -1;
    }

    r_ids.clear();

    // For each variable check for a geometry attribute
    // If it has geometry attribute, then check the associated variable ID
    std::string buf;
    for (int itr = 0; itr < nvars; itr++)
    {
        if (NCDFGetAttr(ncid, itr, CF_SG_GEOMETRY, buf) != CE_None)
        {
            continue;
        }

        int varID;
        if (nc_inq_varid(ncid, buf.c_str(), &varID) != NC_NOERR)
        {
            continue;
        }

        // Now have variable ID. See if vector contains it, and if not
        // insert
        r_ids.insert(varID);  // It's a set. No big deal sets don't allow
                              // duplicates anyways
    }

    return 0;
}
}  // namespace nccfdriver

Coverage Report

Created: 2026-08-14 09:29