ncdas.cc

// -*- mode: c++; c-basic-offset:4 -*-
 
// This file is part of nc_handler, a data handler for the OPeNDAP data
// server.
 
// Copyright (c) 2002,2003 OPeNDAP, Inc.
// Author: James Gallagher <jgallagher@opendap.org>
//
// This is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License as published by the Free
// Software Foundation; either version 2.1 of the License, or (at your
// option) any later version.
//
// This software is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
// License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
//
// You can contact OPeNDAP, Inc. at PO Box 112, Saunderstown, RI. 02874-0112.
 
// (c) COPYRIGHT URI/MIT 1994-1996
// Please read the full copyright statement in the file COPYRIGHT.
//
// Authors:
//      reza            Reza Nekovei (reza@intcomm.net)
 
// This file contains functions which read the variables and their attributes
// from a netcdf file and build the in-memory DAS. These functions form the
// core of the server-side software necessary to extract the DAS from a
// netcdf data file.
//
// It also contains test code which will print the in-memory DAS to
// stdout. It uses both the DAS class and the netcdf library.
// In addition, parts of these functions were taken from the netcdf program
// ncdump, from the netcdf standard distribution (ver 2.3.2)
//
// jhrg 9/23/94
 
#include "config_nc.h"
 
#include <iostream>
#include <string>
#include <sstream>
#include <iomanip>
#include <vector>
 
#include <cmath>
 
#include <netcdf.h>
 
#include <libdap/util.h>
#include <libdap/escaping.h>
#include <libdap/DAS.h>
 
#include <BESDebug.h>
 
#include "NCRequestHandler.h"
#include "nc_util.h"
 
#define ATTR_STRING_QUOTE_FIX 1
#define STOP_ESCAPING_STRING_ATTRS 0    // Never do this. jhrg 2/28/22
 
#define MODULE "nc"
#define prolog std::string("ncdas::").append(__func__).append("() - ")
 
using namespace libdap;

static string print_attr(nc_type type, int loc, void *vals)
{
    ostringstream rep;
    union {
        char *cp;
        char **stringp;
        int16_t *sp;
        uint16_t *usp;
        int32_t *i;
        uint32_t *ui;
        float *fp;
        double *dp;
    } gp;
 
    switch (type) {
    case NC_UBYTE:
        unsigned char uc;
        gp.cp = (char *) vals;
 
        uc = *(gp.cp + loc);
        rep << (int) uc;
        return rep.str();
 
    case NC_BYTE:
        if (NCRequestHandler::get_promote_byte_to_short()) {
            signed char sc;
            gp.cp = (char *) vals;
 
            sc = *(gp.cp + loc);
            rep << (int) sc;
            return rep.str();
        }
        else {
            unsigned char uc;
            gp.cp = (char *) vals;
 
            uc = *(gp.cp + loc);
            rep << (int) uc;
            return rep.str();
        }
 
    case NC_CHAR:
    {
        // Don't escape the special characters. 
        // Special characters will be handled in libdap4.
        // KY 2022-08-25
#if 0
        return escattr(static_cast<const char*>(vals));
#endif
        string tmp_str = static_cast<const char*>(vals);
        return tmp_str;
    }
 
    case NC_STRING:
        gp.stringp = (char **) vals;
        rep << *(gp.stringp + loc);
        return rep.str();
 
    case NC_SHORT:
        gp.sp = (short *) vals;
        rep << *(gp.sp + loc);
        return rep.str();
 
    case NC_USHORT:
        gp.usp = (uint16_t *) vals;
        rep << *(gp.usp + loc);
        return rep.str();
 
    case NC_INT:
        gp.i = (int32_t *) vals; // warning: long int format, int arg (arg 3)
        rep << *(gp.i + loc);
        return rep.str();
 
    case NC_UINT:
        gp.ui = (uint32_t *) vals;
        rep << *(gp.ui + loc);
        return rep.str();
 
    case NC_FLOAT: {
        gp.fp = (float *) vals;
        float valAtLoc = *(gp.fp + loc);
 
        rep << std::showpoint;
        rep << std::setprecision(9);
 
        if (isnan(valAtLoc)) {
            rep << "NaN";
        }
        else {
            rep << valAtLoc;
        }
        // If there's no decimal point and the rep does not use scientific
        // notation, add a decimal point. This little jaunt was taken because
        // this code is modeled after older code and that's what it did. I'm
        // trying to keep the same behavior as the old code without its
        // problems. jhrg 8/11/2006
        string tmp_value = rep.str();
        if (tmp_value.find('.') == string::npos && tmp_value.find('e') == string::npos && tmp_value.find('E') == string::npos
            && tmp_value.find("nan") == string::npos && tmp_value.find("NaN") == string::npos && tmp_value.find("NAN") == string::npos) rep << ".";
        return rep.str();
    }
 
    case NC_DOUBLE: {
        gp.dp = (double *) vals;
        double valAtLoc = *(gp.dp + loc);
 
        rep << std::showpoint;
        rep << std::setprecision(16);
 
        if (std::isnan(valAtLoc)) {
            rep << "NaN";
        }
        else {
            rep << valAtLoc;
        }
        string tmp_value = rep.str();
        if (tmp_value.find('.') == string::npos && tmp_value.find('e') == string::npos && tmp_value.find('E') == string::npos
            && tmp_value.find("nan") == string::npos && tmp_value.find("NaN") == string::npos && tmp_value.find("NAN") == string::npos) rep << ".";
        return rep.str();
    }
 
    default:
        if (NCRequestHandler::get_ignore_unknown_types())
            cerr << "The netcdf handler tried to print an attribute that has an unrecognized type. (1)" << endl;
        else
            throw InternalErr(__FILE__, __LINE__, "The netcdf handler tried to print an attribute that has an unrecognized type. (1)");
        break;
    }
 
    return "";
}

