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

void ConstantSet::netcdf_define(NcFile &nc, std::string const &vname)
{
	// Create the variable to store our constants
	NcVar *ncvar = giss::get_var_safe(nc, vname.c_str(), false);
	if (!ncvar) {
		auto oneDim = giss::get_or_add_dim(nc, "one", 1);
		ncvar = nc.add_var(vname.c_str(), ncDouble, oneDim);
	}

	// Store the constants as attributes
	for (int i=0; i<fields.size_withunit(); ++i) {
		CoupledField const &field(fields.field(i));
		ncvar->add_att(field.name.c_str(), vals[i]);
		ncvar->add_att((field.name + "_units").c_str(), field.units.c_str());
		ncvar->add_att((field.name + "_description").c_str(), field.description.c_str());
	}
}

void PointCounter::output(const string &fileName) const
{
    NcFile file(fileName.c_str(), NcFile::Replace);
    if (!file.is_valid()) {
        string message = "Failed to open file "+fileName+".";
        REPORT_ERROR(message.c_str());
    }

    NcDim *numLonDim = file.add_dim("lon", counters.extent(0));
    NcDim *numLatDim = file.add_dim("lat", counters.extent(1));
    NcDim *numBndsDim = file.add_dim("bnds", 2);
    //NcDim *numLevDim = file.add_dim("lev", counters.extent(2));

    NcVar *lonBndsVar = file.add_var("lon_bnds", ncDouble, numLonDim, numBndsDim);
    NcVar *latBndsVar = file.add_var("lat_bnds", ncDouble, numLatDim, numBndsDim);
    //NcVar *levVar = file.add_var("lev", ncDouble, numLevDim);

    lonBndsVar->add_att("long_name", "longitude bounds");
    lonBndsVar->add_att("units", "degree_east");
    latBndsVar->add_att("long_name", "latitude bounds");
    latBndsVar->add_att("units", "degree_north");

    int numLonBnds = mesh[0].getNumLon()-2, numLatBnds = mesh[0].getNumLat()-1;
    double lonBnds[numLonBnds][2], latBnds[numLatBnds][2];
    for (int i = 0; i < numLonBnds; ++i) {
        lonBnds[i][0] = mesh[0].lon(i)*Rad2Deg;
        lonBnds[i][1] = mesh[0].lon(i+1)*Rad2Deg;
    }
    for (int j = 0; j < numLatBnds; ++j) {
        latBnds[j][0] = mesh[0].lat(j)*Rad2Deg;
        latBnds[j][1] = mesh[0].lat(j+1)*Rad2Deg;
    }

    lonBndsVar->put(&lonBnds[0][0], numLonBnds, 2);
    latBndsVar->put(&latBnds[0][0], numLatBnds, 2);

    NcVar *countersVar = file.add_var("counters", ncInt, numLatDim, numLonDim);

    countersVar->add_att("long_name", "contained point numbers in grid boxes");

    int counters[this->counters.extent(1)][this->counters.extent(0)];
    for (int j = 0; j < this->counters.extent(1); ++j)
        for (int i = 0; i < this->counters.extent(0); ++i)
            counters[j][i] = this->counters(i, j, 0);
    countersVar->put(&counters[0][0], this->counters.extent(1), this->counters.extent(0));
}

