nc_file.h

Go to the documentation of this file.

#pragma once
 
#include <filesystem>
#include <functional>
#include <list>
#include "../dg/blas.h"
#include "../dg/backend/memory.h"
#include "nc_error.h"
#include "easy_atts.h"
#include "easy_output.h"
#include "easy_input.h"
 
namespace dg
{
namespace file
{
enum NcFileMode
{
    nc_nowrite,  
    nc_write,    
    nc_clobber, 
    nc_noclobber, 
};
 
struct SerialNcFile
{
    using Hyperslab = NcHyperslab;
    // ///////////////////////////// CONSTRUCTORS/DESTRUCTOR /////////
    SerialNcFile () = default;
    SerialNcFile(const std::filesystem::path& filename,
            enum NcFileMode mode = nc_nowrite)
    {
        open( filename, mode);
    }
    SerialNcFile(const SerialNcFile& rhs) = delete;
    SerialNcFile& operator =(const SerialNcFile & rhs) = delete;
 
    SerialNcFile(SerialNcFile&& rhs) = default;
    SerialNcFile& operator =(SerialNcFile && rhs) = default;
 
    ~SerialNcFile()
    {
        try
        {
            close();
        }
        catch (NC_Error &e)
        {
            std::cerr << e.what() << std::endl;
        }
    }
    // /////////////////// open/close /////////
 
    void open(const std::filesystem::path& filename,
            enum NcFileMode mode = nc_nowrite)
    {
        // Like a std::fstream opening fails if file already associated
        if( m_open)
            throw NC_Error( 1002);
 
        // TODO Test the pathing on Windows
        NC_Error_Handle err;
        switch (mode)
        {
            case nc_nowrite:
                err = nc_open( filename.string().c_str(), NC_NETCDF4 |
                        NC_NOWRITE, &m_ncid);
                break;
            case nc_write:
                err = nc_open( filename.string().c_str(), NC_NETCDF4 |
                        NC_WRITE, &m_ncid);
                break;
            case nc_noclobber:
                err = nc_create( filename.string().c_str(), NC_NETCDF4 |
                        NC_NOCLOBBER, &m_ncid);
                break;
            case nc_clobber:
                err = nc_create( filename.string().c_str(), NC_NETCDF4 |
                        NC_CLOBBER, &m_ncid);
                break;
        }
        m_open = true;
        m_grp = m_ncid;
    }
    bool is_open() const noexcept
    {
        // is_closed() == not is_open()
        return m_open;
    }
 
    void close()
    {
        // when the file is closed the ncid may be assigned to a new file
        if (m_open)
        {
            NC_Error_Handle err;
            err = nc_close(m_ncid);
        }
 
        m_open = false;
        m_grp = m_ncid = 0;
    }
 
    void sync()
    {
        check_open();
        NC_Error_Handle err;
        err = nc_sync( m_ncid);
    }
    int get_ncid() const noexcept{ return m_ncid;}
 
    // ///////////// Groups /////////////////
    void def_grp( std::string name)
    {
        check_open();
        int new_grp_id = 0;
        NC_Error_Handle err;
        err = nc_def_grp( m_grp, name.c_str(), &new_grp_id);
        // we just forget the id => always need to ask netcdf for id
        // Is no performance hit as cached
    }
 
    void def_grp_p( std::filesystem::path path)
    {
        check_open();
        if( not path.has_root_path()) // it is a relative path
        {
            auto current = get_current_path();
            path = (current / path);
        }
        int groupid = m_ncid;
        auto rel_path = path.relative_path();
        NC_Error_Handle err;
        for( auto it = rel_path.begin(); it != rel_path.end(); it++)
        {
            std::string grp = *it;
            int new_grpid;
            int retval = nc_inq_ncid( groupid, grp.c_str(), &new_grpid);
            if( retval != NC_NOERR)
                err = nc_def_grp( groupid, grp.c_str(), &groupid);
            else
                groupid = new_grpid;
        }
    }
    bool grp_is_defined( std::filesystem::path path) const
    {
        check_open();
        std::string name = path.generic_string();
        if( not path.has_root_path()) // it is a relative path
        {
            auto current = get_current_path();
            name = (current / path ).generic_string();
        }
        int grpid=0;
        int retval = nc_inq_grp_full_ncid( m_ncid, name.c_str(), &grpid);
        return retval == NC_NOERR;
    }
    void set_grp( std::filesystem::path path = "")
    {
        check_open();
        NC_Error_Handle err;
        std::string name = path.generic_string();
        if ( name == ".")
            return;
        if ( name == "" or name == "/")
        {
            m_grp = m_ncid;
        }
        else if ( name == "..")
        {
            if( m_grp == m_ncid)
                return;
            err = nc_inq_grp_parent( m_grp, &m_grp);
        }
        else
        {
            if( not path.has_root_path()) // it is a relative path
            {
                auto current = get_current_path();
                name = (current / path ).generic_string();
            }
            err = nc_inq_grp_full_ncid( m_ncid, name.c_str(), &m_grp);
        }
    }
    void rename_grp( std::string old_name, std::string new_name)
    {
        check_open();
        NC_Error_Handle err;
        int old_grp;
        err = nc_inq_grp_ncid( m_grp, old_name.c_str(), &old_grp);
        err = nc_rename_grp( old_grp, new_name.c_str());
    }
 
