C++ NcAtt类代码示例

NcAtt* NcVar::get_att( NcToken aname ) const
{
    NcAtt* att = new NcAtt(the_file, this, aname);
    if (! att->is_valid()) {
	delete att;
	return 0;
    }
    return att;
}

void DumpableNcFile::dumpgatts( void )
{
    NcAtt* ap;
    for(int n = 0; ap = get_att(n); n++) {
	cout << "\t\t" << ":" << ap->name() << " = " ;
	NcValues* vals = ap->values();
	cout << *vals << " ;" << endl ;
	delete vals;
	delete ap;
    }
}

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;
    }
}

string GrabAttribute(const NcVar* dataVar, const int attIndex)
{
    if (NULL == dataVar)
    {
        cerr << "ERROR: Invalid NcVar!\n";
        return("");
    }

    NcAtt* dataAttrib = dataVar->get_att(attIndex);

    if (NULL == dataAttrib)
    {
        cerr << "ERROR: Could not obtain data attribute: index #" << attIndex << endl;
        return("");
    }

    const string returnVal = dataAttrib->as_string(0);

    delete dataAttrib;

    return(returnVal);
}

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");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    {
        cerr << "ERROR: invalid data file.  No variable called 'lat'!\n";
        radarFile.close();
        return(inputData);
    }

    long latCnt = latVar->num_vals();
    double* latVals = new double[latCnt];
    latVar->get(latVals, latCnt);
    inputData.latUnits = GrabAttribute(latVar, 0);
    inputData.latSpacing = strtod(GrabAttribute(latVar, 1).c_str(), NULL);

    const long minLatIndex = lower_bound(latVals, latmin, latCnt);
    const long maxLatIndex = upper_bound(latVals, latmax, latCnt);

    delete latVals;
    latCnt = (maxLatIndex - minLatIndex) + 1;
    latVar->set_cur(minLatIndex);
    inputData.latVals = new double[latCnt];
    latVar->get(inputData.latVals, latCnt);


    NcVar* lonVar = radarFile.get_var("lon");

    if (NULL == lonVar)
    {
        cerr << "ERROR: invalid data file.  No variable called 'lon'!\n";
        radarFile.close();
        return(inputData);
    }

    long lonCnt = lonVar->num_vals();
    double* lonVals = new double[lonCnt];
    lonVar->get(lonVals, lonCnt);
    inputData.lonUnits = GrabAttribute(lonVar, 0);
    inputData.lonSpacing = strtod(GrabAttribute(lonVar, 1).c_str(), NULL);

    const long minLonIndex = lower_bound(lonVals, lonmin, lonCnt);
    const long maxLonIndex = upper_bound(lonVals, lonmax, lonCnt);

    delete lonVals;
    lonCnt = (maxLonIndex - minLonIndex) + 1;
    lonVar->set_cur(minLonIndex);
    inputData.lonVals = new double[lonCnt];
    lonVar->get(inputData.lonVals, lonCnt);

    NcVar* reflectVar = NULL;
    reflectVar = radarFile.get_var("value");

    if ( reflectVar == NULL )
    {
        // Try this variable name
        reflectVar = radarFile.get_var("Reflectivity");
    }

    if (reflectVar == NULL)
    {
        cerr << "ERROR: invalid data file.  No variable called 'value'!\n";
        radarFile.close();
        return(inputData);
    }


    inputData.dataEdges = reflectVar->edges();	// [0] - time, [1] - lat, [2] - lon
    inputData.dataEdges[1] = latCnt;
    inputData.dataEdges[2] = lonCnt;
    inputData.dataVals = new double[inputData.dataEdges[0] * inputData.dataEdges[1] * inputData.dataEdges[2]];
    reflectVar->set_cur(0, minLatIndex, minLonIndex);
    reflectVar->get(inputData.dataVals, inputData.dataEdges);

    inputData.var_LongName = GrabAttribute(reflectVar, 0);
    inputData.var_Units = "dBZ";//GrabAttribute(reflectVar, 1);


    NcVar* timeVar = radarFile.get_var("time");

    if (NULL == timeVar)
    {
        cerr << "ERROR: invalid data file.  No variable called 'time'!\n";
        radarFile.close();
        return(inputData);
    }

    inputData.scanTime = timeVar->as_long(0);
    inputData.timeUnits = GrabAttribute(timeVar, 0);

    NcAtt* titleAttrib = radarFile.get_att("title");

    inputData.fileTitle = (NULL == titleAttrib ? "" : titleAttrib->as_string(0));

    delete titleAttrib;

    radarFile.close();

    // Success!
    inputData.inputFilename = filename;


    return(inputData);
}

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,
//.........这里部分代码省略.........

