WriteNC.cpp

Go to the documentation of this file.

#include "WriteNC.hpp"
#include "moab/CN.hpp"
#include "MBTagConventions.hpp"
#include "MBParallelConventions.h"
#include "moab/Interface.hpp"
#include "moab/Range.hpp"
#include "moab/WriteUtilIface.hpp"
#include "moab/FileOptions.hpp"
#include "NCWriteHelper.hpp"
 
#include <fstream>
#include <map>
#include <set>
 
#include <iostream>
#include <sstream>
 
#ifdef WIN32
#ifdef size_t
#undef size_t
#endif
#endif
 
namespace moab
{
 
WriterIface* WriteNC::factory( Interface* iface )
{
 return new WriteNC( iface );
}
 
WriteNC::WriteNC( Interface* impl )
 : mbImpl( impl ), dbgOut( stderr ),
#ifdef MOAB_HAVE_MPI
 myPcomm( NULL ),
#endif
 noMesh( false ), noVars( false ), append( false ), mGlobalIdTag( 0 ), isParallel( false ), myHelper( NULL )
{
 assert( impl != NULL );
 impl->query_interface( mWriteIface );
}
 
WriteNC::~WriteNC()
{
 mbImpl->release_interface( mWriteIface );
 if( myHelper != NULL ) delete myHelper;
}
 
//! Writes out a file
ErrorCode WriteNC::write_file( const char* file_name,
 const bool overwrite,
 const FileOptions& options,
 const EntityHandle* file_set,
 const int num_set,
 const std::vector< std::string >&,
 const Tag*,
 int,
 int )
{
 ErrorCode rval;
 // See if opts has variable(s) specified
 std::vector< std::string > var_names;
 std::vector< std::string > desired_names;
 std::vector< int > tstep_nums;
 std::vector< double > tstep_vals;
 
 // Get and cache predefined tag handles
 mGlobalIdTag = mbImpl->globalId_tag();
 
 // num set has to be 1, we will write only one set, the original file set used to load
 if( num_set != 1 ) MB_SET_ERR( MB_FAILURE, "We should write only one set (the file set used to read data into)" );
 
 rval = parse_options( options, var_names, desired_names, tstep_nums, tstep_vals );MB_CHK_SET_ERR( rval, "Trouble parsing option string" );
 
 // Important to create some data that will be used to write the file; dimensions, variables, etc
 // new variables still need to have some way of defining their dimensions
 // maybe it will be passed as write options
 rval = process_conventional_tags( *file_set );MB_CHK_SET_ERR( rval, "Trouble processing conventional tags" );
 
 // Create or append the file
 if( append )
 dbgOut.tprintf( 1, "opening file %s for appending \n", file_name );
 else
 dbgOut.tprintf( 1, "creating file %s\n", file_name );
 fileName = file_name;
 int success;
 
 if( append )
 {
 int omode = NC_WRITE;
#ifdef MOAB_HAVE_PNETCDF
 if( isParallel )
 success = NCFUNC( open )( myPcomm->proc_config().proc_comm(), file_name, omode, MPI_INFO_NULL, &fileId );
 else
 success = NCFUNC( open )( MPI_COMM_SELF, file_name, omode, MPI_INFO_NULL, &fileId );
#else
 // This is a regular netcdf file, open in write mode
 success = NCFUNC( open )( file_name, omode, &fileId );
#endif
 if( success ) MB_SET_ERR( MB_FAILURE, "Trouble opening file " << file_name << " for appending" );
 }
 else
 { // Case when the file is new, will be overwritten, most likely
 int cmode = overwrite ? NC_CLOBBER : NC_NOCLOBBER;
#ifdef MOAB_HAVE_PNETCDF
 if( isParallel )
 success = NCFUNC( create )( myPcomm->proc_config().proc_comm(), file_name, cmode, MPI_INFO_NULL, &fileId );
 else
 success = NCFUNC( create )( MPI_COMM_SELF, file_name, cmode, MPI_INFO_NULL, &fileId );
#else
 // This is a regular netcdf file
 success = NCFUNC( create )( file_name, cmode, &fileId );
#endif
 if( success ) MB_SET_ERR( MB_FAILURE, "Trouble creating file " << file_name << " for writing" );
 }
 
 if( NULL != myHelper ) delete myHelper;
 
 // Get appropriate helper instance for WriteNC class based on some info in the file set
 myHelper = NCWriteHelper::get_nc_helper( this, fileId, options, *file_set );
 if( NULL == myHelper )
 {
 MB_SET_ERR( MB_FAILURE, "Failed to get NCWriteHelper class instance" );
 }
 
 rval = myHelper->collect_mesh_info();MB_CHK_SET_ERR( rval, "Trouble collecting mesh information" );
 
 rval = myHelper->collect_variable_data( var_names, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble collecting variable data" );
 
 rval = myHelper->init_file( var_names, desired_names, append );MB_CHK_SET_ERR( rval, "Trouble initializing file" );
 
 rval = myHelper->write_values( var_names, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble writing values to file" );
 
 success = NCFUNC( close )( fileId );
 if( success ) MB_SET_ERR( MB_FAILURE, "Trouble closing file" );
 
 return MB_SUCCESS;
}
 
ErrorCode WriteNC::parse_options( const FileOptions& opts,
 std::vector< std::string >& var_names,
 std::vector< std::string >& desired_names,
 std::vector< int >& tstep_nums,
 std::vector< double >& tstep_vals )
{
 int tmpval;
 if( MB_SUCCESS == opts.get_int_option( "DEBUG_IO", 1, tmpval ) )
 {
 dbgOut.set_verbosity( tmpval );
 dbgOut.set_prefix( "NCWrite" );
 }
 
 ErrorCode rval = opts.get_strs_option( "VARIABLE", var_names );
 if( MB_TYPE_OUT_OF_RANGE == rval )
 noVars = true;
 else
 noVars = false;
 
 rval = opts.get_strs_option( "RENAME", desired_names );
 if( MB_ENTITY_NOT_FOUND == rval )
 {
 if( !noVars )
 {
 desired_names.resize( var_names.size() );
 std::copy( var_names.begin(), var_names.end(), desired_names.begin() );
 }
 }
 // Either way
 assert( desired_names.size() == var_names.size() );
 
 opts.get_ints_option( "TIMESTEP", tstep_nums );
 opts.get_reals_option( "TIMEVAL", tstep_vals );
 rval = opts.get_null_option( "NOMESH" );
 if( MB_SUCCESS == rval ) noMesh = true;
 
 rval = opts.get_null_option( "APPEND" );
 if( MB_SUCCESS == rval ) append = true;
 
 if( 2 <= dbgOut.get_verbosity() )
 {
 if( !var_names.empty() )
 {
 std::cerr << "Variables requested: ";
 for( unsigned int i = 0; i < var_names.size(); i++ )
 std::cerr << var_names[i];
 std::cerr << std::endl;
 }
 if( !tstep_nums.empty() )
 {
 std::cerr << "Timesteps requested: ";
 for( unsigned int i = 0; i < tstep_nums.size(); i++ )
 std::cerr << tstep_nums[i];
 std::cerr << std::endl;
 }
 if( !tstep_vals.empty() )
 {
 std::cerr << "Time vals requested: ";
 for( unsigned int i = 0; i < tstep_vals.size(); i++ )
 std::cerr << tstep_vals[i];
 std::cerr << std::endl;
 }
 }
 
// FIXME: copied from ReadNC, may need revise
#ifdef MOAB_HAVE_MPI
 isParallel = ( opts.match_option( "PARALLEL", "WRITE_PART" ) != MB_ENTITY_NOT_FOUND );
 
 if( !isParallel )
 // Return success here, since rval still has _NOT_FOUND from not finding option
 // in this case, myPcomm will be NULL, so it can never be used; always check for isParallel
 // before any use for myPcomm
 return MB_SUCCESS;
 
 int pcomm_no = 0;
 rval = opts.get_int_option( "PARALLEL_COMM", pcomm_no );
 if( MB_TYPE_OUT_OF_RANGE == rval )
 {
 MB_SET_ERR( rval, "Invalid value for PARALLEL_COMM option" );
 }
 
 myPcomm = ParallelComm::get_pcomm( mbImpl, pcomm_no );
 if( 0 == myPcomm )
 {
 myPcomm = new ParallelComm( mbImpl, MPI_COMM_WORLD );
 }
 
#ifndef MOAB_HAVE_PNETCDF
 const int procs = myPcomm->proc_config().proc_size();
 if( procs > 1 )
 {
 MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Attempt to launch NC writer in parallel without pnetcdf support" );
 }
#endif
 
 const int rank = myPcomm->proc_config().proc_rank();
 dbgOut.set_rank( rank );
#endif
 
 return MB_SUCCESS;
}
 
// This is the inverse process to create conventional tags
// Will look at <pargal_source>/src/core/fileinfo.cpp, init dim, vars, atts
ErrorCode WriteNC::process_conventional_tags( EntityHandle fileSet )
{
 ErrorCode rval;
 
 // Start copy
 Tag dimNamesTag = 0;
 std::string tag_name = "__DIM_NAMES";
 const void* data = NULL;
 int dimNamesSz = 0;
 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, dimNamesTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
 rval = mbImpl->tag_get_by_ptr( dimNamesTag, &fileSet, 1, &data, &dimNamesSz );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
 const char* p = static_cast< const char* >( data );
 dbgOut.tprintf( 1, "__DIM_NAMES tag has string length %d\n", dimNamesSz );
 
 std::size_t start = 0;
 
 Tag dimLensTag = 0;
 tag_name = "__DIM_LENS";
 data = NULL;
 int dimLensSz = 0;
 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_INTEGER, dimLensTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
 rval = mbImpl->tag_get_by_ptr( dimLensTag, &fileSet, 1, &data, &dimLensSz );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
 const int* int_p = static_cast< const int* >( data );
 dbgOut.tprintf( 1, "__DIM_LENS tag has %d values\n", dimLensSz );
 
 int idxDim = 0;
 // Dim names are separated by '\0' in the string of __DIM_NAMES tag
 for( std::size_t i = 0; i != static_cast< std::size_t >( dimNamesSz ); i++ )
 {
 if( p[i] == '\0' )
 {
 std::string dim_name( &p[start], i - start );
 int len = int_p[idxDim];
 dimNames.push_back( dim_name );
 dimLens.push_back( len );
 dbgOut.tprintf( 2, "Dimension %s has length %d\n", dim_name.c_str(), len );
 // FIXME: Need info from moab to set unlimited dimension
 // Currently assume each file has the same number of time dimensions
 /*if ((dim_name == "time") || (dim_name == "Time"))
 insert(dim_name, *(new pcdim(dim_name, len * m_file_names.size(), true)));
 else
 insert(dim_name, *(new pcdim(dim_name, len)));*/
 start = i + 1;
 idxDim++;
 }
 }
 
 Tag meshTypeTag = 0;
 tag_name = "__MESH_TYPE";
 data = NULL;
 int meshTypeSz = 0;
 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, meshTypeTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
 rval = mbImpl->tag_get_by_ptr( meshTypeTag, &fileSet, 1, &data, &meshTypeSz );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
 p = static_cast< const char* >( data );
 grid_type = std::string( &p[0], meshTypeSz );
 dbgOut.tprintf( 2, "Mesh type: %s\n", grid_type.c_str() );
 
 // Read <__VAR_NAMES_LOCATIONS> tag
 Tag varNamesLocsTag = 0;
 tag_name = "__VAR_NAMES_LOCATIONS";
 data = NULL;
 int varNamesLocsSz = 0;
 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_INTEGER, varNamesLocsTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
 rval = mbImpl->tag_get_by_ptr( varNamesLocsTag, &fileSet, 1, &data, &varNamesLocsSz );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
 int_p = static_cast< const int* >( data );
 std::vector< int > varNamesLocs( varNamesLocsSz );
 std::copy( int_p, int_p + varNamesLocsSz, varNamesLocs.begin() );
 
 Tag varNamesTag = 0;
 tag_name = "__VAR_NAMES";
 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, varNamesTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
 data = NULL;
 int varNamesSz = 0;
 rval = mbImpl->tag_get_by_ptr( varNamesTag, &fileSet, 1, &data, &varNamesSz );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
 dbgOut.tprintf( 2, "__VAR_NAMES tag has string length %d\n", varNamesSz );
 p = static_cast< const char* >( data );
 
 start = 0;
 int idxVar = 0;
 int sz;
 // Var names are separated by '\0' in the string of __VAR_NAMES tag
 for( std::size_t i = 0; i != static_cast< std::size_t >( varNamesSz ); i++ )
 {
 if( p[i] == '\0' )
 {
 std::string var_name( &p[start], i - start );
 
 dbgOut.tprintf( 2, "var name: %s index %d \n", var_name.c_str(), idxVar );
 // Process var name:
 // This will create/initiate map; we will populate variableDataStruct with info about
 // dims, tags, etc reference & is important; otherwise variableDataStruct will go out of
 // scope, and deleted :(
 VarData& variableDataStruct = varInfo[var_name];
 variableDataStruct.varName = var_name;
 variableDataStruct.entLoc = varNamesLocs[idxVar];
 
 dbgOut.tprintf( 2, "at var name %s varInfo size %d \n", var_name.c_str(), (int)varInfo.size() );
 
 sz = 0;
 Tag dims_tag = 0;
 std::string dim_names = "__" + var_name + "_DIMS";
 rval = mbImpl->tag_get_handle( dim_names.c_str(), 0, MB_TYPE_OPAQUE, dims_tag, MB_TAG_ANY );
 if( MB_SUCCESS != rval )
 {
 if( MB_TAG_NOT_FOUND == rval )
 {
 dbgOut.tprintf( 2, "tag : %s not found, continue \n", dim_names.c_str() );
 start = i + 1;
 idxVar++;
 continue;
 }
 MB_SET_ERR( rval, "Trouble getting conventional tag " << dim_names );
 }
 rval = mbImpl->tag_get_length( dims_tag, sz );MB_CHK_SET_ERR( rval, "Trouble getting size of dimensions for variable " << var_name );
 sz /= sizeof( Tag ); // The type is MB_TYPE_OPAQUE, but it is a list of tags, so we
 // need to divide by the size of Tag
 // sz is used for number of dimension tags in this list
 dbgOut.tprintf( 2, "var name: %s has %d dimensions \n", var_name.c_str(), sz );
 
 variableDataStruct.varDims.resize( sz );
 const void* ptr = NULL;
 rval = mbImpl->tag_get_by_ptr( dims_tag, &fileSet, 1, &ptr );
 
 const Tag* ptags = static_cast< const moab::Tag* >( ptr );
 for( std::size_t j = 0; j != static_cast< std::size_t >( sz ); j++ )
 {
 std::string dim_name;
 rval = mbImpl->tag_get_name( ptags[j], dim_name );MB_CHK_SET_ERR( rval, "Trouble getting dimension of variable " << var_name );
 dbgOut.tprintf( 2, "var name: %s has %s as dimension \n", var_name.c_str(), dim_name.c_str() );
 std::vector< std::string >::iterator vit = std::find( dimNames.begin(), dimNames.end(), dim_name );
 if( vit == dimNames.end() )
 MB_SET_ERR( MB_FAILURE, "Dimension " << dim_name << " not found for variable " << var_name );
 variableDataStruct.varDims[j] = (int)( vit - dimNames.begin() ); // Will be used for writing
 // This will have to change to actual file dimension, for writing
 }
 
 // Attributes for this variable
 std::stringstream ssTagName;
 ssTagName << "__" << var_name << "_ATTRIBS";
 tag_name = ssTagName.str();
 Tag varAttTag = 0;
 rval =
 mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, varAttTag, MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
 const void* varAttPtr = NULL;
 int varAttSz = 0;
 rval = mbImpl->tag_get_by_ptr( varAttTag, &fileSet, 1, &varAttPtr, &varAttSz );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
 if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Tag retrieved for variable %s\n", tag_name.c_str() );
 
 std::string attribString( (char*)varAttPtr, (char*)varAttPtr + varAttSz );
 if( attribString == "NO_ATTRIBS" )
 {
 // This variable has no attributes
 variableDataStruct.numAtts = 0;
 }
 else if( attribString == "DUMMY_VAR" )
 {
 // This variable is a dummy coordinate variable
 variableDataStruct.numAtts = 0;
 dummyVarNames.insert( variableDataStruct.varName );
 }
 else
 {
 ssTagName << "_LEN";
 tag_name = ssTagName.str();
 Tag varAttLenTag = 0;
 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_INTEGER, varAttLenTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
 int varAttLenSz = 0;
 rval = mbImpl->tag_get_length( varAttLenTag, varAttLenSz );MB_CHK_SET_ERR( rval, "Trouble getting length of conventional tag " << tag_name );
 std::vector< int > varAttLen( varAttLenSz );
 rval = mbImpl->tag_get_data( varAttLenTag, &fileSet, 1, &varAttLen[0] );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
 
 rval = process_concatenated_attribute( varAttPtr, varAttSz, varAttLen, variableDataStruct.varAtts );MB_CHK_SET_ERR( rval, "Trouble processing attributes of variable " << var_name );
 
 if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Tag metadata for variable %s\n", tag_name.c_str() );
 }
 // End attribute
 
 start = i + 1;
 idxVar++;
 } // if (p[i] == '\0')
 }
 
 // Global attributes
 tag_name = "__GLOBAL_ATTRIBS";
 Tag globalAttTag = 0;
 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, globalAttTag, MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
 std::vector< int > gattLen;
 
 const void* gattptr = NULL;
 int globalAttSz = 0;
 rval = mbImpl->tag_get_by_ptr( globalAttTag, &fileSet, 1, &gattptr, &globalAttSz );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
 
 if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Tag value retrieved for %s size %d\n", tag_name.c_str(), globalAttSz );
 
 // <__GLOBAL_ATTRIBS_LEN>
 tag_name = "__GLOBAL_ATTRIBS_LEN";
 Tag globalAttLenTag = 0;
 
 rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_INTEGER, globalAttLenTag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting conventional tag " << tag_name );
 int sizeGAtt = 0;
 rval = mbImpl->tag_get_length( globalAttLenTag, sizeGAtt );MB_CHK_SET_ERR( rval, "Trouble getting length of conventional tag " << tag_name );
 gattLen.resize( sizeGAtt );
 rval = mbImpl->tag_get_data( globalAttLenTag, &fileSet, 1, &gattLen[0] );MB_CHK_SET_ERR( rval, "Trouble getting data of conventional tag " << tag_name );
 if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Tag retrieved for variable %s\n", tag_name.c_str() );
 
 rval = process_concatenated_attribute( gattptr, globalAttSz, gattLen, globalAtts );MB_CHK_SET_ERR( rval, "Trouble processing global attributes" );
 
 return MB_SUCCESS;
}
 
// Reverse process from create_attrib_string
ErrorCode WriteNC::process_concatenated_attribute( const void* attPtr,
 int attSz,
 std::vector< int >& attLen,
 std::map< std::string, AttData >& attributes )
{
 std::size_t start = 0;
 std::size_t att_counter = 0;
 std::string concatString( (char*)attPtr, (char*)attPtr + attSz );
 
 for( std::size_t i = 0; i != (size_t)attSz; i++ )
 {
 if( concatString[i] == '\0' )
 {
 std::string att_name( &concatString[start], i - start );
 start = i + 1;
 while( concatString[i] != ';' )
 ++i;
 std::string data_type( &concatString[start], i - start );
 ++i;
 start = i;
 i = attLen[att_counter];
 if( concatString[i] != ';' ) MB_SET_ERR( MB_FAILURE, "Error parsing attributes" );
 
 std::string data_val( &concatString[start], i - start );
 start = i + 1;
 
 AttData& attrib = attributes[att_name];
 attrib.attValue = data_val;
 attrib.attLen = data_val.size();
 
 if( data_type == "char" )
 attrib.attDataType = NC_CHAR;
 else if( data_type == "double" )
 attrib.attDataType = NC_DOUBLE;
 else if( data_type == "float" )
 attrib.attDataType = NC_FLOAT;
 else if( data_type == "int" )
 attrib.attDataType = NC_INT;
 else if( data_type == "short" )
 attrib.attDataType = NC_SHORT;
 
 ++att_counter;
 dbgOut.tprintf( 2, " Process attribute %s with value %s \n", att_name.c_str(), data_val.c_str() );
 }
 }
 
 return MB_SUCCESS;
}
 
} // namespace moab