bool OutputClusters(const string &filename, const vector<Cluster> &theClusters,
                    const RadarData_t &radarInfo, const ClustParams_t &clustParam)
{

    NcFile clusterFile(filename.c_str(), NcFile::Replace);

    if (!clusterFile.is_valid())
    {
        cerr << "Could not create file: " << filename << "\n";
        return(false);
    }

    size_t pixelCnt = 0;
    for (vector<Cluster>::const_iterator aClust = theClusters.begin();
            aClust != theClusters.end();
            aClust++)
    {
        pixelCnt += aClust->size();
    }

    NcDim* pixelDim = clusterFile.add_dim("pixel_index", pixelCnt);
    NcDim* latDim = clusterFile.add_dim("lat", radarInfo.dataEdges[1]);
    NcDim* lonDim = clusterFile.add_dim("lon", radarInfo.dataEdges[2]);

    NcVar* clustMember = clusterFile.add_var("clusterIndex", ncInt, pixelDim);
    clustMember->add_att("Description", "pixel cluster membership");
    NcVar* xLoc = clusterFile.add_var("pixel_xLoc", ncInt, pixelDim);
    xLoc->add_att("Description", "x-index of pixel in domain (lat,lon)");
    NcVar* yLoc = clusterFile.add_var("pixel_yLoc", ncInt, pixelDim);
    yLoc->add_att("Description", "y-index of pixel in domain (lat,lon)");


    NcVar* lats = clusterFile.add_var("lat", ncDouble, latDim);
    lats->add_att("units", radarInfo.latUnits.c_str());
    lats->add_att("spacing", radarInfo.latSpacing);

    NcVar* lons = clusterFile.add_var("lon", ncDouble, lonDim);
    lons->add_att("units", radarInfo.lonUnits.c_str());
    lons->add_att("spacing", radarInfo.lonSpacing);

    clusterFile.add_att("title", "SPA Cluster Results of Radar Reflectivities");
    clusterFile.add_att("data_source", radarInfo.inputFilename.c_str());
    clusterFile.add_att("time", radarInfo.scanTime);

    // Need to do type-casting to double because the netcdf libraries can't seem to properly save a float.
    // Also, the truncf() is used to -- somewhat -- handle mantisa issues.
    clusterFile.add_att("Upper_Sensitivity", (double) (truncf(100.0 * clustParam.upperSensitivity)/100.0));
    clusterFile.add_att("Lower_Sensitivity", (double) (truncf(100.0 * clustParam.lowerSensitivity)/100.0));
    clusterFile.add_att("Padding_Level", (double) (truncf(100.0 * clustParam.paddingLevel)/100.0));
    clusterFile.add_att("Reach", (double) (truncf(100.0 * clustParam.reach)/100.0));
    clusterFile.add_att("Subcluster_Depth", clustParam.subClustDepth);

    int* clustIndicies = new int[pixelCnt];
    int* xPixels = new int[pixelCnt];
    int* yPixels = new int[pixelCnt];

    size_t pixelIndex = 0;
    for (size_t clustIndex = 0; clustIndex < theClusters.size(); clustIndex++)
    {
        for (Cluster::const_iterator aMember = theClusters[clustIndex].begin();
                aMember != theClusters[clustIndex].end();
                aMember++, pixelIndex++)
        {
            clustIndicies[pixelIndex] = clustIndex;
            xPixels[pixelIndex] = aMember->XLoc;
            yPixels[pixelIndex] = aMember->YLoc;
        }
    }

    clustMember->put(clustIndicies, pixelCnt);
    xLoc->put(xPixels, pixelCnt);
    yLoc->put(yPixels, pixelCnt);

    lats->put(radarInfo.latVals, radarInfo.dataEdges[1]);
    lons->put(radarInfo.lonVals, radarInfo.dataEdges[2]);

    clusterFile.close();

    delete [] clustIndicies;
    delete [] xPixels;
    delete [] yPixels;

    return(true);
}