void CopyNcVarAttributes(
	NcVar * varIn,
	NcVar * varOut
) {
	for (int a = 0; a < varIn->num_atts(); a++) {
		NcAtt * att = varIn->get_att(a);
		long num_vals = att->num_vals();

		NcValues * pValues = att->values();

		if (att->type() == ncByte) {
			varOut->add_att(att->name(), num_vals,
				(const ncbyte*)(pValues->base()));

		} else if (att->type() == ncChar) {
			varOut->add_att(att->name(), num_vals,
				(const char*)(pValues->base()));

		} else if (att->type() == ncShort) {
			varOut->add_att(att->name(), num_vals,
				(const short*)(pValues->base()));

		} else if (att->type() == ncInt) {
			varOut->add_att(att->name(), num_vals,
				(const int*)(pValues->base()));

		} else if (att->type() == ncFloat) {
			varOut->add_att(att->name(), num_vals,
				(const float*)(pValues->base()));

		} else if (att->type() == ncDouble) {
			varOut->add_att(att->name(), num_vals,
				(const double*)(pValues->base()));

		} else {
			_EXCEPTIONT("Invalid attribute type");
		}

		delete pValues;
	}
}

//
// Reads variable data from dump file
//
bool DataVar::initFromFile(const string& filename, const_DomainChunk_ptr dom)
{
    cleanup();
    
#if ESYS_HAVE_NETCDF
    NcError ncerr(NcError::silent_nonfatal);    
    NcFile* input = new NcFile(filename.c_str());
    if (!input->is_valid()) {
        cerr << "Could not open input file " << filename << "." << endl;
        delete input;
        return false;
    }

    NcDim* dim;
    NcAtt* att;

    att = input->get_att("type_id");
    int typeID = att->as_int(0);
    if (typeID != 2) {
        cerr << "WARNING: Only expanded data supported!" << endl;
        delete input;
        return false;
    }

    att = input->get_att("rank");
    rank = att->as_int(0);

    dim = input->get_dim("num_data_points_per_sample");
    ptsPerSample = dim->size();

    att = input->get_att("function_space_type");
    funcSpace = att->as_int(0);

    centering = dom->getCenteringForFunctionSpace(funcSpace);

    dim = input->get_dim("num_samples");
    numSamples = dim->size();

#ifdef _DEBUG
    cout << varName << ":\t" << numSamples << " samples,  "
        << ptsPerSample << " pts/s,  rank: " << rank << endl;
#endif

    domain = dom;
    NodeData_ptr nodes = domain->getMeshForFunctionSpace(funcSpace);
    if (nodes == NULL) {
        delete input;
        return false;
    }

    meshName = nodes->getName();
    siloMeshName = nodes->getFullSiloName();
    initialized = true;

    size_t dimSize = 1;
    vector<long> counts;

    if (rank > 0) {
        dim = input->get_dim("d0");
        int d = dim->size();
        shape.push_back(d);
        counts.push_back(d);
        dimSize *= d;
    }
    if (rank > 1) {
        dim = input->get_dim("d1");
        int d = dim->size();
        shape.push_back(d);
        counts.push_back(d);
        dimSize *= d;
    }
    if (rank > 2) {
        cerr << "WARNING: Rank " << rank << " data is not supported!\n";
        initialized = false;
    }
 
    if (initialized && numSamples > 0) {
        sampleID.insert(sampleID.end(), numSamples, 0);
        NcVar* var = input->get_var("id");
        var->get(&sampleID[0], numSamples);

        size_t dataSize = dimSize*numSamples*ptsPerSample;
        counts.push_back(ptsPerSample);
        counts.push_back(numSamples);
        float* tempData = new float[dataSize];
        var = input->get_var("data");
        var->get(tempData, &counts[0]);

        const float* srcPtr = tempData;
        for (size_t i=0; i < dimSize; i++, srcPtr++) {
            float* c = averageData(srcPtr, dimSize);
            dataArray.push_back(c);
        }
        delete[] tempData;

        initialized = reorderSamples();
    }
//.........这里部分代码省略.........

int main(void)
{
   // These will hold our pressure and temperature data.
   float presIn[NLAT][NLON];
   float tempIn[NLAT][NLON];

   // These will hold our latitudes and longitudes.
   float latsIn[NLAT];
   float lonsIn[NLON];
  
   // Change the error behavior of the netCDF C++ API by creating an
   // NcError object. Until it is destroyed, this NcError object will
   // ensure that the netCDF C++ API silently returns error codes on
   // any failure, and leaves any other error handling to the calling
   // program. In the case of this example, we just exit with an
   // NC_ERR error code.
   NcError err(NcError::silent_nonfatal);

   // Open the file and check to make sure it's valid.
   NcFile dataFile("sfc_pres_temp.nc", NcFile::ReadOnly);
   if(!dataFile.is_valid())
      return NC_ERR;
  
   // There are a number of inquiry functions in netCDF which can be
   // used to learn about an unknown netCDF file. In this case we know
   // that there are 2 netCDF dimensions, 4 netCDF variables, no
   // global attributes, and no unlimited dimension.
   if (dataFile.num_dims() != 2 || dataFile.num_vars() != 4 || 
       dataFile.num_atts() != 0 || dataFile.rec_dim() != 0)  
      return NC_ERR;
     
   // We get back a pointer to each NcVar we request. Get the
   // latitude and longitude coordinate variables.
   NcVar *latVar, *lonVar;
   if (!(latVar = dataFile.get_var("latitude")))
      return NC_ERR;
   if (!(lonVar = dataFile.get_var("longitude")))
      return NC_ERR;
       
   // Read the latitude and longitude coordinate variables into arrays
   // latsIn and lonsIn.
   if (!latVar->get(latsIn, NLAT))
      return NC_ERR;
   if (!lonVar->get(lonsIn, NLON))
      return NC_ERR;
       
   // Check the coordinate variable data. 
   for(int lat = 0; lat < NLAT; lat++)
      if (latsIn[lat] != START_LAT + 5. * lat)
	 return NC_ERR;
   
   // Check longitude values.
   for (int lon = 0; lon < NLON; lon++)
      if (lonsIn[lon] != START_LON + 5. * lon)
	 return NC_ERR;
   
   // We get back a pointer to each NcVar we request. 
   NcVar *presVar, *tempVar;
   if (!(presVar = dataFile.get_var("pressure")))
      return NC_ERR;
   if (!(tempVar = dataFile.get_var("temperature")))
      return NC_ERR;
       
   // Read the data. Since we know the contents of the file we know
   // that the data arrays in this program are the correct size to
   // hold all the data.
   if (!presVar->get(&presIn[0][0], NLAT, NLON))
      return NC_ERR;
   if (!tempVar->get(&tempIn[0][0], NLAT, NLON))
      return NC_ERR;
       
   // Check the data. 
   for (int lat = 0; lat < NLAT; lat++)
      for (int lon = 0; lon < NLON; lon++)
	 if (presIn[lat][lon] != SAMPLE_PRESSURE + (lon * NLAT + lat)
	     || tempIn[lat][lon] != SAMPLE_TEMP + .25 * (lon * NLAT + lat))
	    return NC_ERR;
   
   // Each of the netCDF variables has a "units" attribute. Let's read
   // them and check them.
   NcAtt *att;
   char *units;
   
   if (!(att = latVar->get_att("units")))
      return NC_ERR;
   units = att->as_string(0);
   if (strncmp(units, "degrees_north", strlen("degrees_north")))
      return NC_ERR;
   // Attributes and attribute values should be deleted by the caller
   // when no longer needed, to prevent memory leaks.
   delete units;
   delete att;

   if (!(att = lonVar->get_att("units")))
      return NC_ERR;
   units = att->as_string(0);
   if (strncmp(units, "degrees_east", strlen("degrees_east")))
      return NC_ERR;
   delete units;
   delete att;
//.........这里部分代码省略.........

bool
RemapWidget::getVolumeInfo(QString volfile,
			   int &skipheaderbytes,
			   uchar &voxelType,
			   int &voxelUnit,
			   float &vx, float &vy, float &vz,
			   QString &description,
			   QList<float> &rawMap,
			   QList<uchar> &pvlMap,
			   int &depth,
			   int &width,
			   int &height)
{
  NcError err(NcError::verbose_nonfatal);

  NcFile pvlFile(volfile.toAscii().data(), NcFile::ReadOnly);

  if (!pvlFile.is_valid())
    {
      QMessageBox::information(0, "Error",
       QString("%1 is not a valid preprocessed volume file").arg(volfile));
      return false;
    }

  int i;
  NcAtt *att;
  char *attval;
  QString pvalue;

  att = pvlFile.get_att("description");
  if (att)
    {
      attval = att->as_string(0);
      description = attval;
      delete [] attval;
    }

  att = pvlFile.get_att("voxeltype");
  if (att)
    {
      attval = att->as_string(0);
      pvalue = attval;
      if (pvalue == "unsigned char")
	voxelType = Raw2Pvl::_UChar;
      if (pvalue == "char")
	voxelType = Raw2Pvl::_Char;
      if (pvalue == "unsigned short")
	voxelType = Raw2Pvl::_UShort;
      if (pvalue == "short")
	voxelType = Raw2Pvl::_Short;
      if (pvalue == "int")
	voxelType = Raw2Pvl::_Int;
      if (pvalue == "float")
	voxelType = Raw2Pvl::_Float;
      delete [] attval;
    }


  att = pvlFile.get_att("voxelunit");
  if (att)
    { 
      attval = att->as_string(0);
      pvalue = attval;
      voxelUnit = Raw2Pvl::_Nounit;
      if (pvalue == "angstrom")
	voxelUnit = Raw2Pvl::_Angstrom;
      else if (pvalue == "nanometer")
	voxelUnit = Raw2Pvl::_Nanometer;
      else if (pvalue == "micron")
	voxelUnit = Raw2Pvl::_Micron;
      else if (pvalue == "millimeter")
	voxelUnit = Raw2Pvl::_Millimeter;
      else if (pvalue == "centimeter")
	voxelUnit = Raw2Pvl::_Centimeter;
      else if (pvalue == "meter")
	voxelUnit = Raw2Pvl::_Meter;
      else if (pvalue == "kilometer")
	voxelUnit = Raw2Pvl::_Kilometer;
      else if (pvalue == "parsec")
	voxelUnit = Raw2Pvl::_Parsec;
      else if (pvalue == "kiloparsec")
	voxelUnit = Raw2Pvl::_Kiloparsec;
      delete [] attval;
    }
  
  
  att = pvlFile.get_att("gridsize");
  if (att)
    {
      depth = att->as_int(0);
      width = att->as_int(1);
      height = att->as_int(2);
    }

  att = pvlFile.get_att("voxelsize");
  if (att)
    {
      vx = att->as_float(0);
      vy = att->as_float(1);
      vz = att->as_float(2);
//.........这里部分代码省略.........

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;
}

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;
		}
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	varOutLon->set_cur((long)0);
	varOutLon->put(&(dLon[0]), nLon);

	CopyNcVarAttributes(varLon, varOutLon);

	// Output topography
	Announce("Topography");
	NcVar * varOutZs = ncdf_out.add_var(
		"Zs", ncDouble, dimOutLat, dimOutLon);

	varOutZs->set_cur((long)0, (long)0);
	varOutZs->put(&(dZs[0][0]), nLat, nLon);

	AnnounceEndBlock("Done");

	// Done
	AnnounceEndBlock("Done");

	// Load all variables
	Announce("Loading variables");

	std::vector<NcVar *> vecNcVar;
	for (int v = 0; v < vecVariableStrings.size(); v++) {
		vecNcVar.push_back(ncdf_in.get_var(vecVariableStrings[v].c_str()));
		if (vecNcVar[v] == NULL) {
			_EXCEPTION1("Unable to load variable \"%s\" from file",
				vecVariableStrings[v].c_str());
		}
	}

	// Physical constants
	Announce("Initializing thermodynamic variables");

	NcAtt * attEarthRadius = ncdf_in.get_att("earth_radius");
	double dEarthRadius = attEarthRadius->as_double(0);

	NcAtt * attRd = ncdf_in.get_att("Rd");
	double dRd = attRd->as_double(0);

	NcAtt * attCp = ncdf_in.get_att("Cp");
	double dCp = attCp->as_double(0);

	double dGamma = dCp / (dCp - dRd);

	NcAtt * attP0 = ncdf_in.get_att("P0");
	double dP0 = attP0->as_double(0);

	double dPressureScaling = dP0 * std::pow(dRd / dP0, dGamma);

	NcAtt * attZtop = ncdf_in.get_att("Ztop");
	double dZtop = attZtop->as_double(0);

	// Input data
	DataArray3D<double> dataIn(nLev, nLat, nLon);
	DataArray3D<double> dataInt(nILev, nLat, nLon);

	// Output data
	DataArray2D<double> dataOut(nLat, nLon);

	// Pressure in column
	DataArray1D<double> dataColumnP(nLev);

	// Height in column
	DataArray1D<double> dataColumnZ(nLev);
	DataArray1D<double> dataColumnIZ(nILev);

void PolygonManager::init()
{
    string fileName;
    ConfigTools::read("parcel_polygon_file", fileName);
    NOTICE("PolygonManager::init", "Reading polygons from \""+fileName+"\" ...");
    NcFile file(fileName.c_str(), NcFile::ReadOnly);
    if (!file.is_valid()) {
        REPORT_ERROR(string("Failed to open file "+fileName+"."))
    }

    NcError ncError(NcError::silent_nonfatal);

    if (TimeManager::onLine()) {
        NcAtt *timeAtt = file.get_att("time");
        NcAtt *timeStepAtt = file.get_att("time_step");
        NcAtt *stepsAtt = file.get_att("steps");
        if (timeAtt != NULL && timeStepAtt != NULL && stepsAtt != NULL) {
            TimeManager::reset();
            double dt = timeStepAtt->as_double(0);
            double second = timeAtt->as_double(0);
            double steps = stepsAtt->as_int(0);
            TimeManager::setClock(dt, second, steps);
        }
    }
    NcDim *numVertexDim = file.get_dim("num_total_vertex");
    if (numVertexDim == NULL) {
        Message message;
        message << "Failed to find \"num_total_vertex\" dimension in file \"";
        message << fileName << "\"!";
        REPORT_ERROR(message.str());
    }
    NcDim *numEdgeDim = file.get_dim("num_total_edge");
    if (numEdgeDim == NULL) {
        Message message;
        message << "Failed to find \"num_total_edge\" dimension in file \"";
        message << fileName << "\"!";
        REPORT_ERROR(message.str());
    }
    NcDim *numPolygonDim = file.get_dim("num_total_polygon");
    if (numPolygonDim == NULL) {
        Message message;
        message << "Failed to find \"num_total_polygon\" dimension in file \"";
        message << fileName << "\"!";
        REPORT_ERROR(message.str());
    }
    int numVertex = static_cast<int>(numVertexDim->size());
    int numEdge = static_cast<int>(numEdgeDim->size());
    int numPolygon = static_cast<int>(numPolygonDim->size());
    // -------------------------------------------------------------------------
    // vertices part
    vertices.create(numVertex);
    double *oldVtxLon = new double[numVertex];
    double *oldVtxLat = new double[numVertex];
    double *oldVtxLev = new double[numVertex];
    double *newVtxLon = new double[numVertex];
    double *newVtxLat = new double[numVertex];
    double *newVtxLev = new double[numVertex];
    file.get_var("old_vertex_lon")->get(oldVtxLon, numVertex);
    file.get_var("old_vertex_lat")->get(oldVtxLat, numVertex);
    file.get_var("new_vertex_lon")->get(newVtxLon, numVertex);
    file.get_var("new_vertex_lat")->get(newVtxLat, numVertex);
    if (file.get_var("old_vertex_lev") != NULL) {
        file.get_var("old_vertex_lev")->get(oldVtxLev, numVertex);
        file.get_var("new_vertex_lev")->get(newVtxLev, numVertex);
    } else {
        for (int i = 0; i < vertices.size(); ++i) {
            oldVtxLev[i] = 0.0;
            newVtxLev[i] = 0.0;
        }
    }
    // change the units from degree to rad
    for (int i = 0; i < numVertex; ++i) {
        oldVtxLon[i] /= Rad2Deg;
        oldVtxLat[i] /= Rad2Deg;
        newVtxLon[i] /= Rad2Deg;
        newVtxLat[i] /= Rad2Deg;
    }
    Vertex *vertexMap[vertices.size()];
    Vertex *vertex = vertices.front();
    for (int i = 0; i < vertices.size(); ++i) {
        vertexMap[i] = vertex;
        vertex->setCoordinate(newVtxLon[i], newVtxLat[i], newVtxLev[i], NewTimeLevel);
        vertex->setCoordinate(oldVtxLon[i], oldVtxLat[i], oldVtxLev[i], OldTimeLevel);
        vertex = vertex->next;
    }
    delete [] newVtxLon;
    delete [] newVtxLat;
    delete [] newVtxLev;
    delete [] oldVtxLon;
    delete [] oldVtxLat;
    delete [] oldVtxLev;
    // -------------------------------------------------------------------------
    // edges part
    edges.create(numEdge);
    int *firstPoint = new int[numEdge];
    int *secondPoint = new int[numEdge];
    file.get_var("first_point_idx")->get(firstPoint, numEdge);
    file.get_var("second_point_idx")->get(secondPoint, numEdge);
    Edge *edgeMap[edges.size()];
    Edge *edge = edges.front();
//.........这里部分代码省略.........