C++ NcVar::get_att方法代码示例

ERMsg C20thReanalysisProject::ReadData( CString filePath, CString varName, CData3Array& data)
{
	//typedef boost::multi_array<int, 2> array_type;


	ERMsg msg;
	
	NcError::set_err( NcError::silent_nonfatal );

  
	//CString filePath = m_path + "cccma_cgcm3_1-20c3m-run1-pr-1961-2000_monthly.nc";//GetFilePath(v, year, m);
	NcFile file(filePath);	//current year file

	if( !file.is_valid() )
	{
		CString err;
		err.FormatMessage(IDS_CMN_UNABLE_OPEN_READ, filePath);
		msg.ajoute(err);
		return msg;
	}

	NcVar* pVarData = file.get_var((LPCTSTR)varName);//is the varaible always at ffirst???
	
	//CString varName = pVarData->name();
	
	size_t sizeTime = pVarData->get_dim(0)->size();
	size_t sizeY = pVarData->get_dim(1)->size();
	size_t sizeX = pVarData->get_dim(2)->size();
	float offset = pVarData->get_att("add_offset")->as_float(0);
	float scaleFactor = pVarData->get_att("scale_factor")->as_float(0);
	
	boost::multi_array<short, 3> tmp(boost::extents[sizeTime][sizeY][sizeX]);

	ENSURE( pVarData->num_dims() == 3);
	
	if( pVarData->get(&(tmp[0][0][0]), sizeTime, sizeY, sizeX) )
	{
		//tmp.extens()
		//data.resize(sizeTime);
		//data.resize(boost::extents[sizeTime][sizeY][sizeX]);
		//apply offset and scale factor
		for(size_t i=0; i<tmp.size(); i++)
			for(size_t j=0; j<tmp[i].size(); j++)
				for(size_t k=0; k<tmp[i][j].size(); k++)
					data[i][j][k] = tmp[i][j][k]*scaleFactor+offset;

		file.close();
	}
	else
	{
		msg.ajoute( "Unable to get NetCDFData");
	}

	return msg;
}

void dumpatts(NcVar& var)
{
    NcToken vname = var.name();
    NcAtt* ap;
    for(int n = 0; ap = var.get_att(n); n++) {
	cout << "\t\t" << vname << ":" << ap->name() << " = " ;
	NcValues* vals = ap->values();
	cout << *vals << " ;" << endl ;
	delete ap;
	delete vals;
    }
}

QMap<QString, QString> DataInterfaceNetCdfVector::metaStrings(const QString& field)
{
  QMap<QString, QString> fieldStrings;
  QString tmpString;
  NcVar *var = netcdf._ncfile->get_var(field.toLatin1().constData());
  for (int i=0; i<var->num_atts(); ++i) {
    NcAtt *att = var->get_att(i);
    // Only handle char/unspecified attributes as fieldStrings, the others as fieldScalars
    if (att->type() == NC_CHAR || att->type() == NC_UNSPECIFIED) {
      fieldStrings[att->name()] = QString(att->values()->as_string(0));
    }
    // qDebug() << att->name() << ": " << att->values()->num() << endl;
  }
  return fieldStrings;
}