    int get_grpid() const noexcept{ return m_grp;}
 
    std::filesystem::path get_current_path( ) const
    {
        check_open();
        return get_grp_path( m_grp);
    }
 
    std::list<std::filesystem::path> get_grps( ) const
    {
        check_open();
        auto grps = get_grps_abs(m_grp);
        std::list<std::filesystem::path> grps_v;
        for( auto grp : grps)
            grps_v.push_back( grp.second);
        return grps_v;
 
    }
    std::list<std::filesystem::path> get_grps_r( ) const
    {
        auto grps = get_grps_abs_r(m_grp);
        std::list<std::filesystem::path> grps_v;
        for( auto grp : grps)
            grps_v.push_back( grp.second);
        return grps_v;
    }
 
    // //////////// Dimensions ////////////////////////
    void def_dim( std::string name, size_t size)
    {
        check_open();
        NC_Error_Handle err;
        int dim_id;
        err = nc_def_dim( m_grp, name.c_str(), size, &dim_id);
    }
    void rename_dim( std::string old_name, std::string new_name)
    {
        check_open();
        int dimid;
        NC_Error_Handle err;
        err = nc_inq_dimid( m_grp, old_name.c_str(), &dimid);
        err = nc_rename_dim( m_grp, dimid, new_name.c_str());
    }
    size_t get_dim_size( std::string name) const
    {
        check_open();
        NC_Error_Handle err;
        int dimid;
        err = nc_inq_dimid( m_grp, name.c_str(), &dimid);
        size_t len;
        err = nc_inq_dimlen( m_grp, dimid, &len);
        return len;
 
    }
 
    std::vector<size_t> get_dims_shape( const std::vector<std::string>& dims) const
    {
        std::vector<size_t> shape( dims.size());
        for( unsigned u=0; u<dims.size(); u++)
            shape[u] = get_dim_size( dims[u]);
        return shape;
    }
    std::vector<std::string> get_dims( bool include_parents = true) const
    {
        check_open();
        NC_Error_Handle err;
        // There is a correct way of getting the dimensions in a group
        // https://github.com/Unidata/netcdf-c/issues/2873
        int ndims;
        err = nc_inq_dimids( m_grp, &ndims, NULL, include_parents);
        if( ndims == 0)
            return {};
        int dimids[ndims];
        err = nc_inq_dimids( m_grp, &ndims, dimids, include_parents);
        // Globally dimension ids are 0, 1, 2, ... in the order in which the
        // dimensions were defined
        std::vector<std::string> dims;
        for ( int i = 0; i < ndims; i++)
        {
            char dimname [ NC_MAX_NAME+1];
            err = nc_inq_dimname( m_grp, dimids[i], dimname);
            if( std::find(dims.begin(), dims.end(), std::string(dimname) ) ==
                    dims.end())
                dims.push_back(dimname);
        }
        return dims;
    }
 
    std::vector<std::string> get_unlim_dims( ) const
    {
        check_open();
        int ndims;
        NC_Error_Handle err;
        err = nc_inq_unlimdims( m_grp, &ndims, NULL);
        if( ndims == 0)
            return {};
        int dimids[ndims];
        // Our tests indicate that this does not return the unlimited dimensions
        // of the parent group even though the documentation says so...
        err = nc_inq_unlimdims( m_grp, &ndims, dimids);
        std::vector<std::string> dims;
        for( int i=0; i<ndims; i++)
        {
            char dimname [ NC_MAX_NAME+1];
            err = nc_inq_dimname( m_grp, dimids[i], dimname);
            if( std::find(dims.begin(), dims.end(), std::string(dimname) ) ==
                    dims.end())
                dims.push_back(dimname);
        }
        return dims;
    }
 