void PolygonManager::output(const string &fileName)
{
    // -------------------------------------------------------------------------
    NcFile file(fileName.c_str(), NcFile::Replace);
    if (!file.is_valid()) {
        string message = "Failed to open file "+fileName+".";
        REPORT_ERROR(message.c_str());
    }
    // -------------------------------------------------------------------------
    // time information
    if (TimeManager::onLine()) {
        file.add_att("time", TimeManager::getSeconds());
        file.add_att("time_step", TimeManager::getTimeStep());
        file.add_att("steps", TimeManager::getSteps());
    }
    // -------------------------------------------------------------------------
    // dimensions
    NcDim *numVertexDim = file.add_dim("num_total_vertex", vertices.size());
    NcDim *numEdgeDim = file.add_dim("num_total_edge", edges.size());
    NcDim *numPolygonDim = file.add_dim("num_total_polygon", polygons.size());
    // -------------------------------------------------------------------------
    // vertices part
    NcVar *oldVtxLonVar = file.add_var("old_vertex_lon", ncDouble, numVertexDim);
    NcVar *oldVtxLatVar = file.add_var("old_vertex_lat", ncDouble, numVertexDim);
    NcVar *newVtxLonVar = file.add_var("new_vertex_lon", ncDouble, numVertexDim);
    NcVar *newVtxLatVar = file.add_var("new_vertex_lat", ncDouble, numVertexDim);
    oldVtxLonVar->add_att("long_name", "old vertex longitude");
    oldVtxLonVar->add_att("units", "degree_east");
    oldVtxLatVar->add_att("long_name", "old vertex latitude");
    oldVtxLatVar->add_att("units", "degree_north");
    newVtxLonVar->add_att("long_name", "new vertex longitude");
    newVtxLonVar->add_att("units", "degree_east");
    newVtxLatVar->add_att("long_name", "new vertex latitude");
    newVtxLatVar->add_att("units", "degree_north");
    double *oldVtxLon = new double[vertices.size()];
    double *oldVtxLat = new double[vertices.size()];
    double *newVtxLon = new double[vertices.size()];
    double *newVtxLat = new double[vertices.size()];
    Vertex *vertex = vertices.front();
    for (int i = 0; i < vertices.size(); ++i) {
        oldVtxLon[i] = vertex->getCoordinate(OldTimeLevel).getLon()*Rad2Deg;
        oldVtxLat[i] = vertex->getCoordinate(OldTimeLevel).getLat()*Rad2Deg;
        newVtxLon[i] = vertex->getCoordinate(NewTimeLevel).getLon()*Rad2Deg;
        newVtxLat[i] = vertex->getCoordinate(NewTimeLevel).getLat()*Rad2Deg;
        vertex = vertex->next;
    }
    oldVtxLonVar->put(oldVtxLon, vertices.size());
    oldVtxLatVar->put(oldVtxLat, vertices.size());
    newVtxLonVar->put(newVtxLon, vertices.size());
    newVtxLatVar->put(newVtxLat, vertices.size());
    delete [] oldVtxLon;
    delete [] oldVtxLat;
    delete [] newVtxLon;
    delete [] newVtxLat;
    // -------------------------------------------------------------------------
    // edges part
    NcVar *firstPointVar = file.add_var("first_point_idx", ncInt, numEdgeDim);
    NcVar *secondPointVar = file.add_var("second_point_idx", ncInt, numEdgeDim);
    firstPointVar->add_att("long_name", "first point index of edge");
    secondPointVar->add_att("long_name", "second point index of edge");
    int *idx1 = new int[edges.size()];
    int *idx2 = new int[edges.size()];
    Edge *edge = edges.front();
    for (int i = 0; i < edges.size(); ++i) {
        idx1[i] = edge->getEndPoint(FirstPoint)->getID();
        idx2[i] = edge->getEndPoint(SecondPoint)->getID();
        edge = edge->next;
    }
    firstPointVar->put(idx1, edges.size());
    secondPointVar->put(idx2, edges.size());
    delete [] idx1;
    delete [] idx2;
    // test point
    NcVar *oldTestLonVar = file.add_var("old_testpoint_lon", ncDouble, numEdgeDim);
    NcVar *oldTestLatVar = file.add_var("old_testpoint_lat", ncDouble, numEdgeDim);
    NcVar *newTestLonVar = file.add_var("new_testpoint_lon", ncDouble, numEdgeDim);
    NcVar *newTestLatVar = file.add_var("new_testpoint_lat", ncDouble, numEdgeDim);
    oldTestLonVar->add_att("long_name", "old test point longitude");
    oldTestLonVar->add_att("units", "degree_east");
    oldTestLatVar->add_att("long_name", "old test point latitude");
    oldTestLatVar->add_att("units", "degree_north");
    newTestLonVar->add_att("long_name", "new test point longitude");
    newTestLonVar->add_att("units", "degree_east");
    newTestLatVar->add_att("long_name", "new test point latitude");
    newTestLatVar->add_att("units", "degree_north");
#if defined TTS_ONLINE || PREPROCESS
    double *oldTestLon = new double[edges.size()];
    double *oldTestLat = new double[edges.size()];
    double *newTestLon = new double[edges.size()];
    double *newTestLat = new double[edges.size()];
    edge = edges.front();
    for (int i = 0; i < edges.size(); ++i) {
        Vertex *testPoint = edge->getTestPoint();
        oldTestLon[i] = testPoint->getCoordinate(OldTimeLevel).getLon()*Rad2Deg;
        oldTestLat[i] = testPoint->getCoordinate(OldTimeLevel).getLat()*Rad2Deg;
        newTestLon[i] = testPoint->getCoordinate(NewTimeLevel).getLon()*Rad2Deg;
        newTestLat[i] = testPoint->getCoordinate(NewTimeLevel).getLat()*Rad2Deg;
        edge = edge->next;
    }
    oldTestLonVar->put(oldTestLon, edges.size());
//.........这里部分代码省略.........

//
// Creates NetCDF product
//
bool NetCDFProduct(MSG_header *PRO_head, MSG_data* PRO_data,
                   int totalsegs, int *segsindexes,
		   MSG_header *header, MSG_data *msgdat)
{
  struct tm *tmtime;
  char NcName[1024];
  char reftime[64];
  char projname[16];
  int wd, hg;
  int bpp;
  int ncal;
  float *cal;
  NcVar *ivar;
  NcVar *tvar;
  NcDim *tdim;
  NcDim *ldim;
  NcDim *cdim;
  NcDim *caldim;

  int npix = header[0].image_structure->number_of_columns;
  int nlin = header[0].image_structure->number_of_lines;
  size_t npixperseg = npix*nlin;

  size_t total_size = totalsegs*npixperseg;
  MSG_SAMPLE *pixels = new MSG_SAMPLE[total_size];
  memset(pixels, 0, total_size*sizeof(MSG_SAMPLE));
  size_t pos = 0;
  for (int i = 0; i < totalsegs; i ++)
  {
    if (segsindexes[i] >= 0)
      memcpy(pixels+pos, msgdat[segsindexes[i]].image->data,
             npixperseg*sizeof(MSG_SAMPLE));
    pos += npixperseg;
  }

  nlin = nlin*totalsegs;

  // Manage subarea
  if (is_subarea)
  {
    if (AreaLinStart < 0                             ||
        AreaLinStart > nlin - AreaNlin ||
        AreaNlin > nlin - AreaLinStart)
    {
      std::cerr << "Wrong Subarea in lines...." << std::endl;
      throw;
    }
    if (AreaPixStart < 0               ||
        AreaPixStart > npix - AreaNpix ||
        AreaNpix > npix - AreaPixStart)
    {
      std::cerr << "Wrong Subarea in Pixels...." << std::endl;
      throw;
    }
    size_t newsize = AreaNpix * AreaNlin;
    MSG_SAMPLE *newpix = new MSG_SAMPLE[newsize];
    memset(newpix, 0, newsize*sizeof(MSG_SAMPLE));
    for (int i = 0; i < AreaNlin; i ++)
      memcpy(newpix + i * AreaNpix,
             pixels + (AreaLinStart + i) * npix + AreaPixStart,
             AreaNpix * sizeof(MSG_SAMPLE));
    delete [ ] pixels;
    pixels = newpix;
    total_size = newsize;
  }
  else
  {
    AreaNpix = npix;
    AreaNlin = nlin;
  }

  tmtime = PRO_data->prologue->image_acquisition.PlannedAquisitionTime.TrueRepeatCycleStart.get_timestruct( );
  t_enum_MSG_spacecraft spc = header[0].segment_id->spacecraft_id;
  uint_1 chn = header[0].segment_id->spectral_channel_id;
  float sublon = header[0].image_navigation->subsatellite_longitude;
  int cfac = header[0].image_navigation->column_scaling_factor;
  int lfac = header[0].image_navigation->line_scaling_factor;
  int coff = header[0].image_navigation->column_offset;
  int loff = header[0].image_navigation->line_offset;
  float sh = header[0].image_navigation->satellite_h;

  char *channelstring = strdup(MSG_channel_name(spc, chn).c_str( ));
  char *channel = chname(channelstring, strlen(channelstring) + 1);

  // Build up output NetCDF file name and open it
  sprintf( NcName, "%s_%4d%02d%02d_%02d%02d.nc", channel,
           tmtime->tm_year + 1900, tmtime->tm_mon + 1, tmtime->tm_mday,
	   tmtime->tm_hour, tmtime->tm_min );
  NcFile ncf ( NcName , NcFile::Replace );
  if (! ncf.is_valid()) return false;

  // Fill arrays on creation
  ncf.set_fill(NcFile::Fill);

  // Add Global Attributes
  if (! ncf.add_att("Satellite", MSG_spacecraft_name(spc).c_str()))
	            return false;
//.........这里部分代码省略.........

void TracerManager::output(const string &fileName)
{
    // -------------------------------------------------------------------------
    // output polygon stuffs
    polygonManager.output(fileName);
#ifdef TTS_REMAP
    // -------------------------------------------------------------------------
    NcFile file(fileName.c_str(), NcFile::Write);
    if (!file.is_valid()) {
        Message message;
        message << "Failed to open tracer output file \"";
        message << fileName << "\" for appending meshed density field!";
        REPORT_ERROR(message.str());
    }
    // -------------------------------------------------------------------------
    // output tracer densities on the polygons
    NcDim *numPolygonDim = file.get_dim("num_total_polygon");
    double q0[polygonManager.polygons.size()];
    for (int l = 0; l < tracerNames.size(); ++l) {
        char varName[30];
        sprintf(varName, "q%d", l);
        NcVar *qVar = file.add_var(varName, ncDouble, numPolygonDim);
        Polygon *polygon = polygonManager.polygons.front();
        for (int i = 0; i < polygonManager.polygons.size(); ++i) {
            q0[i] = polygon->tracers[l].getDensity();
            polygon = polygon->next;
        }
        qVar->put(q0, polygonManager.polygons.size());
    }
    // -------------------------------------------------------------------------
    // output tracer densities on the mesh
    const RLLMesh &mesh = tracerDensities[0].getMesh();
    int numLon = mesh.getNumLon()-2;
    int numLat = mesh.getNumLat();
    double lon[numLon], lat[numLat];
    for (int i = 0; i < numLon; ++i)
        lon[i] = mesh.lon(i+1)*Rad2Deg;
    for (int j = 0; j < numLat; ++j)
        lat[j] = mesh.lat(j)*Rad2Deg;
    NcDim *lonDim = file.add_dim("lon", numLon);
    NcDim *latDim = file.add_dim("lat", numLat);
    NcVar *lonVar = file.add_var("lon", ncDouble, lonDim);
    lonVar->add_att("long_name", "longitude");
    lonVar->add_att("units", "degrees_east");
    lonVar->put(lon, numLon);
    NcVar *latVar = file.add_var("lat", ncDouble, latDim);
    latVar->add_att("long_name", "latitude");
    latVar->add_att("units", "degrees_north");
    latVar->put(lat, numLat);
    NcVar *areaVar = file.add_var("area_mesh", ncDouble, latDim, lonDim);
    areaVar->add_att("long_name", "area of fixed mesh cell");
    areaVar->add_att("units", "m2");
    double area[numLat][numLon];
    for (int i = 0; i < numLon; ++i)
        for (int j = 0; j < numLat; ++j)
            area[j][i] = tracerDensities[0].getMesh(Field::Bound).area(i, j);
    areaVar->put(&area[0][0], numLat, numLon);
    double q[numLat][numLon];
    for (int l = 0; l < tracerNames.size(); ++l) {
        char varName[30];
        sprintf(varName, "q%d_mesh", l);
        NcVar *qVar = file.add_var(varName, ncDouble, latDim, lonDim);
        qVar->add_att("long_name", tracerNames[l].c_str());
        for (int i = 0; i < numLon; ++i)
            for (int j = 0; j < numLat; ++j)
                q[j][i] = tracerDensities[l].values(i+1, j, 0).getNew();
        qVar->put(&q[0][0], numLat, numLon);
    }
    // -------------------------------------------------------------------------
    file.close();
#endif
    NOTICE("TracerManager", fileName+" is generated.");
}

bool OutputManagerReference::OpenFile(
	const std::string & strFileName
) {
#ifdef TEMPEST_NETCDF
	// Determine processor rank; only proceed if root node
	int nRank = 0;

#ifdef TEMPEST_MPIOMP
	MPI_Comm_rank(MPI_COMM_WORLD, &nRank);
#endif

	// The active model
	const Model & model = m_grid.GetModel();

	// Open NetCDF file on root process
	if (nRank == 0) {

		// Check for existing NetCDF file
		if (m_pActiveNcOutput != NULL) {
			_EXCEPTIONT("NetCDF file already open");
		}

		// Append .nc extension to file
		std::string strNcFileName = strFileName + ".nc";

		// Open new NetCDF file
		m_pActiveNcOutput = new NcFile(strNcFileName.c_str(), NcFile::Replace);
		if (m_pActiveNcOutput == NULL) {
			_EXCEPTION1("Error opening NetCDF file \"%s\"",
				strNcFileName.c_str());
		}
		if (!m_pActiveNcOutput->is_valid()) {
			_EXCEPTION1("Error opening NetCDF file \"%s\"",
				strNcFileName.c_str());
		}

		// Create nodal time dimension
		NcDim * dimTime = m_pActiveNcOutput->add_dim("time");
		if (dimTime == NULL) {
			_EXCEPTIONT("Error creating \"time\" dimension");
		}

		m_varTime = m_pActiveNcOutput->add_var("time", ncDouble, dimTime);
		if (m_varTime == NULL) {
			_EXCEPTIONT("Error creating \"time\" variable");
		}

		std::string strUnits =
			"days since " + model.GetStartTime().ToDateString();

		std::string strCalendarName = model.GetStartTime().GetCalendarName();

		m_varTime->add_att("long_name", "time");
		m_varTime->add_att("units", strUnits.c_str());
		m_varTime->add_att("calendar", strCalendarName.c_str());
		m_varTime->add_att("bounds", "time_bnds");

		// Create levels dimension
		NcDim * dimLev =
			m_pActiveNcOutput->add_dim("lev", m_dREtaCoord.GetRows());

		// Create interfaces dimension
		NcDim * dimILev =
			m_pActiveNcOutput->add_dim("ilev", m_grid.GetRElements()+1);

		// Create latitude dimension
		NcDim * dimLat =
			m_pActiveNcOutput->add_dim("lat", m_nYReference);

		// Create longitude dimension
		NcDim * dimLon =
			m_pActiveNcOutput->add_dim("lon", m_nXReference);

		// Output physical constants
		const PhysicalConstants & phys = model.GetPhysicalConstants();

		m_pActiveNcOutput->add_att("earth_radius", phys.GetEarthRadius());
		m_pActiveNcOutput->add_att("g", phys.GetG());
		m_pActiveNcOutput->add_att("omega", phys.GetOmega());
		m_pActiveNcOutput->add_att("alpha", phys.GetAlpha());
		m_pActiveNcOutput->add_att("Rd", phys.GetR());
		m_pActiveNcOutput->add_att("Cp", phys.GetCp());
		m_pActiveNcOutput->add_att("T0", phys.GetT0());
		m_pActiveNcOutput->add_att("P0", phys.GetP0());
		m_pActiveNcOutput->add_att("rho_water", phys.GetRhoWater());
		m_pActiveNcOutput->add_att("Rvap", phys.GetRvap());
		m_pActiveNcOutput->add_att("Mvap", phys.GetMvap());
		m_pActiveNcOutput->add_att("Lvap", phys.GetLvap());

		// Output grid parameters
		m_pActiveNcOutput->add_att("Ztop", m_grid.GetZtop());

		// Output equation set
		const EquationSet & eqn = model.GetEquationSet();

		m_pActiveNcOutput->add_att("equation_set", eqn.GetName().c_str());

		// Create variables
		for (int c = 0; c < eqn.GetComponents(); c++) {
			if ((m_fOutputAllVarsOnNodes) ||
//.........这里部分代码省略.........

//.........这里部分代码省略.........
					}
					fWhitespace = false;
				}
			}

			// Latitude and longitude index
			int iLon =
				static_cast<int>(static_cast<double>(nLon)
					* (dLon - dLonBegin) / (dLonEnd - dLonBegin));
			int iLat =
				static_cast<int>(static_cast<double>(nLat)
					* (dLat - dLatBegin) / (dLatEnd - dLatBegin));

			if (iLon == (-1)) {
				iLon = 0;
			}
			if (iLon == nLon) {
				iLon = nLon - 1;
			}
			if (iLat == (-1)) {
				iLat = 0;
			}
			if (iLat == nLat) {
				iLat = nLat - 1;
			}

			if ((iLat < 0) || (iLat >= nLat)) {
				_EXCEPTION1("Latitude index (%i) out of range", iLat);
			}
			if ((iLon < 0) || (iLon >= nLon)) {
				_EXCEPTION1("Longitude index (%i) out of range", iLon);
			}

			nCounts[iLat][iLon]++;
		}

		fclose(fp);
	}

	AnnounceEndBlock("Done");

	// Output results
	AnnounceStartBlock("Output results");

	// Load the netcdf output file
	NcFile ncOutput(strOutputFile.c_str(), NcFile::Replace);
	if (!ncOutput.is_valid()) {
		_EXCEPTION1("Unable to open output file \"%s\"",
			strOutputFile.c_str());
	}

	// Create output
	NcDim * dimLat = ncOutput.add_dim("lat", nLat);
	NcDim * dimLon = ncOutput.add_dim("lon", nLon);

	NcVar * varLat = ncOutput.add_var("lat", ncDouble, dimLat);
	NcVar * varLon = ncOutput.add_var("lon", ncDouble, dimLon);

	varLat->add_att("units", "degrees_north");
	varLon->add_att("units", "degrees_east");

	DataVector<double> dLat(nLat);
	DataVector<double> dLon(nLon);

	for (int j = 0; j < nLat; j++) {
		dLat[j] = dLatBegin
			+ (dLatEnd - dLatBegin)
				* (static_cast<double>(j) + 0.5)
				/ static_cast<double>(nLat);
	}
	for (int i = 0; i < nLon; i++) {
		dLon[i] = dLonBegin
			+ (dLonEnd - dLonBegin)
			* (static_cast<double>(i) + 0.5)
			/ static_cast<double>(nLon);
	}

	varLat->put(&(dLat[0]), nLat);
	varLon->put(&(dLon[0]), nLon);

	// Output counts
	NcVar * varCount =
		ncOutput.add_var(
			strOutputVariable.c_str(),
			ncInt,
			dimLat,
			dimLon);

	varCount->put(&(nCounts[0][0]), nLat, nLon);

	ncOutput.close();

	AnnounceEndBlock("Done");

	AnnounceBanner();

} catch(Exception & e) {
	Announce(e.ToString().c_str());
}
}