void ReadCFTimeDataFromNcFile(
	NcFile * ncfile,
	const std::string & strFilename,
	std::vector<Time> & vecTimes,
	bool fWarnOnMissingCalendar
) {
	// Empty existing Time vector
	vecTimes.clear();

	// Get time dimension
	NcDim * dimTime = ncfile->get_dim("time");
	if (dimTime == NULL) {
		_EXCEPTION1("Dimension \"time\" not found in file \"%s\"",
			strFilename.c_str());
	}

	// Get time variable
	NcVar * varTime = ncfile->get_var("time");
	if (varTime == NULL) {
		_EXCEPTION1("Variable \"time\" not found in file \"%s\"",
			strFilename.c_str());
	}
	if (varTime->num_dims() != 1) {
		_EXCEPTION1("Variable \"time\" has more than one dimension in file \"%s\"",
			strFilename.c_str());
	}
	if (strcmp(varTime->get_dim(0)->name(), "time") != 0) {
		_EXCEPTION1("Variable \"time\" does not have dimension \"time\" in file \"%s\"",
			strFilename.c_str());
	}

	// Calendar attribute
	NcAtt * attTimeCal = varTime->get_att("calendar");
	std::string strCalendar;
	if (attTimeCal == NULL) {
		if (fWarnOnMissingCalendar) {
			Announce("WARNING: Variable \"time\" is missing \"calendar\" attribute; assuming \"standard\"");
		}
		strCalendar = "standard";
	} else {
		strCalendar = attTimeCal->as_string(0);
	}
	Time::CalendarType eCalendarType =
		Time::CalendarTypeFromString(strCalendar);

	// Units attribute
	NcAtt * attTimeUnits = varTime->get_att("units");
	if (attTimeUnits == NULL) {
		_EXCEPTION1("Variable \"time\" is missing \"units\" attribute in file \"%s\"",
			strFilename.c_str());
	}
	std::string strTimeUnits = attTimeUnits->as_string(0);

	// Load in time data
	DataVector<int> vecTimeInt;
	DataVector<float> vecTimeFloat;
	DataVector<double> vecTimeDouble;
	DataVector<ncint64> vecTimeInt64;

	if (varTime->type() == ncInt) {
		vecTimeInt.Initialize(dimTime->size());
		varTime->set_cur((long)0);
		varTime->get(&(vecTimeInt[0]), dimTime->size());

	} else if (varTime->type() == ncFloat) {
		vecTimeFloat.Initialize(dimTime->size());
		varTime->set_cur((long)0);
		varTime->get(&(vecTimeFloat[0]), dimTime->size());

	} else if (varTime->type() == ncDouble) {
		vecTimeDouble.Initialize(dimTime->size());
		varTime->set_cur((long)0);
		varTime->get(&(vecTimeDouble[0]), dimTime->size());

	} else if (varTime->type() == ncInt64) {
		vecTimeInt64.Initialize(dimTime->size());
		varTime->set_cur((long)0);
		varTime->get(&(vecTimeInt64[0]), dimTime->size());

	} else {
		_EXCEPTION1("Variable \"time\" has invalid type "
			"(expected \"int\", \"int64\", \"float\" or \"double\")"
			" in file \"%s\"", strFilename.c_str());
	}

	for (int t = 0; t < dimTime->size(); t++) {
		Time time(eCalendarType);
		if (varTime->type() == ncInt) {
			time.FromCFCompliantUnitsOffsetInt(
				strTimeUnits,
				vecTimeInt[t]);

		} else if (varTime->type() == ncFloat) {
			time.FromCFCompliantUnitsOffsetDouble(
				strTimeUnits,
				static_cast<double>(vecTimeFloat[t]));

		} else if (varTime->type() == ncDouble) {
			time.FromCFCompliantUnitsOffsetDouble(
				strTimeUnits,
//.........这里部分代码省略.........

int main(int argc, char** argv) {

	NcError error(NcError::silent_nonfatal);

try {
	// Input filename
	std::string strInputFile;

	// Output mesh filename
	std::string strOutputFile;

	// Polynomial degree per element
	int nP = 2;

	// Parse the command line
	BeginCommandLine()
		CommandLineString(strInputFile, "in", "");
		CommandLineString(strOutputFile, "out", "");
		//CommandLineInt(nP, "np", 2);
		//CommandLineBool(fCGLL, "cgll");

		ParseCommandLine(argc, argv);
	EndCommandLine(argv)

	// Check file names
	if (strInputFile == "") {
		std::cout << "ERROR: No input file specified" << std::endl;
		return (-1);
	}
	if (strOutputFile == "") {
		std::cout << "ERROR: No output file specified" << std::endl;
		return (-1);
	}
	if (nP < 1) {
		std::cout << "ERROR: --np must be >= 2" << std::endl;
		return (-1);
	}

	AnnounceBanner();

	// Load input mesh
	AnnounceStartBlock("Loading input mesh");

	Mesh meshIn(strInputFile);
	meshIn.RemoveZeroEdges();

	AnnounceEndBlock("Done");

	// Construct edge map
	AnnounceStartBlock("Constructing edge map");

	meshIn.ConstructEdgeMap();

	AnnounceEndBlock("Done");

	// Build connectivity vector using edge map
	AnnounceStartBlock("Constructing connectivity");

    std::vector< std::set<int> > vecConnectivity;
    int err = GenerateConnectivityData(meshIn, vecConnectivity);
    if (err) return err;

	AnnounceEndBlock("Done");

	// Open output file
	AnnounceStartBlock("Writing connectivity file");

	NcFile ncmesh(strInputFile.c_str(), NcFile::ReadOnly);

	NcVar * varLat = ncmesh.get_var("grid_center_lat");
	NcVar * varLon = ncmesh.get_var("grid_center_lon");

	// Check if center latitudes and longitudes are already available
	DataArray1D<double> dAllLats;
	DataArray1D<double> dAllLons;

	bool fConvertLatToDegrees = true;
	bool fConvertLonToDegrees = true;

	if ((varLat == NULL) || (varLon == NULL)) {
		Announce("grid_center_lat not found, recalculating face centers");
	} else {
		Announce("grid_center_lat found in file, loading values");

		if (varLat->get_dim(0)->size() != vecConnectivity.size()) {
			_EXCEPTIONT("grid_center_lat dimension mismatch");
		}
		if (varLon->get_dim(0)->size() != vecConnectivity.size()) {
			_EXCEPTIONT("grid_center_lon dimension mismatch");
		}

		dAllLats.Allocate(vecConnectivity.size());
		varLat->set_cur((long)0);
		varLat->get(dAllLats, vecConnectivity.size());

		NcAtt * attLatUnits = varLat->get_att("units");
		std::string strLatUnits = attLatUnits->as_string(0);
		if (strLatUnits == "degrees") {
			fConvertLatToDegrees = false;
		}
//.........这里部分代码省略.........

ERMsg C20thReanalysisProject::ExtractTopo(const CString& filePath, CSCCallBack& callback)
{
	ERMsg msg;

	

	CString outputFilePath(filePath);
	UtilWin::SetFileExtension( outputFilePath, ".csv");
	
	CString outputFilePathGrid(filePath);
	UtilWin::SetFileExtension( outputFilePathGrid, ".tif");
	
	CString filePathIn;
	if(m_type==_2X2)
		filePathIn.Format("%shgt.2x2.nc", m_path);
	else filePathIn.Format("%shgt.gauss.nc", m_path);

	int sizeX=0;
	int sizeY=0;
	int nbBand=0;
	CGeoRectWP rect;
	float noData=0;
	double cellSizeX=0;
	double cellSizeY=0;

	msg = GetGridInfo(filePathIn, sizeX, sizeY, nbBand, rect, noData, cellSizeX, cellSizeY);
	CMFCGDALDataset dataset;

	rect.m_xMax+=cellSizeX;
	dataset.Create(outputFilePathGrid, sizeX+1, sizeY, rect, "GTiff", GDT_Int16, 1, -9999);
	

	NcFile file(filePathIn);	
	
	int nbDim = file.num_dims();
	int nbAtt = file.num_atts();
	int nbVar = file.num_vars();

	//NcDim* dim0 = file.get_dim(0);
	//NcDim* dim1 = file.get_dim(1);
	//NcDim* dim2 = file.get_dim(2);

	if( !file.is_valid() )
	{
		CString err;
		err.FormatMessage(IDS_CMN_UNABLE_OPEN_READ, filePath);
		msg.ajoute(err);
	}

	callback.SetCurrentDescription("Create Topo");
	callback.SetCurrentStepRange(0, 17*28, 1);
	callback.SetStartNewStep();

	CStdioFileEx fileOut;


	msg += fileOut.Open( outputFilePath, CFile::modeCreate|CFile::modeWrite);
	if(!msg)
		return msg;

	fileOut.WriteString("Name,ID,Latitude,Longitude,Elevation\n");
	
//	NcVar* pVarX = file.get_var("xc");
	//NcVar* pVarY = file.get_var("yc");
//	NcDim* pDim = file.get_dim("lat");
	//int size = pDim->size();
	NcDim* pDimX = file.get_dim("lon");
	NcDim* pDimY = file.get_dim("lat");
	NcDim* pDimTime = file.get_dim("time");
	

		
	NcVar* pVarLat = file.get_var("lat");
	NcVar* pVarLon = file.get_var("lon");
	NcVar* pVarElev = file.get_var("hgt");
	float offset = pVarElev->get_att("add_offset")->as_float(0);
	float scaleFactor = pVarElev->get_att("scale_factor")->as_float(0);
	//int sizeX2 = pDimX->size();
	//int sizeY2 = pDimY->size();
	
	//typedef boost::multi_array<short, 3> CData3Array;
	//typedef boost::multi_array<float, 2> CDataFloatArray;
	vector<float> X(sizeX);
	vector<float> Y(sizeY);
//	CDataFloatArray X(boost::extents[180]);
//	CDataFloatArray Y(boost::extents[172]);
	//CData2Array lat(boost::extents[28][17]);
	//CData2Array lon(boost::extents[28][17]);
	boost::multi_array<short, 2> elev(boost::extents[sizeY][sizeX]);
	
	
//	pVarX->get(&(X[0]), 180);
//	pVarY->get(&(Y[0]), 172);
	pVarLon->get(&(X[0]), sizeX);
	pVarLat->get(&(Y[0]), sizeY);
	pVarElev->get(&(elev[0][0]), 1, sizeY, sizeX);

//	CProjection prj = GetDataGrid().GetPrj();

	CGridLine profile(sizeX+1);
//.........这里部分代码省略.........

vtkDataArray *
avtS3DFileFormat::GetVar(int timeState, int domain, const char *varname)
{
    debug5 << "avtS3DFileFormat::GetVar( timeState=" << timeState << ", domain="
        << domain << ", varname=" << varname << ")" << endl;

    // Calculate the timestep directory that the data lives in.
    char *pathcopy = strdup(mainFilename);
    string dir = parse_dirname(pathcopy);
    string timestepDir = CreateStringFromDouble(fileTimes[timeState]);
    debug4 << "Timestep directory is <" << timestepDir <<  ">" << endl;
    
    // Figure out how big this piece is.
    CalculateSubpiece(domain);

    // Open up the NetCDF file.
    char path[256];
    SNPRINTF(path,256,"%s%s%s%sfield.%05d",dir.c_str(),VISIT_SLASH_STRING, timestepDir.c_str(), VISIT_SLASH_STRING, domain);
    debug5 << "avtS3DFileFormat::GetVar: Full path to data file is " << path << endl;

    NcFile nf(path);
    if (!nf.is_valid())
    {
        debug1 << nc_strerror(NcError().get_err()) << endl;
        EXCEPTION1(InvalidFilesException, path);
    }
    debug5 << "avtS3DFileFormat::GetVar: Got valid file." << endl;

    // Pull out the appropriate variable.
    NcVar *v = nf.get_var(varname);
    if (!v)
    {
        debug1 << nc_strerror(NcError().get_err()) << endl;
        EXCEPTION1(InvalidVariableException, varname);
    }

    // Check if it fits the size of the mesh.  Always node-centered, remember.
    int ntuples = localDims[0] * localDims[1] * localDims[2];
    debug5 << "ntuples:" << ntuples << endl;
    int nvals = v->num_vals();
    if (ntuples != nvals)
    {
        debug1 << "The variable " << v->name() <<
                  " does not conform to its mesh (" << nvals << " != " << 
                  ntuples << ")" << endl;
        EXCEPTION1(InvalidVariableException, v->name());
    }

    // Set up the VTK dataset.
    vtkFloatArray *rv = vtkFloatArray::New();
    rv->SetNumberOfTuples(ntuples);
    float *p = (float*)rv->GetVoidPointer(0);
    NcValues *input = v->values();
    if (!input)
    {
        debug1 << nc_strerror(NcError().get_err()) << endl;
        EXCEPTION1(InvalidVariableException, v->name());
    }

    // Get the scaling factor.
    NcAtt *scaling = v->get_att(NcToken("scale_factor"));
    float scaling_factor = 1;
    if (scaling)
    {
        scaling_factor = scaling->as_float(0);
        debug5 << "avtS3DFileFormat::GetVar: Set the scaling factor as " << scaling_factor << endl;
    }

    // Process the variable into the returned data.
    float *base = (float*)input->base();
    for(int i=0;i<ntuples;i++)
    {
        p[i] = *(base + i) * scaling_factor;
    }

    return rv;
}

void
avtS3DFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md,
                                           int timeState)
{
    debug5 << "avtS3DFileFormat::PopulateDatabaseMetaData" << endl;
    // Get the metadata from the log file first.
    OpenLogFile();

    // Mesh
    avtMeshMetaData *mesh = new avtMeshMetaData;
    mesh->name = "mesh";
    mesh->meshType = AVT_RECTILINEAR_MESH;
    mesh->numBlocks = procs[0] * procs[1] * procs[2];
    mesh->blockOrigin = 1;
    mesh->cellOrigin = 0;
    mesh->spatialDimension = 3;
    mesh->topologicalDimension = 3;
    mesh->blockTitle = "blocks";
    mesh->blockPieceName = "block";
    mesh->hasSpatialExtents = false;
    mesh->xUnits = "mm";
    mesh->yUnits = "mm";
    mesh->zUnits = "mm";
    md->Add(mesh);

    //
    // Look in the NetCDF file for the first block for the list of variables.
    //

    // Calculate the timestep directory that the data lives in.
    char *pathcopy = strdup(mainFilename);
    string dir = parse_dirname(pathcopy);
    string timestepDir = CreateStringFromDouble(fileTimes[timeState]);

    char path[256];
    SNPRINTF(path,256,"%s%s%s%sfield.00000",dir.c_str(),VISIT_SLASH_STRING, timestepDir.c_str(), VISIT_SLASH_STRING);

    NcError err(NcError::verbose_nonfatal);
 
    NcFile nf(path);
    if (!nf.is_valid())
    {
        EXCEPTION1(InvalidFilesException, path);
    }
    debug5 << "avtS3DFileFormat::PopulateDatabaseMetaData: Got valid file" << endl;

    int nvars = nf.num_vars();
    debug5 << "avtS3DFileFormat::PopulateDatabaseMetaData: Found " << nvars << " variables" << endl;
    for (int i=0 ; i<nvars; i++)
    {
        NcVar *v = nf.get_var(i);
        if (!v)
            continue;
        debug4 << "Found variable " << v->name() << endl;

        // Check dimensionality
        int nvals = v->num_vals();
        if (nvals != 1) // Single scalars are useless.
        {
            avtScalarMetaData *scalar = new avtScalarMetaData();
            scalar->name = v->name();
            scalar->meshName = "mesh";
            scalar->centering = AVT_NODECENT;
            scalar->hasDataExtents = false;
            scalar->treatAsASCII = false;

            NcAtt *units = v->get_att(NcToken("units"));
            if (units)
            {
                long nv = units->num_vals();
                if (nv == 0)
                {
                    scalar->hasUnits = false;
                } else {
                    char *unitString = units->as_string(0);
                    scalar->units = unitString;
                    scalar->hasUnits = true;
                }
            } else
                scalar->hasUnits = false;

            md->Add(scalar);
        } else {
            debug4 << "Unable to process variable " << v->name() <<
                      " since it is a single scalar" << endl;

        }
    }

#if 0
    // Expressions
    Expression tempGradient_expr;
    tempGradient_expr.SetName("Temperature_gradient");
    tempGradient_expr.SetDefinition("gradient(Temperature)");
    tempGradient_expr.SetType(Expression::VectorMeshVar);
    tempGradient_expr.SetHidden(true);
    md->AddExpression(&tempGradient_expr);

    Expression tempUnit_expr;
    tempUnit_expr.SetName("Temperature_grad_unit");
//.........这里部分代码省略.........