static string print_type(nc_type datatype)
{
    switch (datatype) {
    case NC_STRING:
    case NC_CHAR:
        return "String";
 
    case NC_UBYTE:
        return "Byte";
 
    case NC_BYTE:
        if (NCRequestHandler::get_promote_byte_to_short()) {
            return "Int16";
        }
        else {
            return "Byte";
        }
 
    case NC_SHORT:
        return "Int16";
 
    case NC_INT:
        return "Int32";
 
    case NC_USHORT:
        return "UInt16";
 
    case NC_UINT:
        return "UInt32";
 
    case NC_FLOAT:
        return "Float32";
 
    case NC_DOUBLE:
        return "Float64";
 
    case NC_COMPOUND:
        return "NC_COMPOUND";
 
        // These are all new netcdf 4 types that we don't support yet
        // as attributes. It's useful to have a print representation for
        // them so that we can return useful information about why some
        // information was elided or an exception thrown.
    case NC_INT64:
        return "NC_INT64";
 
    case NC_UINT64:
        return "NC_UINT64";
 
    case NC_VLEN:
        return "NC_VLEN";
    case NC_OPAQUE:
        return "NC_OPAQUE";
    case NC_ENUM:
        return "NC_ENUM";
 
    default:
        if (NCRequestHandler::get_ignore_unknown_types())
            cerr << "The netcdf handler tried to print an attribute that has an unrecognized type. (2)" << endl;
        else
            throw InternalErr(__FILE__, __LINE__, "The netcdf handler tried to print an attribute that has an unrecognized type. (2)");
        break;
    }
 
    return "";
}

static void append_values(int ncid, int v, int len, nc_type datatype, char *attrname, AttrTable *at)
{
    size_t size;
    int errstat = nc_inq_type(ncid, datatype, 0, &size);
    if (errstat != NC_NOERR) throw Error(errstat, "Could not get the size for the type.");
 
    vector<char> value((len + 1) * size);
    errstat = nc_get_att(ncid, v, attrname, value.data());
    if (errstat != NC_NOERR) {
        throw Error(errstat, string("Could not get the value for attribute '") + attrname + string("'"));
    }
 
    // If the datatype is NC_CHAR then we have a string. netCDF 3
    // represents strings as arrays of char, but we represent them as X
    // strings. So... Add the null and set the length to 1
    if (datatype == NC_CHAR) {
        value[len] = '\0';
        len = 1;
    }
 
    // add all the attributes in the array
    for (int loc = 0; loc < len; loc++) {
        string print_rep = print_attr(datatype, loc, value.data());
        at->append_attr(attrname, print_type(datatype), print_rep);
    }
}

