<p><b>cahrens@lanl.gov</b> 2010-09-21 19:01:10 -0600 (Tue, 21 Sep 2010)</p><p>New version of translator that creates either sphere or latlon view and either single or multilayer view. When compiled, it is called mpasnc2vtu. Run with no arguments to see the usage information. Writes out data in XML and binary.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/graphics/paraview/translator/Makefile
===================================================================
--- branches/ocean_projects/graphics/paraview/translator/Makefile 2010-09-17 17:59:55 UTC (rev 506)
+++ branches/ocean_projects/graphics/paraview/translator/Makefile 2010-09-22 01:01:10 UTC (rev 507)
@@ -10,8 +10,8 @@
# -lstdc++ \
# -L/opt/Intel/cce/10.0.023/lib/ -lirc
-transdual: nc2vtk_dual.cpp
- $(CC) $(INCLUDES) -o $@ $< $(LIBS)
+mpasnc2vtu: mpasnc2vtu.cpp EncodeBase64.o
+ $(CC) $(INCLUDES) -o $@ $< EncodeBase64.o $(LIBS)
.cpp.o:
$(CC) $(MYDEBUG) $(INCLUDES) -c $<
@@ -22,5 +22,8 @@
transdualbin: nc2vtu_dualbin.cpp EncodeBase64.o
$(CC) $(INCLUDES) -o $@ $< EncodeBase64.o $(LIBS)
+transdual: nc2vtu_dual.cpp
+ $(CC) $(INCLUDES) -o $@ $< $(LIBS)
+
clean:
- rm transdual transdualbin *.o
+ rm transdual transdualbin mpasnc2vtu *.o
Added: branches/ocean_projects/graphics/paraview/translator/mpasnc2vtu.cpp
===================================================================
--- branches/ocean_projects/graphics/paraview/translator/mpasnc2vtu.cpp (rev 0)
+++ branches/ocean_projects/graphics/paraview/translator/mpasnc2vtu.cpp 2010-09-22 01:01:10 UTC (rev 507)
@@ -0,0 +1,1645 @@
+////////////////////////////////////////////////////////////////////////////////
+// mpasnc2vtu.cpp
+//
+// This program translates a MPAS netCDF data file to a dual grid in ParaView
+// .vtu format. Can do either sphere view or lat/lon view. Can show single
+// or multi-layers. Can translate to xml with either ascii or binary data.
+//
+// Assumes little-endian platform.
+// Assume all variables are of interest.
+// Assume variable data type is double.
+// Assume variable dims are (Time, nCells|nVertices, nVertLevels|nVertLevelsP1)
+// Assume no more than 100 vars each for cell and point data
+// Does not deal with edge data.
+// Only vis up to nVertLevels, not nVertLevelsP1.
+// Doubles converted to be valid follow these rules:
+// if a NaN, then become -DBL_MAX.
+// if positive infinity, then become DBL_MAX
+// if negative infinity, then become -DBL_MAX
+//
+// Version 1.5
+// Christine Ahrens
+// 9/21/2010
+////////////////////////////////////////////////////////////////////////////////
+
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include "stdlib.h"
+#include "vtkCellType.h"
+#include "netcdfcpp.h"
+#include <string>
+#include <cmath>
+//#include <math>
+#include <cfloat>
+#include "EncodeBase64.h"
+
+using namespace std;
+
+struct params {
+ bool projectLatLon;
+ bool writeBinary;
+ bool singleLayer;
+ bool isAtmosphere;
+ bool includeTopography;
+ bool doBugFix;
+ int outputVertLevel;
+ char* inputFileName;
+ char* outputFileName;
+ double layerThickness;
+ bool isZeroCentered;
+ double centerLon;
+ double centerRad;
+ int startTimeStep;
+ int finishTimeStep;
+};
+
+// points and cells refer to DUAL points and cells
+struct geometry {
+ int nVertLevels;
+ int maxNVertLevels;
+ int Time;
+ int numCells;
+ int numPoints;
+ int cellOffset;
+ int pointOffset;
+ int pointsPerCell;
+ int currentExtraPoint; // current extra point
+ int currentExtraCell; // current extra cell
+ double* pointX; // x coord of point
+ double* pointY; // y coord of point
+ double* pointZ; // z coord of point
+ int modNumPoints;
+ int modNumCells;
+ int* origConnections; // original connections
+ int* modConnections; // modified connections
+ int* cellMap; // maps from added cell to original cell #
+ int* pointMap; // maps from added point to original point #
+ int* maxLevelPoint; //
+ int maxCells; // max cells
+ int maxPoints; // max points
+ int verticalIndex; // for singleLayer, which vertical level
+};
+
+#define MAX_VARS 100
+struct vars {
+ NcVar* cellVars[MAX_VARS];
+ NcVar* pointVars[MAX_VARS];
+ int numCellVars;
+ int numPointVars;
+};
+
+#define DEFAULT_LAYER_THICKNESS 10
+double PI = 3.141592;
+#define BUFF_SIZE 2048
+
+#define CHECK_MALLOC(ptr) \
+ if (ptr == NULL) { \
+ cerr << "malloc failed!</font>
<font color="blue">"; \
+ exit(1); \
+ }
+
+
+
+// append inFile to outFile
+int AppendFile(string *inFile, string *outFile)
+{
+ char buff[BUFF_SIZE];
+ int readBytes = 1;
+
+ ifstream inFileStream(inFile->c_str(), ios::in|ios::binary);
+ if(!inFileStream)
+ {
+ return -1;
+ }
+
+ ofstream outFileStream(outFile->c_str(), ios::out|ios::app|ios::binary);
+ if(!outFileStream)
+ {
+ return -1;
+ }
+
+ while(readBytes != 0)
+ {
+ inFileStream.read((char*)buff, BUFF_SIZE);
+ readBytes = inFileStream.gcount();
+ outFileStream.write((char*)buff, readBytes);
+ }
+ return 0;
+}
+
+int IsInf_f(float x){
+ if (isinf(x)) {
+ return (x < 0.0 ? -1 : 1);
+ } else return 0;
+}
+
+unsigned int PadSize(unsigned int binarySize) {
+ switch ((binarySize*4)%3) {
+ case 1:
+ binarySize = binarySize + 2;
+ break;
+ case 2:
+ binarySize++;
+ break;
+ case 0:
+ break;
+ }
+ return binarySize;
+}
+
+float ConvertDouble2ValidFloat(double inputData) {
+
+ // check for NaN
+ if (inputData != inputData) {
+ //cerr << "found NaN!" << endl;
+ return -FLT_MAX;
+ }
+
+ // check for infinity
+ int retval = IsInf_f((float)inputData);
+ if (retval < 0) {
+ return -FLT_MAX;
+ } else if (retval > 0) {
+ return FLT_MAX;
+ }
+
+ // check number too small for float
+ if (abs(inputData) < 1e-126) return 0.0;
+
+ if ((float)inputData == 0) return 0.0;
+
+ if ((abs(inputData) > 0) && (abs(inputData) < FLT_MIN)) {
+ if (inputData < 0) return -FLT_MIN; else return FLT_MIN;
+ }
+
+ if (abs(inputData) > FLT_MAX) {
+ if (inputData < 0) return -FLT_MAX; else return FLT_MAX;
+ }
+
+ return (float)inputData;
+}
+
+double ConvertDouble2ValidDouble(double inputData) {
+
+ // check for NaN
+ if (isnan(inputData)) {
+ cerr << "found NaN!" << endl;
+ return -DBL_MAX;
+ }
+
+ // check for infinity
+ int retval = isinf(inputData);
+ if (retval < 0) {
+ return -DBL_MAX;
+ } else if (retval > 0) {
+ return DBL_MAX;
+ }
+
+ return inputData;
+}
+
+int CartesianToSpherical(double x, double y, double z, double* rho,
+ double* phi, double* theta)
+{
+ double trho, ttheta, tphi;
+
+ trho = sqrt( (x*x) + (y*y) + (z*z));
+ ttheta = atan2(y, x);
+ tphi = acos(z/(trho));
+ if (isnan(trho) || isnan(ttheta) || isnan(tphi)) {
+ cerr << "CartesianToSpherical on x: " << x << " y: " << y << " z: " << z << " produced nan!" << endl;
+ return -1;
+ }
+ *rho = trho;
+ *theta = ttheta;
+ *phi = tphi;
+ return 0;
+}
+
+int SphericalToCartesian(double rho, double phi, double theta, double* x,
+ double* y, double* z)
+{
+ double tx, ty, tz;
+ tx = rho* sin(phi) * cos(theta);
+ ty = rho* sin(phi) * sin(theta);
+ tz = rho* cos(phi);
+ if (isnan(tx) || isnan(ty) || isnan(tz)) {
+ cerr << "SphericalToCartesian produced nan!" << endl;
+ return -1;
+ }
+
+ *x = tx;
+ *y = ty;
+ *z = tz;
+ return 0;
+}
+
+void WriteDouble(EncodeBase64 &b64, ofstream &file, double data, bool writeBinary) {
+ double validData = ConvertDouble2ValidDouble (data);
+ if (writeBinary) {
+ b64.Write((const unsigned char*)&validData, sizeof(validData));
+ } else {
+ file << validData << " ";
+ }
+}
+
+void WriteFloat(EncodeBase64 &b64, ofstream &file, double data, bool writeBinary) {
+ float validData = ConvertDouble2ValidFloat (data);
+ if (writeBinary) {
+ b64.Write((const unsigned char*)&validData, sizeof(validData));
+ } else {
+ file << validData << " ";
+ }
+}
+
+void WriteInt(EncodeBase64 &b64, ofstream &file, int data, bool writeBinary) {
+ if (writeBinary) {
+ b64.Write((const unsigned char*)&data, sizeof(data));
+ } else {
+ file << data << " ";
+ }
+}
+
+void WriteChar(EncodeBase64 &b64, ofstream &file, char data, bool writeBinary) {
+ if (writeBinary) {
+ b64.Write((const unsigned char*)&data, sizeof(data));
+ } else {
+ file << (int)data << " ";
+ }
+}
+
+void PrintUsage() {
+ cerr << "USAGE: mpasnc2vtu [OPTIONS] infile.nc outfile.vtu" << endl;
+ cerr << "Create view of MPAS-style data from input file infile.nc " << endl;
+ cerr << "and write to XML binary .vtu format for use in ParaView." << endl;
+ cerr << "A series of vtu files are created, one file for each time step." << endl;
+ cerr << "with time prepended to the file name (e.g. 0outfile.vtu)." << endl;
+ cerr << "OPTIONS:" << endl;
+ cerr << " -l intval " << endl;
+ cerr << " vertical level, in range 0 to nVertLevels-1, default is 0 " << endl;
+ cerr << " -m " << endl;
+ cerr << " multilayer view -- default is single layer view" << endl;
+ cerr << " -p " << endl;
+ cerr << " project to lat/lon -- default is sphere view" << endl;
+ cerr << " -t floatval " << endl;
+ cerr << " layer_thickness -- use negative value for atmosphere, default is 10" << endl;
+ cerr << " -c floatval " << endl;
+ cerr << " display image with center_longitude in degrees, default is 180" << endl;
+ cerr << " -z " << endl;
+ cerr << " input grid has center longitude at 0, the default is 180 degree or PI radians" << endl;
+ cerr << " -a " << endl;
+ cerr << " atmosphere data (longitude data is in range -PI to PI) " << endl;
+ cerr << " The -a flag also converts the multilayer view to have the " << endl;
+ cerr << " vertical levels stack in the positive Z direction, rather " << endl;
+ cerr << " the negative Z direction, which is the default." << endl;
+ cerr << " -s " << endl;
+ cerr << " timestep to start on (defaults to 0) " << endl;
+ cerr << " -f " << endl;
+ cerr << " timestep to finish on (defaults to last timestep available) " << endl;
+ cerr << " -g " << endl;
+ cerr << " include topography (assumes you have maxLevelPoint variable) " << endl;
+ cerr << " -b " << endl;
+ cerr << " do bugfix to ensure no connections to faraway points " << endl;
+ cerr << " -x " << endl;
+ cerr << " write data in ascii, not binary " << endl;
+
+ cerr << "Variables with time and vertical level are written out." << endl;
+ cerr << "[vtk datafile Version 2.0, projlatlonbin Version 1.3]" << endl;
+ exit(1);
+}
+
+void GetArgs(int argc, char* argv[], struct params &p)
+{
+ int c;
+ optind = 1;
+
+ //set default params
+ p.writeBinary = true;
+ p.singleLayer = true;
+ p.isAtmosphere = false;
+ p.includeTopography = false;
+ p.doBugFix = false;
+ p.outputVertLevel = 0;
+ p.inputFileName = NULL;
+ p.outputFileName = NULL;
+ p.layerThickness = DEFAULT_LAYER_THICKNESS;
+ p.isZeroCentered = false;
+ p.centerLon = 180.0;
+ p.centerRad = p.centerLon * PI / 180.0;
+ p.startTimeStep = 0;
+ p.finishTimeStep = -1;
+ p.projectLatLon = false;
+
+ bool lflag = false;
+ bool tflag = false;
+ bool mflag = false;
+
+ // -l vertical level (if not specified, assume multilayer)
+ // -t layer thickness (defaults to 1/10 * radius/#layers)
+ // -c center longitude in degrees for latlon projection
+ while ((c=getopt(argc, argv, "l:t:c:mazpxgbs:f:")) != -1) {
+ switch(c) {
+ case 'l':
+ lflag = true;
+ if (tflag) {
+ cerr << "Can't specify -l and -t at same time." << endl;
+ exit(1);
+ }
+ p.outputVertLevel = atoi(optarg);
+ p.singleLayer = true;
+ break;
+ case 't':
+ tflag = true;
+ if (lflag) {
+ cerr << "Can't specify -l and -t at same time." << endl;
+ exit(1);
+ }
+ p.layerThickness = atof(optarg);
+ p.singleLayer = false;
+ break;
+ case 'c':
+ p.centerLon = (double)atof(optarg);
+ p.centerRad = p.centerLon * PI / 180.0;
+ break;
+ case 'm':
+ mflag = true;
+ if (lflag) {
+ cerr << "Can't specify -l and -m at same time." << endl;
+ exit(1);
+ }
+ p.singleLayer = false;
+ break;
+ case 'a':
+ p.isAtmosphere = true;
+ break;
+ case 'x':
+ p.writeBinary = false;
+ break;
+ case 's':
+ p.startTimeStep = atoi(optarg);
+ break;
+ case 'f':
+ p.finishTimeStep = atoi(optarg);
+ break;
+ case 'z':
+ p.isZeroCentered = true;
+ break;
+ case 'p':
+ p.projectLatLon = true;
+ break;
+ case 'g':
+ p.includeTopography = true;
+ break;
+ case 'b':
+ p.doBugFix = true;
+ break;
+ default:
+ cerr << "Unrecognized option: -" << c << endl;
+ PrintUsage();
+ exit(1);
+ break;
+ }
+ }
+
+ p.inputFileName = argv[optind++];
+ p.outputFileName = argv[optind];
+
+}
+
+// Check if there is a variable by that name
+bool isNcVar(NcFile &ncFile, NcToken name) {
+ int num_vars = ncFile.num_vars();
+ for (int i = 0; i < num_vars; i++) {
+ NcVar* ncVar = ncFile.get_var(i);
+ if ((strcmp(ncVar->name(), name)) == 0) {
+ // we have a match, so return
+ return true;
+ }
+ }
+ return false;
+}
+
+// Check if there is a dimension by that name
+bool isNcDim(NcFile &ncFile, NcToken name) {
+ int num_dims = ncFile.num_dims();
+ //cerr << "looking for: " << name << endl;
+ for (int i = 0; i < num_dims; i++) {
+ NcDim* ncDim = ncFile.get_dim(i);
+ //cerr << "checking " << ncDim->name() << endl;
+ if ((strcmp(ncDim->name(), name)) == 0) {
+ // we have a match, so return
+ return true;
+ }
+ }
+ return false;
+}
+
+
+#define CHECK_DIM(ncFile, name) \
+ if (!isNcDim(ncFile, name)) { \
+ cerr << "Cannot find dimension: " << name << endl; \
+ exit(1); \
+ }
+
+#define CHECK_VAR(ncFile, name) \
+ if (!isNcVar(ncFile, name)) { \
+ cerr << "Cannot find variable: " << name << endl; \
+ exit(1); \
+ }
+
+
+// put dims into dual grid geometry
+void GetNcDims(NcFile &ncFile, geometry &g)
+{
+ CHECK_DIM(ncFile, "nCells");
+ NcDim* nCells = ncFile.get_dim("nCells");
+ g.numPoints = nCells->size();
+ g.pointOffset = 1;
+
+ CHECK_DIM(ncFile, "nVertices");
+ NcDim* nVertices = ncFile.get_dim("nVertices");
+ g.numCells = nVertices->size();
+ g.cellOffset = 0;
+
+ CHECK_DIM(ncFile, "vertexDegree");
+ NcDim* vertexDegree = ncFile.get_dim("vertexDegree");
+ g.pointsPerCell = vertexDegree->size();
+
+ CHECK_DIM(ncFile, "Time");
+ NcDim* Time = ncFile.get_dim("Time");
+ g.Time = Time->size();
+
+ CHECK_DIM(ncFile, "nVertLevels");
+ NcDim* nVertLevels = ncFile.get_dim("nVertLevels");
+ g.nVertLevels = nVertLevels->size();
+
+ g.maxNVertLevels = g.nVertLevels;
+}
+
+void CheckArgs(params &p, geometry &g)
+{
+
+ if ((g.pointsPerCell != 3) && (g.pointsPerCell != 4)) {
+ cerr << "This code is only for hexagonal or quad grid projection" << endl;
+ exit(1);
+ }
+
+ // double-check we can do multilayer
+ if ((!p.singleLayer) && (g.maxNVertLevels == 1)) {
+ p.singleLayer = true;
+ }
+
+ if (p.singleLayer && (p.outputVertLevel > (g.maxNVertLevels-1))) {
+ cerr << "There is no data for level " << p.outputVertLevel << "." << endl;
+ exit(1);
+ }
+
+ if (p.singleLayer && (p.outputVertLevel < 0)) {
+ cerr << "Level number must be 0 or greater." << endl;
+ exit(1);
+ }
+
+ if ((p.centerLon < 0) || (p.centerLon > 360)) {
+ cerr << "Center longitude must be between 0 and 360." << endl;
+ exit(1);
+ }
+
+ if (p.singleLayer) {
+ g.maxNVertLevels = 1;
+ }
+
+ if (p.isAtmosphere) {
+ p.layerThickness = -p.layerThickness;
+ }
+
+ int availTimeSteps = g.Time;
+
+ if (p.startTimeStep > (availTimeSteps-1)) {
+ cerr << "Start timestep is too large, last timestep is: " << availTimeSteps-1 << endl;
+ exit(1);
+ }
+
+ if (p.startTimeStep < 0) {
+ cerr << "Start timestep is too small, must be at least 0." << endl;
+ exit(1);
+ }
+
+ if (p.finishTimeStep > (availTimeSteps-1)) {
+ cerr << "Finish timestep is too large, last timestep is: " << availTimeSteps-1 << endl;
+ exit(1);
+ }
+
+ if (p.finishTimeStep < 0) {
+ p.finishTimeStep = availTimeSteps-1;
+ }
+
+ cerr << "centerLon: " << p.centerLon << " singleLayer: " << p.singleLayer << " outputVertLevel: " << p.outputVertLevel << " start: " << p.startTimeStep << " finish: " << p.finishTimeStep << " isAtmosphere: " << p.isAtmosphere << " includeTopography: " << p.includeTopography << " doBugFix: " << p.doBugFix << " layerThickness: " << p.layerThickness << endl;
+
+}
+
+
+void GetNcVars (NcFile &ncFile, vars &v, const char* cellDimName,
+ const char* pointDimName)
+{
+ int cellVarIndex = -1;
+ int pointVarIndex = -1;
+
+ int numVars = ncFile.num_vars();
+
+ for (int i = 0; i < numVars; i++) {
+ NcVar* aVar = ncFile.get_var(i);
+
+ // must have 3 dims
+ // (Time, nCells | nVertices, nVertLevels | nVertLevelsP1)
+
+ int numDims = aVar->num_dims();
+ //cout << "Num Dims of var: " << aVar->name() << " is " << numDims << endl;
+ if (numDims != 3) continue;
+
+
+ // TODO, check if it is a double
+ // assume a double for now
+
+ // check for Time dim 0
+ NcToken dim0Name = aVar->get_dim(0)->name();
+ if (strcmp(dim0Name, "Time"))
+ continue;
+
+ // check for dim 1 being cell or point
+ bool isCellData = false;
+ bool isPointData = false;
+ NcToken dim1Name = aVar->get_dim(1)->name();
+ if (!strcmp(dim1Name, cellDimName))
+ isCellData = true;
+ else if (!strcmp(dim1Name, pointDimName))
+ isPointData = true;
+ else continue;
+
+ // check if dim 2 is nVertLevels or nVertLevelsP1, too
+ NcToken dim2Name = aVar->get_dim(2)->name();
+ if ((strcmp(dim2Name, "nVertLevels"))
+ && (strcmp(dim2Name, "nVertLevelsP1"))) {
+ continue;
+ }
+
+ // Add to cell or point var array
+ if (isCellData) { // means it is cell data
+ cellVarIndex++;
+ if (cellVarIndex > MAX_VARS-1) {
+ cerr << "Exceeded number of cell vars." << endl;
+ exit(1);
+ }
+ v.cellVars[cellVarIndex] = aVar;
+ //cout << "Adding var " << aVar->name() << " to cellVars" << endl;
+ } else if (isPointData)
+ {
+ pointVarIndex++;
+ if (pointVarIndex > MAX_VARS-1) {
+ cerr << "Exceeded number of point vars." << endl;
+ exit(1);
+ }
+ v.pointVars[pointVarIndex] = aVar;
+ //cout << "Adding var " << aVar->name() << " to pointVars" << endl;
+ }
+ }
+
+ v.numPointVars = pointVarIndex+1;
+ v.numCellVars = cellVarIndex+1;
+
+}
+
+// allocate into sphere view of dual geometry
+void AllocSphereGeometry(NcFile &ncFile, params &p, geometry &g)
+{
+ CHECK_VAR(ncFile, "xCell");
+ g.pointX = (double*)malloc((g.numPoints + g.pointOffset) * sizeof(double));
+ CHECK_MALLOC(g.pointX);
+ NcVar* xCellVar = ncFile.get_var("xCell");
+ xCellVar->get(g.pointX + g.pointOffset, g.numPoints);
+ // point 0 is 0.0
+ g.pointX[0] = 0.0L;
+
+ CHECK_VAR(ncFile, "yCell");
+ g.pointY = (double*)malloc((g.numPoints + g.pointOffset) * sizeof(double));
+ CHECK_MALLOC(g.pointY);
+ NcVar* yCellVar = ncFile.get_var("yCell");
+ yCellVar->get(g.pointY+g.pointOffset, g.numPoints);
+ // point 0 is 0.0
+ g.pointY[0] = 0.0L;
+
+ CHECK_VAR(ncFile, "zCell");
+ g.pointZ = (double*)malloc((g.numPoints + g.pointOffset) * sizeof(double));
+ CHECK_MALLOC(g.pointZ);
+ NcVar* zCellVar = ncFile.get_var("zCell");
+ zCellVar->get(g.pointZ+g.pointOffset, g.numPoints);
+ // point 0 is 0.0
+ g.pointZ[0] = 0.0L;
+
+ CHECK_VAR(ncFile, "cellsOnVertex");
+ g.origConnections = (int *) malloc(g.numCells * g.pointsPerCell *
+ sizeof(int));
+ CHECK_MALLOC(g.origConnections);
+ NcVar *connectionsVar = ncFile.get_var("cellsOnVertex");
+ connectionsVar->get(g.origConnections, g.numCells, g.pointsPerCell);
+
+ if (p.includeTopography) {
+ CHECK_VAR(ncFile, "maxLevelCell");
+ g.maxLevelPoint = (int*)malloc((g.numPoints + g.numPoints) * sizeof(int));
+ CHECK_MALLOC(g.maxLevelPoint);
+ NcVar *maxLevelPointVar = ncFile.get_var("maxLevelCell");
+ maxLevelPointVar->get(g.maxLevelPoint + g.pointOffset, g.numPoints);
+ }
+
+ g.currentExtraPoint = g.numPoints + g.pointOffset;
+ g.currentExtraCell = g.numCells + g.cellOffset;
+
+ if (p.singleLayer) {
+ g.maxCells = g.currentExtraCell;
+ g.maxPoints = g.currentExtraPoint;
+ } else {
+ g.maxCells = g.currentExtraCell*g.maxNVertLevels;
+ g.maxPoints = g.currentExtraPoint*(g.maxNVertLevels+1);
+ }
+}
+
+// allocate into lat/lon projection of dual geometry
+void AllocLatLonGeometry(NcFile &ncFile, params &p, geometry &g)
+{
+ const float BLOATFACTOR = .5;
+ g.modNumPoints = (int)floor(g.numPoints*(1.0 + BLOATFACTOR));
+ g.modNumCells = (int)floor(g.numCells*(1.0 + BLOATFACTOR))+1;
+
+ CHECK_VAR(ncFile, "lonCell");
+ g.pointX = (double*)malloc(g.modNumPoints * sizeof(double));
+ CHECK_MALLOC(g.pointX);
+ NcVar* xCellVar = ncFile.get_var("lonCell");
+ xCellVar->get(g.pointX + g.pointOffset, g.numPoints);
+ // point 0 is 0.0
+ *g.pointX = 0.0;
+
+ CHECK_VAR(ncFile, "latCell");
+ g.pointY = (double*)malloc(g.modNumPoints * sizeof(double));
+ CHECK_MALLOC(g.pointY);
+ NcVar* yCellVar = ncFile.get_var("latCell");
+ yCellVar->get(g.pointY+g.pointOffset, g.numPoints);
+ // point 0 is 0.0
+ *g.pointY = 0.0;
+
+ CHECK_VAR(ncFile, "cellsOnVertex");
+ g.origConnections = (int *) malloc(g.numCells * g.pointsPerCell *
+ sizeof(int));
+ CHECK_MALLOC(g.origConnections);
+ NcVar *connectionsVar = ncFile.get_var("cellsOnVertex");
+ connectionsVar->get(g.origConnections, g.numCells, g.pointsPerCell);
+
+ // create my own list to include modified origConnections (due to
+ // eliminating wraparound in the lat/lon projection) plus additional
+ // cells added when mirroring cells that had previously wrapped around
+
+ g.modConnections = (int *) malloc(g.modNumCells * g.pointsPerCell
+ * sizeof(int));
+ CHECK_MALLOC(g.modConnections);
+
+
+ // allocate an array to map the extra points and cells to the original
+ // so that when obtaining data, we know where to get it
+ g.pointMap = (int*)malloc((int)floor(g.numPoints*BLOATFACTOR)
+ * sizeof(int));
+ CHECK_MALLOC(g.pointMap);
+ g.cellMap = (int*)malloc((int)floor(g.numCells*BLOATFACTOR)
+ * sizeof(int));
+ CHECK_MALLOC(g.cellMap);
+
+ g.currentExtraPoint = g.numPoints + g.pointOffset;
+ g.currentExtraCell = g.numCells + g.cellOffset;
+
+ if (p.includeTopography) {
+ CHECK_VAR(ncFile, "maxLevelCell");
+ g.maxLevelPoint = (int*)malloc((g.numPoints + g.numPoints) * sizeof(int));
+ CHECK_MALLOC(g.maxLevelPoint);
+ NcVar *maxLevelPointVar = ncFile.get_var("maxLevelCell");
+ maxLevelPointVar->get(g.maxLevelPoint + g.pointOffset, g.numPoints);
+ }
+}
+
+void ShiftLonData(params &p, geometry &g)
+{
+ // if atmospheric data, or zero centered, set center to 180 instead of 0
+ if (p.isAtmosphere || p.isZeroCentered) {
+ for (int j = g.pointOffset; j < g.numPoints + g.pointOffset; j++) {
+ // need to shift over the point so center is at PI
+ if (g.pointX[j] < 0) {
+ g.pointX[j] += 2*PI;
+ }
+ }
+ }
+
+ if (p.centerLon != 180) {
+ for (int j = g.pointOffset; j < g.numPoints + g.pointOffset; j++) {
+ // need to shift over the point if centerLon dictates
+ if (p.centerRad < PI) {
+ if (g.pointX[j] > (p.centerRad + PI)) {
+ g.pointX[j] = -((2*PI) - g.pointX[j]);
+ }
+ } else if (p.centerRad > PI) {
+ if (g.pointX[j] < (p.centerRad - PI)) {
+ g.pointX[j] += 2*PI;
+ }
+ }
+ }
+ }
+}
+
+
+int AddMirrorPoint(int index, double dividerX, geometry &g)
+{
+ double X = g.pointX[index];
+ double Y = g.pointY[index];
+
+ // add on east
+ if (X < dividerX) {
+ X += 2*PI;
+ } else {
+ // add on west
+ X -= 2*PI;
+ }
+
+ g.pointX[g.currentExtraPoint] = X;
+ g.pointY[g.currentExtraPoint] = Y;
+
+ int mirrorPoint = g.currentExtraPoint;
+
+ // record mapping
+ *(g.pointMap + (g.currentExtraPoint - g.numPoints - g.pointOffset)) = index;
+ g.currentExtraPoint++;
+
+ return mirrorPoint;
+}
+
+// Check for out-of-range values and do bugfix
+void FixPoints(params &p, geometry &g)
+{
+
+ for (int j = g.cellOffset; j < g.numCells + g.cellOffset; j++ ) {
+ int *conns = g.origConnections + (j * g.pointsPerCell);
+
+ // go through and make sure none of the referenced points are
+ // out of range
+ // if so, set all to point 0
+ for (int k = 0; k < g.pointsPerCell; k++) {
+ if ((conns[k] <= 0) || (conns[k] > g.numPoints)) {
+ for (int m = 0; m < g.pointsPerCell; m++) {
+ conns[m] = 0;
+ }
+ break;
+ }
+ }
+
+ int lastk;
+
+ if (p.doBugFix) {
+ //BUG FIX for problem where cells are stretching to a faraway point
+ lastk = g.pointsPerCell-1;
+ const double thresh = .06981317007977; // 4 degrees
+ for (int k = 0; k < g.pointsPerCell; k++) {
+ double xdiff = abs(g.pointX[conns[k]]
+ - g.pointX[conns[lastk]]);
+ double ydiff = abs(g.pointY[conns[k]]
+ - g.pointY[conns[lastk]]);
+ // Don't look at cells at map border
+ if (ydiff > thresh) {
+ for (int m = 0; m < g.pointsPerCell; m++) {
+ conns[m] = 0;
+ }
+ break;
+ }
+ }
+ }
+ }
+}
+
+
+// Eliminate wraparound at east/west edges of lat/lon projection
+void EliminateXWrap (params &p, geometry &g)
+{
+
+ // For each cell, examine vertices
+ // Add new points and cells where needed to account for wraparound.
+
+ int mirrorcell;
+ int apoint;
+ int mirrorpoint;
+
+
+ for (int j = g.cellOffset; j < g.numCells + g.cellOffset; j++ ) {
+ int *conns = g.origConnections + (j * g.pointsPerCell);
+ int *modConns = g.modConnections + (j * g.pointsPerCell);
+ int lastk;
+
+ // Determine if we are wrapping in X direction
+ lastk = g.pointsPerCell-1;
+ bool xWrap = false;
+ for (int k = 0; k < g.pointsPerCell; k++) {
+ if (abs(g.pointX[conns[k]]
+ - g.pointX[conns[lastk]]) > 5.5) xWrap = true;
+ lastk = k;
+ }
+
+ // If we wrapped in X direction, modify cell and add mirror cell
+ if (xWrap) {
+
+ // first point is anchor it doesn't move
+ double anchorX = g.pointX[conns[0]];
+ double anchorY = g.pointY[conns[0]];
+ modConns[0] = conns[0];
+
+ // modify existing cell, so it doesn't wrap
+ // move points to one side
+ for (int k = 1; k < g.pointsPerCell; k++) {
+ int neigh = conns[k];
+
+ // add a new point, figure out east or west
+ if (abs(g.pointX[neigh] - anchorX) > 5.5) {
+ modConns[k] = AddMirrorPoint(neigh, anchorX, g);
+ } else {
+ // use existing kth point
+ modConns[k] = neigh;
+ }
+ }
+
+ // move addedConns to modConnections extra cells area
+ int* addedConns = g.modConnections
+ + (g.currentExtraCell * g.pointsPerCell);
+
+ // add a mirroring cell to other side
+
+ // add mirrored anchor first
+ addedConns[0] = AddMirrorPoint(conns[0], p.centerRad, g);
+ anchorX = g.pointX[addedConns[0]];
+
+ // add mirror cell points if needed
+ for (int k = 1; k < g.pointsPerCell; k++) {
+ int neigh = conns[k];
+
+ // add a new point for neighbor, figure out east or west
+ if (abs(g.pointX[neigh] - anchorX) > 5.5) {
+ addedConns[k] = AddMirrorPoint(neigh, anchorX, g);
+ } else {
+ // use existing kth point
+ addedConns[k] = neigh;
+ }
+ }
+ *(g.cellMap + (g.currentExtraCell - g.numCells - g.cellOffset)) = j;
+ g.currentExtraCell++;
+ } else {
+
+ // just add cell "as is" to modConnections
+ for (int k=0; k< g.pointsPerCell; k++) {
+ modConns[k] = conns[k];
+ }
+ }
+ if (g.currentExtraCell > g.modNumCells) {
+ cerr << "Exceeded storage for extra cells!" << endl;
+ exit(1);
+ }
+ if (g.currentExtraPoint > g.modNumPoints) {
+ cerr << "Exceeded storage for extra points!" << endl;
+ exit(1);
+ }
+ }
+
+ if (p.singleLayer) {
+ g.maxCells = g.currentExtraCell;
+ g.maxPoints = g.currentExtraPoint;
+ } else {
+ g.maxCells = g.currentExtraCell*g.maxNVertLevels;
+ g.maxPoints = g.currentExtraPoint*(g.maxNVertLevels+1);
+ }
+}
+
+void WritePoints(ofstream &geoFile, EncodeBase64 &b64, params &p, geometry &g)
+{
+ // write points (the points are the primal-grid cell centers)
+
+ // write the point at each vertical level, plus one level for last layer
+
+ geoFile << " <Points>" << endl;
+
+ string formatString;
+
+ if (p.writeBinary) {
+ formatString = "binary";
+ } else {
+ formatString = "ascii";
+ }
+
+ // see if we can change to double
+
+ //geoFile << " <DataArray type=\"Float32\" Name=\"Points\" NumberOfComponents=\"3\" format=\"" << formatString << "\">" << endl;
+ geoFile << " <DataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"3\" format=\"" << formatString << "\">" << endl;
+
+ geoFile.precision(16);
+
+ geoFile << " ";
+
+ if (p.writeBinary) {
+ b64.StartWriting();
+ //int binarySize = g.maxPoints*sizeof(float)*3;
+ int binarySize = g.maxPoints*sizeof(double)*3;
+ binarySize = PadSize(binarySize);
+ b64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ //cout << "Writing points at each level" << endl;
+
+ for (int j = 0; j < g.currentExtraPoint; j++ )
+ {
+
+ //float x, y, z;
+ double x, y, z;
+
+ if (p.projectLatLon) {
+ x = g.pointX[j] * 180.0 / PI;
+ y = g.pointY[j] * 180.0 / PI;
+ z = 0.0;
+ } else {
+ x = g.pointX[j];
+ y = g.pointY[j];
+ z = g.pointZ[j];
+ }
+
+ if (p.singleLayer) {
+ if (p.writeBinary) {
+ //b64.Write((const unsigned char*)&x, sizeof(float));
+ b64.Write((const unsigned char*)&x, sizeof(double));
+ //b64.Write((const unsigned char*)&y, sizeof(float));
+ b64.Write((const unsigned char*)&y, sizeof(double));
+ //b64.Write((const unsigned char*)&z, sizeof(float));
+ b64.Write((const unsigned char*)&z, sizeof(double));
+ } else {
+ geoFile << " " << x << " " << y << " " << z << endl;
+ }
+ } else {
+ double rho, rholevel, theta, phi;
+ int retval = -1;
+
+ if (!p.projectLatLon) {
+ if ((x != 0.0) || (y != 0.0) || (z != 0.0)) {
+ retval = CartesianToSpherical(x, y, z, &rho, &phi, &theta);
+ }
+ }
+
+ for (int levelNum = 0; levelNum < g.maxNVertLevels+1; levelNum++) {
+ if (p.projectLatLon){
+ z = -(((double)levelNum)*p.layerThickness);
+ } else {
+ if (!retval && ((x != 0.0) || (y != 0.0) || (z != 0.0))) {
+ rholevel = rho - (p.layerThickness * levelNum);
+ SphericalToCartesian(rholevel, phi, theta, &x, &y, &z);
+ }
+ }
+
+ if (p.writeBinary) {
+ //b64.Write((const unsigned char*)&x, sizeof(float));
+ b64.Write((const unsigned char*)&x, sizeof(double));
+ //b64.Write((const unsigned char*)&y, sizeof(float));
+ b64.Write((const unsigned char*)&y, sizeof(double));
+ //b64.Write((const unsigned char*)&z, sizeof(float));
+ b64.Write((const unsigned char*)&z, sizeof(double));
+ } else {
+ geoFile << " " << x << " " << y << " " << z << endl;
+ }
+ }
+ }
+ }
+
+ if (p.writeBinary) {
+ b64.EndWriting();
+ geoFile << endl;
+ }
+
+ geoFile << " </DataArray>" << endl;
+ geoFile << " </Points>" << endl;
+}
+
+
+void WriteConnectivity(ofstream &geoFile, EncodeBase64 &b64, params &p,
+ geometry &g)
+{
+
+ // Write cell connectivity
+
+ string formatString;
+ if (p.writeBinary) {
+ formatString = "binary";
+ } else {
+ formatString = "ascii";
+ }
+
+ geoFile << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"" << formatString << "\">" << endl;
+
+ // for each cell, write the points for each
+
+
+ if (p.writeBinary) {
+ geoFile << " ";
+ b64.StartWriting();
+ int binarySize;
+ if (p.singleLayer) {
+ binarySize = g.maxCells*g.pointsPerCell*sizeof(int);
+ } else {
+ binarySize = g.maxCells*g.pointsPerCell*2*sizeof(int);
+ }
+ binarySize = PadSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ b64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ // for each cell, write number of points for each
+ // and then list the points by number
+
+ //cout << "Writing connectivity" << endl;
+
+ unsigned int val;
+
+ for (int j = 0; j < g.currentExtraCell ; j++) {
+
+ int* conns;
+ if (p.projectLatLon) {
+ conns = g.modConnections + (j * g.pointsPerCell);
+ } else {
+ conns = g.origConnections + (j * g.pointsPerCell);
+ }
+
+ int minLevel= 0;
+
+ if (p.includeTopography) {
+ int* connections;
+
+ //check if it is a mirror cell, if so, get original
+ if (j >= g.numCells + g.cellOffset) {
+ //cout << "setting origCell" << endl;
+ int origCellNum = *(g.cellMap + (j - g.numCells - g.cellOffset));
+ //cout << "mirror cell: " <<j<< " origCell: "<< origCell << endl;
+ conns = g.origConnections + (origCellNum*g.pointsPerCell);
+ } else connections = g.origConnections + (j * g.pointsPerCell);
+
+ //cout << "calc minLevel" << endl;
+ minLevel = g.maxLevelPoint[conns[0]];
+ //cout << "initial minLevel:" << minLevel << endl;
+ // Take the min of the maxLevelPoint of each point
+ for (int k = 1; k < g.pointsPerCell; k++) {
+ minLevel = min(minLevel, g.maxLevelPoint[connections[k]]);
+ }
+ //cout << endl;
+ }
+
+ if (p.singleLayer) {
+ if (!p.writeBinary) geoFile << " ";
+
+ // If that min is greater than or equal to this output level,
+ // include the cell, otherwise set all points to zero.
+
+ if (p.includeTopography && ((minLevel-1) < p.outputVertLevel)) {
+ //cerr << "Setting all points to zero" << endl;
+ val = 0;
+ for (int k = 0; k < g.pointsPerCell; k++) {
+ WriteInt(b64, geoFile, val, p.writeBinary);
+ }
+ } else {
+ for (int k = 0; k < g.pointsPerCell; k++) {
+ WriteInt(b64, geoFile, conns[k], p.writeBinary);
+ }
+ }
+
+ if (!p.writeBinary) geoFile << endl;
+ } else { // multilayer
+
+ // for each level, write the cell
+ for (int levelNum = 0; levelNum < g.maxNVertLevels; levelNum++) {
+ if(!p.writeBinary) geoFile << " ";
+
+ if (p.includeTopography && ((minLevel-1) < levelNum)) {
+ // setting all points to zero
+ val = 0;
+ for (int m = 0; m < g.pointsPerCell*2; m++) {
+ WriteInt(b64, geoFile, val, p.writeBinary);
+ }
+
+ } else {
+
+ for (int k = 0; k < g.pointsPerCell; k++)
+ {
+ val = (conns[k]*(g.maxNVertLevels+1)) + levelNum;
+ WriteInt(b64, geoFile, val, p.writeBinary);
+ }
+
+ for (int k = 0; k < g.pointsPerCell; k++)
+ {
+ val = (conns[k]*(g.maxNVertLevels+1)) + levelNum + 1;
+ WriteInt(b64, geoFile, val, p.writeBinary);
+ }
+ }
+ if (!p.writeBinary) geoFile << endl;
+ }
+ }
+ }
+
+ if (p.writeBinary) {
+ b64.EndWriting();
+ geoFile << endl;
+ }
+
+ geoFile << " </DataArray>" << endl;
+
+}
+
+
+void WriteOffsets(ofstream &geoFile, EncodeBase64 &b64, params &p, geometry &g)
+{
+ string formatString;
+ if (p.writeBinary) {
+ formatString = "binary";
+ } else {
+ formatString = "ascii";
+ }
+
+ // Specify the offset into connectivity array for the end of each cell
+ geoFile << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"" << formatString << "\">";
+
+ if (p.writeBinary) {
+ geoFile << endl << " ";
+ b64.StartWriting();
+ int binarySize = g.maxCells*sizeof(int);
+ binarySize = PadSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ b64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ int counter = 0;
+ unsigned int offset = 0;
+
+ for (int j = 0; j < g.maxCells ; j++) {
+ if ((!p.writeBinary) && (counter%6 == 0)) geoFile << endl << " ";
+ if (p.singleLayer) {
+ offset += g.pointsPerCell;
+ } else {
+ offset += 2*g.pointsPerCell;
+ }
+ WriteInt(b64, geoFile, offset, p.writeBinary);
+ counter++;
+ }
+
+ if (p.writeBinary) {
+ b64.EndWriting();
+ }
+
+ geoFile << endl;
+ geoFile << " </DataArray>" << endl;
+}
+
+
+void WriteCellTypes(ofstream &geoFile, EncodeBase64 &b64,params &p, geometry &g)
+{
+ // write cell types
+ unsigned char cellType;
+ switch (g.pointsPerCell) {
+ case 3:
+ if (p.singleLayer) {
+ cellType = VTK_TRIANGLE;
+ } else {
+ cellType = VTK_WEDGE;
+ }
+ break;
+ case 4:
+ if (p.singleLayer) {
+ cellType = VTK_QUAD;
+ } else {
+ cellType = VTK_HEXAHEDRON;
+ }
+ break;
+ default:
+ break;
+ }
+
+ string formatString;
+ if (p.writeBinary) {
+ formatString = "binary";
+ } else {
+ formatString = "ascii";
+ }
+
+ geoFile << " <DataArray type=\"UInt8\" Name=\"types\" format=\"" <<
+ formatString << "\">";
+
+ if (p.writeBinary) {
+ geoFile << endl << " ";
+ b64.StartWriting();
+ int binarySize = g.maxCells*sizeof(char);
+ binarySize = PadSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ b64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ for (int j = 0; j < g.maxCells ; j++) {
+ if ((!p.writeBinary) && (j%6 == 0)) {
+ geoFile << endl << " ";
+ }
+ WriteChar(b64, geoFile, cellType, p.writeBinary);
+ }
+
+ if (p.writeBinary) b64.EndWriting();
+
+ geoFile << endl;
+ geoFile << " </DataArray>" << endl;
+}
+
+
+void WriteGeometry(ostringstream &geoFileName, params &p, geometry &g)
+{
+ geoFileName << "geo_" << p.outputFileName;
+
+ ofstream geoFile(geoFileName.str().c_str(), ios::out);
+ if (!geoFile)
+ {
+ cerr << "output file " << geoFileName << " could not be opened" <<endl;
+ exit(1);
+ }
+
+ EncodeBase64 b64(&geoFile);
+
+ WritePoints(geoFile, b64, p, g);
+
+ geoFile << " <Cells>" << endl;
+
+ WriteConnectivity(geoFile, b64, p, g);
+ WriteOffsets(geoFile, b64, p, g);
+ WriteCellTypes(geoFile, b64, p, g);
+
+ geoFile << " </Cells>" << endl;
+ geoFile << " </Piece>" << endl;
+ geoFile << " </UnstructuredGrid>" << endl;
+ geoFile << "</VTKFile>" << endl;
+ geoFile.close();
+}
+
+
+void DeallocBasicGeometry(params &p, geometry &g)
+{
+ free (g.pointX);
+ free (g.pointY);
+ free (g.origConnections);
+
+ if (p.projectLatLon) {
+ free(g.modConnections);
+ } else {
+ free (g.pointZ);
+ }
+
+ if (p.includeTopography) {
+ free(g.maxLevelPoint);
+ }
+}
+
+
+void WritePointVarData(ofstream &vtkFile, EncodeBase64 &vtkb64, int t,
+ params &p, geometry &g, vars &v)
+{
+ // If by point, write out point data
+ vtkFile << " <PointData>" << endl;
+
+ double* pointVarData;
+ pointVarData = (double*)malloc((sizeof(double)) * g.currentExtraPoint
+ * g.maxNVertLevels);
+ CHECK_MALLOC(pointVarData);
+
+ int vertLevelIndex = 0;
+ if (p.singleLayer) vertLevelIndex = p.outputVertLevel;
+
+ string formatString;
+
+ if (p.writeBinary) {
+ formatString = "binary";
+ } else {
+ formatString = "ascii";
+ }
+
+ for (int i = 0; i < v.numPointVars; i++) {
+
+ // Read variable number i data for that timestep
+
+ v.pointVars[i]->set_cur(t, 0, vertLevelIndex);
+ if (p.singleLayer) {
+ v.pointVars[i]->get(pointVarData+g.pointOffset, 1, g.numPoints, 1);
+ } else {
+ v.pointVars[i]->get(pointVarData +
+ (g.maxNVertLevels * g.pointOffset), 1,
+ g.numPoints, g.maxNVertLevels);
+ }
+
+ // Write variable number i data for that timestep
+
+ //vtkFile << " <DataArray type=\"Float32\" Name=\"" <<
+ vtkFile << " <DataArray type=\"Float64\" Name=\"" <<
+ v.pointVars[i]->name();
+ vtkFile << "\" format=\"" << formatString << "\"> " << endl;
+
+
+ vtkFile << " ";
+
+ double *var_target;
+ //float validData;
+ double validData;
+
+ if (p.writeBinary) {
+ vtkb64.StartWriting();
+ //int binarySize = g.maxPoints*sizeof(float);
+ int binarySize = g.maxPoints*sizeof(double);
+ binarySize = PadSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ vtkb64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ // if (g.pointOffset == 1)
+ //write dummy
+ // 0th pointVarData value is a filler, zero, which is not necessarily in
+ // the range of the other data, so we look at the 1st value instead for
+ // our dummy
+ if (p.singleLayer) {
+ var_target = pointVarData + 1;
+ //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
+ WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
+ } else {
+ var_target = pointVarData + g.maxNVertLevels;
+ for (int levelNum = 0; levelNum < g.maxNVertLevels; levelNum++) {
+ //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
+ WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
+ var_target++;
+ }
+
+ // write highest level dummy point (duplicate of last level)
+ var_target--;
+ //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
+ WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
+
+ // Now set the var_target to the 1st data value
+ var_target = pointVarData + g.maxNVertLevels;
+ }
+
+ for (int j = g.pointOffset; j < g.numPoints + g.pointOffset; j++) {
+
+ if (p.singleLayer) {
+ //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
+ WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
+ var_target++;
+ } else {
+
+ // write data for one point -- lowest level to highest
+ for (int levelNum = 0; levelNum < g.maxNVertLevels; levelNum++) {
+ //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
+ WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
+ var_target++;
+ }
+
+ // for last layer of points, repeat last level's values
+ // Need Mark's input on this one
+ var_target--;
+ //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
+ WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
+ var_target++;
+ }
+ }
+
+ // put out data for extra points
+ for (int j = g.pointOffset + g.numPoints; j < g.currentExtraPoint; j++) {
+ // use map to find out what point data we are using
+ var_target = pointVarData + ((*(g.pointMap + j - g.numPoints - g.pointOffset))*g.maxNVertLevels);
+
+ if (p.singleLayer) {
+ //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
+ WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
+ //validData = ConvertDouble2ValidFloat (*var_target);
+ validData = ConvertDouble2ValidDouble (*var_target);
+ var_target++;
+ } else {
+ // write data for one point -- lowest level to highest
+ for (int levelNum = 0; levelNum < g.maxNVertLevels; levelNum++) {
+ //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
+ WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
+ var_target++;
+ }
+
+ // for last layer of points, repeat last level's values
+ // Need Mark's input on this one
+ var_target--;
+ //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
+ WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
+ }
+ }
+ if (p.writeBinary) vtkb64.EndWriting();
+ vtkFile << endl << " </DataArray>" << endl;
+
+ }
+
+ free (pointVarData);
+
+ vtkFile << " </PointData>" << endl;
+
+}
+
+
+void WriteCellVarData(ofstream &vtkFile, EncodeBase64 &vtkb64, int t,
+ params &p, geometry &g, vars &v)
+{
+ // if by cell, then write out cell data
+ double* cellVarData;
+ // use maxCell instead?
+ cellVarData = (double*)malloc((sizeof(double)) * g.currentExtraCell
+ * g.maxNVertLevels);
+ CHECK_MALLOC(cellVarData);
+
+ vtkFile << " <CellData>" << endl;
+
+ int vertLevelIndex = 0;
+ if (p.singleLayer) vertLevelIndex = p.outputVertLevel;
+
+ string formatString;
+
+ if (p.writeBinary) {
+ formatString = "binary";
+ } else {
+ formatString = "ascii";
+ }
+
+ for (int i = 0; i < v.numCellVars; i++) {
+
+ // Write variable number v data for that timestep
+ //vtkFile << " <DataArray type=\"Float32\" Name=\"" <<
+ vtkFile << " <DataArray type=\"Float64\" Name=\"" <<
+ v.cellVars[i]->name() << "\" format=\"" << formatString << "\"> " ;
+
+ // Read variable number v data for that timestep and level
+ v.cellVars[i]->set_cur(t, 0, vertLevelIndex);
+ if (p.singleLayer) {
+ v.cellVars[i]->get(cellVarData + g.cellOffset, 1, g.numCells,1);
+ } else {
+ v.cellVars[i]->get(cellVarData + (g.cellOffset * g.maxNVertLevels),
+ 1, g.numCells, g.maxNVertLevels);
+ }
+
+ if (p.writeBinary) {
+ vtkFile << endl << " ";
+ vtkb64.StartWriting();
+ //int binarySize = g.maxCells*sizeof(float);
+ int binarySize = g.maxCells*sizeof(double);
+ binarySize = PadSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ vtkb64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ double *var_target = cellVarData;
+ int counter = 0;
+
+ for (int j = g.cellOffset; j < g.numCells + g.cellOffset; j++) {
+
+ if (p.singleLayer) {
+ if ((!p.writeBinary) && (j%6 == 0)) vtkFile << endl
+ << " ";
+ //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
+ WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
+ var_target++;
+ } else {
+ for (int levelNum = 0; levelNum < g.maxNVertLevels; levelNum++)
+ {
+ if ((!p.writeBinary) && (counter%6 == 0)) vtkFile << endl
+ << " ";
+ //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
+ WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
+ var_target++;
+ counter++;
+ }
+ }
+ }
+
+ // put out for extra cells
+ counter = 0;
+ for (int j = g.numCells + g.cellOffset; j < g.currentExtraCell; j++) {
+
+ var_target = cellVarData
+ + (*(g.cellMap + j - g.numCells - g.cellOffset))
+ *g.maxNVertLevels;
+ if (p.singleLayer) {
+ if ((!p.writeBinary) && (j%6 == 0))
+ vtkFile << endl << " ";
+ //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
+ WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
+ } else {
+ for (int levelNum = 0; levelNum < g.maxNVertLevels; levelNum++)
+ {
+ if ((!p.writeBinary) && (counter%6 == 0))
+ vtkFile << endl << " ";
+ //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
+ WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
+ var_target++;
+ counter++;
+ }
+ }
+ }
+
+ if (p.writeBinary) {
+ vtkb64.EndWriting();
+ }
+ vtkFile << endl << " </DataArray>" << endl;
+
+ }
+
+ vtkFile << " </CellData>" << endl;
+
+ free(cellVarData);
+}
+
+
+void WriteVarData(ostringstream &vtkFileName, int t, params &p, geometry &g,
+ vars &v)
+{
+ vtkFileName << t << p.outputFileName;
+ ofstream vtkFile(vtkFileName.str().c_str(), ios::out);
+ if (!vtkFile)
+ {
+ cerr << "vtk output file could not be opened" <<endl;
+ exit(1);
+ }
+
+ vtkFile.precision(16);
+
+ vtkFile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << endl;
+ vtkFile << " <UnstructuredGrid>" << endl;
+
+ vtkFile << " <Piece NumberOfPoints=\"" << g.maxPoints << "\" NumberOfCells=\"" << g.maxCells << "\">" << endl;
+
+ EncodeBase64 vtkb64(&vtkFile);
+
+ WritePointVarData(vtkFile, vtkb64, t, p, g, v);
+ WriteCellVarData(vtkFile, vtkb64, t, p, g, v);
+
+ vtkFile.close();
+}
+
+
+int main(int argc, char* argv[])
+{
+ if (argc < 3) {
+ PrintUsage();
+ exit(1);
+ }
+
+ struct params p;
+ GetArgs(argc, argv, p);
+
+ NcFile ncFile(p.inputFileName);
+ if (!ncFile.is_valid())
+ {
+ cerr << "Couldn't open file: " << p.inputFileName << endl;
+ exit(1);
+ }
+
+ // figure out geometry of a dual-grid lat/lon projection
+ geometry g;
+ GetNcDims(ncFile, g);
+ CheckArgs(p, g);
+ if (p.projectLatLon) {
+ AllocLatLonGeometry(ncFile, p, g);
+ ShiftLonData(p, g);
+ FixPoints(p, g);
+ EliminateXWrap(p, g);
+ } else {
+ AllocSphereGeometry(ncFile, p, g);
+ }
+
+ // write a file with the geometry.
+ ostringstream geoFileName;
+ WriteGeometry(geoFileName, p, g);
+ DeallocBasicGeometry(p, g);
+
+ // figure out what variables to visualize
+ vars v;
+ GetNcVars(ncFile, v, "nVertices", "nCells");
+
+ // For each timestep, write data for each variable
+ for (int t = p.startTimeStep; t <= p.finishTimeStep; t++) {
+
+ ostringstream vtkFileName;
+
+ WriteVarData(vtkFileName, t, p, g, v);
+
+ // append geoFile to vtkFile
+ string geoFileNameString = geoFileName.str();
+ string vtkFileNameString = vtkFileName.str();
+ int appRetVal = AppendFile(&geoFileNameString, &vtkFileNameString);
+ }
+
+ return(0);
+}
+
Property changes on: branches/ocean_projects/graphics/paraview/translator/mpasnc2vtu.cpp
___________________________________________________________________
Name: svn:mime-type
+ text/x-c++src
Name: svn:keywords
+ Author Date ID Revision
Name: svn:eol-style
+ native
</font>
</pre>