void gen(const char* path, NcFile::FileFormat format)		// Generate a netCDF file
{

    NcFile nc(path, NcFile::Replace, NULL, 0, format); // Create, leave in define mode

    // Check if the file was opened successfully
    if (! nc.is_valid()) {
	cerr << "can't create netCDF file " << path << "\n";
	return;
    }

    // Create dimensions
    const int NLATS = 4;
    const int NLONS = 3;
    const int NFRTIMES = 2;
    const int TIMESTRINGLEN = 20;
    NcDim* latd = nc.add_dim("lat", NLATS);
    NcDim* lond = nc.add_dim("lon", NLONS);
    NcDim* frtimed = nc.add_dim("frtime"); // unlimited dimension
    NcDim* timelend = nc.add_dim("timelen", TIMESTRINGLEN); 

    // Create variables and their attributes
    NcVar* P = nc.add_var("P", ncFloat, frtimed, latd, lond);
    P->add_att("long_name", "pressure at maximum wind");
    P->add_att("units", "hectopascals");
    static float range[] = {0., 1500.};
    P->add_att("valid_range", 2, range);
    P->add_att("_FillValue", -9999.0f);

    NcVar* lat = nc.add_var("lat", ncFloat, latd);
    lat->add_att("long_name", "latitude");
    lat->add_att("units", "degrees_north");

    NcVar* lon = nc.add_var("lon", ncFloat, lond);
    lon->add_att("long_name", "longitude");
    lon->add_att("units", "degrees_east");

    NcVar* frtime = nc.add_var("frtime", ncLong, frtimed);
    frtime->add_att("long_name", "forecast time");
    frtime->add_att("units", "hours");

    NcVar* reftime = nc.add_var("reftime",ncChar,timelend);
    reftime->add_att("long_name", "reference time");
    reftime->add_att("units", "text_time");

    NcVar* scalar = nc.add_var("scalarv", ncInt);
    scalar->add_att("scalar_att", 1);

    // Global attributes
    nc.add_att("history", "created by Unidata LDM from NPS broadcast");
    nc.add_att("title", "NMC Global Product Set: Pressure at Maximum Wind");

    // Start writing data, implictly leaves define mode

    static float lats[NLATS] = {-90, -87.5, -85, -82.5};
    lat->put(lats, NLATS);

    static float lons[NLONS] = {-180, -175, -170};
    lon->put(lons, NLONS);

    static int frtimes[NFRTIMES] = {12, 18};
    frtime->put(frtimes, NFRTIMES);

    static const char* s = "1992-3-21 12:00" ;
    reftime->put(s, strlen(s));

    static float P_data[2][4][3] = {
	{{950, 951, 952}, {953, 954, 955}, {956, 957, 958}, {959, 960, 961}},
	{{962, 963, 964}, {965, 966, 967}, {968, 969, 970}, {971, 972, 973}}
      };
    // We could write all P data at once with P->put(&P_data[0][0][0], P->edges()),
    // but instead we write one record at a time, to show use of setcur().
    long rec = 0;                                      // start at zero-th
    const long nrecs = 1;		               // # records to write
    P->put(&P_data[0][0][0], nrecs, NLATS, NLONS);           // write zero-th record
    P->set_cur(++rec);		                       // set to next record
    P->put(&P_data[1][0][0], nrecs, NLATS, NLONS); // write next record

    // close of nc takes place in destructor
}