static void read_attributes_netcdf4(int ncid, int varid, int natts, AttrTable *at)
{
    BESDEBUG(MODULE, prolog << "In read_attributes_netcdf4" << endl);
 
    for (int attr_num = 0; attr_num < natts; ++attr_num) {
        int errstat = NC_NOERR;
        // Get the attribute name
        char attrname[MAX_NC_NAME];
        errstat = nc_inq_attname(ncid, varid, attr_num, attrname);
        if (errstat != NC_NOERR) throw Error(errstat, "Could not get the name for attribute " + long_to_string(attr_num));
 
        // Get datatype and len; len is the number of values.
        nc_type datatype;
        size_t len;
        errstat = nc_inq_att(ncid, varid, attrname, &datatype, &len);
        if (errstat != NC_NOERR) throw Error(errstat, "Could not get the name for attribute '" + string(attrname) + "'");
 
        BESDEBUG(MODULE, prolog << "nc_inq_att returned datatype = " << datatype << " for '" << attrname << "'" << endl);
 
        if (datatype >= NC_FIRSTUSERTYPEID) {
            char type_name[NC_MAX_NAME + 1];
            size_t size;
            nc_type base_type;
            size_t nfields;
            int class_type;
            errstat = nc_inq_user_type(ncid, datatype, type_name, &size, &base_type, &nfields, &class_type);
            if (errstat != NC_NOERR)
                throw(InternalErr(__FILE__, __LINE__, "Could not get information about a user-defined type (" + long_to_string(errstat) + ")."));
 
            BESDEBUG(MODULE, prolog << "Before switch(class_type)" << endl);
            switch (class_type) {
            case NC_COMPOUND: {
                // Make recursive attrs work?
                vector<unsigned char> values((len + 1) * size);
 
                int errstat = nc_get_att(ncid, varid, attrname, values.data());
                if (errstat != NC_NOERR) throw Error(errstat, string("Could not get the value for attribute '") + attrname + string("'"));
 
                for (size_t i = 0; i < nfields; ++i) {
                    char field_name[NC_MAX_NAME + 1];
                    nc_type field_typeid;
                    size_t field_offset;
                    nc_inq_compound_field(ncid, datatype, i, field_name, &field_offset, &field_typeid, 0, 0);
 
                    at->append_attr(field_name, print_type(field_typeid), print_attr(field_typeid, 0, values.data() + field_offset));
                }
                break;
            }
 
            case NC_VLEN:
                if (NCRequestHandler::get_ignore_unknown_types())
                    cerr << "in build_user_defined; found a vlen." << endl;
                else
                    throw Error("The netCDF handler does not yet support the NC_VLEN type.");
                break;
 
            case NC_OPAQUE: {
                vector<unsigned char> values((len + 1) * size);
 
                int errstat = nc_get_att(ncid, varid, attrname, values.data());
                if (errstat != NC_NOERR) throw Error(errstat, string("Could not get the value for attribute '") + attrname + string("'"));
 
                for (size_t i = 0; i < size; ++i)
                    at->append_attr(attrname, print_type(NC_BYTE), print_attr(NC_BYTE, i, values.data()));
 
                break;
            }
 
            case NC_ENUM: {
#if 0
                nc_type basetype;
                size_t base_size, num_members;
                errstat = nc_inq_enum(ncid, datatype, 0/*char *name*/, &basetype, &base_size, &num_members);
                if (errstat != NC_NOERR) throw Error(errstat, string("Could not get the size of the enum base type for '") + attrname + string("'"));
#endif
                vector<unsigned char> values((len + 1) * size);
 
                int errstat = nc_get_att(ncid, varid, attrname, values.data());
                if (errstat != NC_NOERR) throw Error(errstat, string("Could not get the value for attribute '") + attrname + string("'"));
 
                for (size_t i = 0; i < len; ++i)
                    at->append_attr(attrname, print_type(base_type), print_attr(base_type, i, values.data()));
 
                break;
            }
 
            default:
                throw InternalErr(__FILE__, __LINE__, "Expected one of NC_COMPOUND, NC_VLEN, NC_OPAQUE or NC_ENUM");
            }
 
            BESDEBUG(MODULE, prolog << "After switch(class-type)" << endl);
        }
        else {
            switch (datatype) {
            case NC_STRING:
            case NC_BYTE:
            case NC_CHAR:
            case NC_SHORT:
            case NC_INT:
            case NC_FLOAT:
            case NC_DOUBLE:
            case NC_UBYTE:
            case NC_USHORT:
            case NC_UINT:
                BESDEBUG(MODULE, prolog << "Before append_values ..." << endl);
                append_values(ncid, varid, len, datatype, attrname, at);
                BESDEBUG(MODULE, prolog << "After append_values ..." << endl);
                break;
 
            case NC_INT64:
            case NC_UINT64: {
                string note = "Attribute elided: Unsupported attribute type ";
                note += "(" + print_type(datatype) + ")";
                at->append_attr(attrname, "String", note);
                break;
            }
 
            case NC_COMPOUND:
            case NC_VLEN:
            case NC_OPAQUE:
            case NC_ENUM:
                throw InternalErr(__FILE__, __LINE__, "user-defined attribute type not recognized as such!");
 
            default:
                throw InternalErr(__FILE__, __LINE__, "Unrecognized attribute type.");
            }
        }
    }
    BESDEBUG(MODULE, prolog << "Exiting read_attributes_netcdf4" << endl);
}