    bool dim_is_defined( std::string name) const
    {
        int dimid=0;
        int retval = nc_inq_dimid( m_grp, name.c_str(), &dimid);
        return retval == NC_NOERR;
    }
 
    void put_att( const std::pair<std::string, nc_att_t>& att, std::string id = "")
    {
        int varid = name2varid( id);
        dg::file::detail::put_att( m_grp, varid, att);
    }
 
    template<class S, class T> // T cannot be nc_att_t
    void put_att( const std::tuple<S,nc_type, T>& att, std::string id = "")
    {
        int varid = name2varid( id);
        dg::file::detail::put_att( m_grp, varid, att);
    }
    template<class Attributes = std::map<std::string, nc_att_t> > // *it must be usable in put_att
    void put_atts(const Attributes& atts, std::string id = "")
    {
        int varid = name2varid( id);
        dg::file::detail::put_atts( m_grp, varid, atts);
    }
 
    // ///////////////// Attribute getters
 
    template<class T>
    T get_att_as(std::string att_name, std::string id = "") const
    {
        int varid = name2varid( id);
        return dg::file::detail::get_att_as<T>( m_grp, varid, att_name);
    }
 
    template<class T>
    std::vector<T> get_att_vec_as(std::string att_name, std::string id = "") const
    {
        return get_att_as<std::vector<T>>( att_name, id);
    }
 
    template<class T>
    std::map<std::string, T> get_atts_as( std::string id = "") const
    {
        int varid = name2varid( id);
        return dg::file::detail::get_atts_as<T>( m_grp, varid);
    }
 
    std::map<std::string, nc_att_t> get_atts( std::string id = "") const
    {
        return get_atts_as<nc_att_t>( id);
    }
 
    void del_att(std::string att_name, std::string id = "")
    {
        int varid = name2varid( id);
        auto name = att_name.c_str();
        NC_Error_Handle err;
        err = nc_del_att( m_grp, varid, name);
    }
    bool att_is_defined(std::string att_name, std::string id = "") const
    {
        int varid = name2varid( id);
        int attid;
        int retval = nc_inq_attid( m_grp, varid, att_name.c_str(), &attid);
        return retval == NC_NOERR;
    }
    void rename_att(std::string old_att_name, std::string
            new_att_name, std::string id = "")
    {
        int varid = name2varid( id);
        auto old_name = old_att_name.c_str();
        auto new_name = new_att_name.c_str();
        NC_Error_Handle err;
        err = nc_rename_att( m_grp, varid, old_name, new_name);
    }
 
    // //////////// Variables ////////////////////////
    template<class T, class Attributes = std::map<std::string, nc_att_t>>
    void def_var_as( std::string name,
        const std::vector<std::string>& dim_names,
        const Attributes& atts = {})
    {
        def_var( name, detail::getNCDataType<T>(), dim_names, atts);
    }
 
    template<class Attributes = std::map<std::string, nc_att_t>>
    void def_var( std::string name, nc_type xtype,
            const std::vector<std::string>& dim_names,
            const Attributes& atts = {})
    {
        file::NC_Error_Handle err;
        std::vector<int> dimids( dim_names.size());
        for( unsigned u=0; u<dim_names.size(); u++)
            err = nc_inq_dimid( m_grp, dim_names[u].c_str(), &dimids[u]);
        int varid;
        err = nc_def_var( m_grp, name.c_str(), xtype, dim_names.size(),
                &dimids[0], &varid);
        put_atts( atts, name);
    }
 