void nc_read_dataset_attributes(DAS &das, const string &filename)
{
    BESDEBUG(MODULE, prolog << "In nc_read_dataset_attributes" << endl);
 
    int ncid, errstat;
    errstat = nc_open(filename.c_str(), NC_NOWRITE, &ncid);
    if (errstat != NC_NOERR) throw Error(errstat, "NetCDF handler: Could not open " + filename + ".");
 
    // how many variables? how many global attributes?
    int nvars, ngatts;
    errstat = nc_inq(ncid, (int *) 0, &nvars, &ngatts, (int *) 0);
    if (errstat != NC_NOERR) throw Error(errstat, "NetCDF handler: Could not inquire about netcdf file: " + path_to_filename(filename) + ".");
 
    // for each variable
    char varname[MAX_NC_NAME];
    int natts = 0;
    nc_type var_type;
    for (int varid = 0; varid < nvars; ++varid) {
        BESDEBUG(MODULE, prolog << "Top of for loop; for each var..." << endl);
 
        errstat = nc_inq_var(ncid, varid, varname, &var_type, (int*) 0, (int*) 0, &natts);
        if (errstat != NC_NOERR) throw Error(errstat, "Could not get information for variable: " + long_to_string(varid));
 
        AttrTable *attr_table_ptr = das.get_table(varname);
        if (!attr_table_ptr) attr_table_ptr = das.add_table(varname, new AttrTable);
 
        read_attributes_netcdf4(ncid, varid, natts, attr_table_ptr);
 
        // Add a special attribute for string lengths
        if (var_type == NC_CHAR) {
            // number of dimensions and size of Nth dimension
            int num_dim;
            int vdimids[MAX_VAR_DIMS]; // variable dimension ids
            errstat = nc_inq_var(ncid, varid, (char *) 0, (nc_type *) 0, &num_dim, vdimids, (int *) 0);
            if (errstat != NC_NOERR)
                throw Error(errstat, string("NetCDF handler: Could not read information about a NC_CHAR variable while building the DAS."));
 
            if (num_dim == 0) {
                // a scalar NC_CHAR is stuffed into a string of length 1
                int size = 1;
                string print_rep = print_attr(NC_INT, 0, (void *) &size);
                attr_table_ptr->append_attr("string_length", print_type(NC_INT), print_rep);
            }
            else {
                // size_t *dim_sizes = new size_t[num_dim];
                vector<size_t> dim_sizes(num_dim);
                for (int i = 0; i < num_dim; ++i) {
                    if ((errstat = nc_inq_dimlen(ncid, vdimids[i], &dim_sizes[i])) != NC_NOERR) {
                        throw Error(errstat,
                            string("NetCDF handler: Could not read dimension information about the variable `") + varname + string("'."));
                    }
                }
 
                // add attribute
                string print_rep = print_attr(NC_INT, 0, (void *) (&dim_sizes[num_dim - 1]));
                attr_table_ptr->append_attr("string_length", print_type(NC_INT), print_rep);
            }
        }
        else if (is_user_defined_type(ncid, var_type)) {
            //var_type >= NC_FIRSTUSERTYPEID) {
            vector<char> name(MAX_NC_NAME + 1);
            int class_type;
            errstat = nc_inq_user_type(ncid, var_type, name.data(), 0, 0, 0, &class_type);
            if (errstat != NC_NOERR)
                throw(InternalErr(__FILE__, __LINE__, "Could not get information about a user-defined type (" + long_to_string(errstat) + ")."));
 
            switch (class_type) {
            case NC_OPAQUE: {
                attr_table_ptr->append_attr("DAP2_OriginalNetCDFBaseType", print_type(NC_STRING), "NC_OPAQUE");
                attr_table_ptr->append_attr("DAP2_OriginalNetCDFTypeName", print_type(NC_STRING), name.data());
                break;
            }
 
            case NC_ENUM: {
                //vector<char> name(MAX_NC_NAME + 1);
                nc_type base_nc_type;
                size_t base_size, num_members;
                errstat = nc_inq_enum(ncid, var_type, 0/*name.data()*/, &base_nc_type, &base_size, &num_members);
                if (errstat != NC_NOERR)
                    throw(InternalErr(__FILE__, __LINE__, "Could not get information about an enum(" + long_to_string(errstat) + ")."));
 
                // If the base type is a 64-bit int, bail with an error or
                // a message about unsupported types
                if (base_nc_type == NC_INT64 || base_nc_type == NC_UINT64) {
                    if (NCRequestHandler::get_ignore_unknown_types())
                        cerr << "An Enum uses 64-bit integers, but this handler does not support that type." << endl;
                    else
                        throw Error("An Enum uses 64-bit integers, but this handler does not support that type.");
                    break;
                }
 
                for (size_t i = 0; i < num_members; ++i) {
                    vector<char> member_name(MAX_NC_NAME + 1);
                    vector<char> member_value(base_size);
                    errstat = nc_inq_enum_member(ncid, var_type, i, member_name.data(), member_value.data());
                    if (errstat != NC_NOERR)
                        throw(InternalErr(__FILE__, __LINE__, "Could not get information about an enum value (" + long_to_string(errstat) + ")."));
                    attr_table_ptr->append_attr("DAP2_EnumValues", print_type(base_nc_type), print_attr(base_nc_type, 0, member_value.data()));
                    attr_table_ptr->append_attr("DAP2_EnumNames", print_type(NC_STRING), member_name.data());
                }
 
                attr_table_ptr->append_attr("DAP2_OriginalNetCDFBaseType", print_type(NC_STRING), "NC_ENUM");
                attr_table_ptr->append_attr("DAP2_OriginalNetCDFTypeName", print_type(NC_STRING), name.data());
 
                break;
            }
 
            default:
                break;
            }
        }
    }
 
    BESDEBUG(MODULE, prolog << "Starting global attributes" << endl);
 
    // global attributes
    if (ngatts > 0) {
        AttrTable *attr_table_ptr = das.add_table("NC_GLOBAL", new AttrTable);
        read_attributes_netcdf4(ncid, NC_GLOBAL, ngatts, attr_table_ptr);
    }
 
    // Add unlimited dimension name in DODS_EXTRA attribute table
    int xdimid;
    char dimname[MAX_NC_NAME];
    nc_type datatype = NC_CHAR;
    if ((errstat = nc_inq(ncid, (int *) 0, (int *) 0, (int *) 0, &xdimid)) != NC_NOERR)
        throw InternalErr(__FILE__, __LINE__, string("NetCDF handler: Could not access variable information: ") + nc_strerror(errstat));
    if (xdimid != -1) {
        if ((errstat = nc_inq_dim(ncid, xdimid, dimname, (size_t *) 0)) != NC_NOERR)
            throw InternalErr(__FILE__, __LINE__, string("NetCDF handler: Could not access dimension information: ") + nc_strerror(errstat));
        string print_rep = print_attr(datatype, 0, dimname);
        AttrTable *attr_table_ptr = das.add_table("DODS_EXTRA", new AttrTable);
        attr_table_ptr->append_attr("Unlimited_Dimension", print_type(datatype), print_rep);
    }
 
    if (nc_close(ncid) != NC_NOERR) throw InternalErr(__FILE__, __LINE__, "NetCDF handler: Could not close the dataset!");
 
    BESDEBUG(MODULE, prolog << "Exiting nc_read_dataset_attributes" << endl);
}