    template<class ContainerType, std::enable_if_t< dg::is_vector_v<
        ContainerType, SharedVectorTag>, bool > = true>
    void put_var( std::string name, const NcHyperslab& slab,
            const ContainerType& data)
    {
        int varid = name2varid( name);
        file::NC_Error_Handle err;
        int ndims;
        err = nc_inq_varndims( m_grp, varid, &ndims);
        assert( (unsigned)ndims == slab.ndim());
        if constexpr ( dg::has_policy_v<ContainerType, dg::CudaTag>)
        {
            using value_type = dg::get_value_type<ContainerType>;
            m_buffer.template set<value_type>( data.size());
            auto& buffer = m_buffer.template get<value_type>( );
            dg::assign ( data, buffer);
            err = detail::put_vara_T( m_grp, varid, slab.startp(), slab.countp(),
                buffer.data());
        }
        else
            err = detail::put_vara_T( m_grp, varid, slab.startp(), slab.countp(),
                thrust::raw_pointer_cast(data.data()));
    }
 
    template<class ContainerType, class Attributes = std::map<std::string, nc_att_t>,
        std::enable_if_t< dg::is_vector_v<ContainerType, SharedVectorTag>, bool > = true>
    void defput_var( std::string name, const std::vector<std::string>& dim_names,
            const Attributes& atts, const NcHyperslab& slab,
            const ContainerType& data)
    {
        def_var_as<dg::get_value_type<ContainerType>>( name, dim_names, atts);
        put_var( name, slab, data);
    }
 
    template<class T, std::enable_if_t< dg::is_scalar_v<T>, bool> = true>
    void put_var( std::string name, const std::vector<size_t>& start, T data)
    {
        int varid = name2varid( name);
        file::NC_Error_Handle err;
        std::vector<size_t> count( start.size(), 1);
        err = detail::put_vara_T( m_grp, varid, &start[0], &count[0], &data);
    }
 
    template<class T, class Attributes = std::map<std::string, nc_att_t>>
    void def_dimvar_as( std::string name, size_t size, const Attributes& atts)
    {
        def_dim( name, size);
        def_var_as<T>( name, {name}, atts);
    }
 
    template<class ContainerType, class Attributes = std::map<std::string, nc_att_t>>
    void defput_dim( std::string name,
            const Attributes& atts,
            const ContainerType& abscissas)
    {
        def_dimvar_as<dg::get_value_type<ContainerType>>( name,
            abscissas.size(), atts);
        put_var( name, {abscissas}, abscissas);
    }
 
    template<class ContainerType, std::enable_if_t< dg::is_vector_v<
        ContainerType, SharedVectorTag>, bool > = true>
    void get_var( std::string name, const NcHyperslab& slab,
            ContainerType& data) const
    {
        int varid = name2varid( name);
        file::NC_Error_Handle err;
        int ndims;
        err = nc_inq_varndims( m_grp, varid, &ndims);
        assert( (unsigned)ndims == slab.ndim());
        unsigned size = 1;
        for( unsigned u=0; u<slab.ndim(); u++)
            size *= slab.count()[u];
        if constexpr ( dg::has_policy_v<ContainerType, dg::CudaTag>)
        {
            using value_type = dg::get_value_type<ContainerType>;
            m_buffer.template set<value_type>( size);
            auto& buffer = m_buffer.template get<value_type>( );
            err = detail::get_vara_T( m_grp, varid, slab.startp(), slab.countp(),
                buffer.data());
            dg::assign ( buffer, data);
        }
        else
        {
            data.resize( size);
            err = detail::get_vara_T( m_grp, varid, slab.startp(), slab.countp(),
                thrust::raw_pointer_cast(data.data()));
        }
    }
 
    template<class T, std::enable_if_t< dg::is_scalar_v<T>, bool> = true>
    void get_var( std::string name, const std::vector<size_t>& start, T& data) const
    {
        int varid = name2varid( name);
        file::NC_Error_Handle err;
        int ndims;
        err = nc_inq_varndims( m_grp, varid, &ndims);
        assert( (unsigned)ndims == start.size());
        if( ndims == 0)
            err = detail::get_vara_T( m_grp, varid, NULL, NULL, &data);
        else
        {
            std::vector<size_t> count( start.size(), 1);
            err = detail::get_vara_T( m_grp, varid, &start[0], &count[0], &data);
        }
    }
 
    bool var_is_defined( std::string name) const
    {
        check_open();
        int varid=0;
        int retval = nc_inq_varid( m_grp, name.c_str(), &varid);
        return retval == NC_NOERR;
    }
 
    nc_type get_var_type(std::string name) const
    {
        int varid = name2varid( name);
        int xtype;
        file::NC_Error_Handle err;
        err = nc_inq_vartype( m_grp, varid, &xtype);
        return xtype;
    }
    std::vector<std::string> get_var_dims(std::string name) const
    {
        int varid = name2varid( name);
        file::NC_Error_Handle err;
 
        int ndims;
        err = nc_inq_varndims( m_grp, varid, &ndims);
        if( ndims == 0)
            return {};
        int dimids[ndims];
        err = nc_inq_vardimid( m_grp, varid, dimids);
 
        std::vector<std::string> dims(ndims);
        for( int i=0; i<ndims; i++)
        {
            char dimname [ NC_MAX_NAME+1];
            err = nc_inq_dimname( m_grp, dimids[i], dimname);
            dims[i] = dimname;
        }
        return dims;
    }
 
    std::list<std::string> get_var_names() const
    {
        check_open();
        return get_var_names_private( m_grp);
    }
 
    private:
    int name2varid( std::string id) const
    {
        // Variable ids are persistent once created
        // Attribute ids can change if one deletes attributes
        // This is fast even for lots variables (1000 vars take <ms)
        check_open();
        NC_Error_Handle err;
        int varid;
        if ( id == "")
            varid = NC_GLOBAL;
        else
        {
            err = nc_inq_varid( m_grp, id.c_str(), &varid);
        }
        return varid;
 
    }
    // Absolute current path
    std::filesystem::path get_grp_path( int ncid ) const
    {
        size_t len;
        NC_Error_Handle err;
        err = nc_inq_grpname_full( ncid, &len, NULL);
        std::string current( len, 'x');
        err = nc_inq_grpname_full( ncid, &len, &current[0]);
        return current;
    }
    // Absolute paths of subgroups
    std::map<int, std::filesystem::path> get_grps_abs( int ncid ) const
    {
        NC_Error_Handle err;
        int num_grps;
        err = nc_inq_grps( ncid, &num_grps, NULL);
        if( num_grps == 0)
            return {};
        std::vector<int> group_ids( num_grps);
        err = nc_inq_grps( ncid, &num_grps, &group_ids[0]);
        std::map<int, std::filesystem::path> groups;
        for( int i=0; i<num_grps; i++)
        {
            size_t len;
            err = nc_inq_grpname_full( group_ids[i], &len, NULL);
            std::string name( len, 'z');
            err = nc_inq_grpname_full( group_ids[i], &len, &name[0]);
            groups[group_ids[i]] = name;
        }
        return groups;
    }
    std::map<int, std::filesystem::path> get_grps_abs_r( int ncid) const
    {
        auto grps = get_grps_abs(ncid);
        for( auto grp : grps)
        {
            auto subgrps = get_grps_abs_r( grp.first);
            grps.merge( subgrps);
        }
        return grps;
    }
    // Need a different name so mpi_file can invoke
    std::list<std::string> get_var_names_private(int grpid) const
    {
        std::list<std::string> vars;
        int num_vars = 0, num_dims;
        file::NC_Error_Handle err;
        err = nc_inq(grpid, &num_dims, &num_vars, NULL, NULL);
        // https://docs.unidata.ucar.edu/netcdf-c/current/reading_unknown.html
        for( int i=0; i<num_vars; i++)
        {
            char name[NC_MAX_NAME+1]; // 256
            err = nc_inq_varname( grpid, i, name);
            vars.push_back( name);
        }
        return vars;
    }
 
 
    void check_open() const
    {
        if( !m_open)
            throw NC_Error( 1000);
    }
 
    bool m_open = false;
    int m_ncid = 0; // ncid can change by closing and re-opening a file
    int m_grp = 0;
    // the currently active group (All group ids in all open files are unique
    // and thus group ids can be different by opening the same file twice),
    // dims can be seen by all child groups
 
    // Buffer for device to host transfer, and dg::assign
    mutable dg::detail::AnyVector<thrust::host_vector> m_buffer;
    // ABANDONED: Variable tracker (persists on closing and opening a different
    // file) The problem with trying to track is
    // 1) do we track every file that is opened?
    // 2) we cannot prevent someone from opening closing a file here and simply
    // opening the same file somewhere else. The info is lost or corrupted then
    //std::map<std::filesystem::path, std::map<std::filesystem::path, NcVarInfo>>
    //    m_vars;
};
 
#ifndef MPI_VERSION
using NcFile = SerialNcFile;
#endif // MPI_VERSION
 
}// namespace file
}// namespace dg