<p><b>ringler@lanl.gov</b> 2009-12-17 08:32:09 -0700 (Thu, 17 Dec 2009)</p><p>a tool to map output.nc files to paraview viz files<br>
</p><hr noshade><pre><font color="gray">Added: branches/ocean/visualization/translator/Makefile
===================================================================
--- branches/ocean/visualization/translator/Makefile (rev 0)
+++ branches/ocean/visualization/translator/Makefile 2009-12-17 15:32:09 UTC (rev 86)
@@ -0,0 +1,12 @@
+# NETCDF should point to the directory holding the netcdf lib and bin subdirs
+
+CC = gcc
+INCLUDES = -I$(NETCDF)/include
+LIBS = -L$(NETCDF)/lib -lnetcdf_c++ -lnetcdf -L/usr/lib -lstdc++
+
+
+transdual: nc2vtk_dual.cpp
+ $(CC) $(INCLUDES) -o $@ $< $(LIBS)
+
+clean:
+ rm transdual
Added: branches/ocean/visualization/translator/README
===================================================================
--- branches/ocean/visualization/translator/README (rev 0)
+++ branches/ocean/visualization/translator/README 2009-12-17 15:32:09 UTC (rev 86)
@@ -0,0 +1,16 @@
+Here are the steps
+
+1. build executable "transdual" by editing Makefile then "make"
+(you will have to insure that the compiler used in the make file was also used to build your netCDF)
+
+2. usage: transdual infile.nc output.vtk
+
+3. download and install Paraview (binaries available online at http://www.paraview.org/paraview/resources/software.html )
+
+4. boot Paraview, open one of the vtk files (open "0" as it is the first in the time series)
+
+5. click "apply" on LHS
+
+6. viz at will!
+
+12-17-2009: We will be adding support for multiple vertical levels soon
Added: branches/ocean/visualization/translator/nc2vtk_dual.cpp
===================================================================
--- branches/ocean/visualization/translator/nc2vtk_dual.cpp (rev 0)
+++ branches/ocean/visualization/translator/nc2vtk_dual.cpp 2009-12-17 15:32:09 UTC (rev 86)
@@ -0,0 +1,328 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+// This program translates a newpop netCDF data file to dual-grid legacy, ascii VTK format
+// The variables that have time dim are automatically written.
+//
+// Assume all variables are of interest.
+// Assume variable data type is double.
+// Assume variable dims are (Time, nCells|nVertices, nVertLevels)
+// Assume no more than 10 vars each for cell and point data
+// Does not deal with tracers.
+// Does not deal with edge data.
+// Does not vis vertical levels
+//
+// Christine Ahrens
+// 10/16/2009
+//
+////////////////////////////////////////////////////////////////////////////////
+
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include "stdlib.h"
+#include "vtkCellType.h"
+#include "netcdfcpp.h"
+#include <string>
+#include <cmath>
+
+using namespace std;
+
+#define CHECK_MALLOC(ptr) \
+ if (ptr == NULL) { \
+ cerr << "malloc failed!</font>
<font color="blue">"; \
+ exit(1); \
+ }
+#define MAX_VARS 10
+
+int main(int argc, char* argv[])
+{
+ if (argc != 3) {
+ cerr << "Usage: transdual infile.nc outfile.vtk" << endl;
+ cerr << "Variables with time dim are written out." << endl;
+ cerr << "A series of files 0outfile.vtk, 1outfile.vtk, etc. will" << endl << "be created, one file for each time step." << endl;
+ exit(1);
+ }
+
+ NcFile ncFile(argv[1]);
+
+ if (!ncFile.is_valid())
+ {
+ cerr << "Couldn't open file: " << argv[1] << endl;
+ exit(1);
+ }
+
+ NcDim* nCells = ncFile.get_dim("nCells");
+ NcDim* nVertices = ncFile.get_dim("nVertices");
+ NcDim* THREE = ncFile.get_dim("THREE");
+ NcDim* Time = ncFile.get_dim("Time");
+ NcDim* nVertLevels = ncFile.get_dim("nVertLevels");
+
+ // figure out what variables to visualize -- array of NcVars?
+ NcVar* cellVars[MAX_VARS];
+ NcVar* pointVars[MAX_VARS];
+ int cellVarIndex = -1;
+ int pointVarIndex = -1;
+ int numCells = nVertices->size();
+ int numPoints = nCells->size()+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)
+
+ int numDims = aVar->num_dims();
+ //cout << "Num Dims of var: " << aVar->name() << " is " << numDims << endl;
+ if (numDims != 3) {
+ continue; // try the next var
+ } else {
+ // TODO, check if it is a double
+ // assume a double for now
+
+ // check for Time dim
+ NcToken aName = aVar->get_dim(0)->name();
+
+ // TODO check if third dim is nVertLevels, too
+
+ // Add to cell or point var array
+ if (!strcmp(aName, "Time")) {
+
+ aName = aVar->get_dim(1)->name();
+ if (!strcmp(aName, "nVertices")) {
+ cellVarIndex++;
+ // TODO check index no higher than MAX_VARS-1
+ cellVars[cellVarIndex] = aVar;
+ //cout << "Adding var " << aVar->name() << " to cellVars" << endl;
+ } else if (!strcmp(aName, "nCells")) {
+ pointVarIndex++;
+ // TODO check index no higher than MAX_VARS-1
+ pointVars[pointVarIndex] = aVar;
+ //cout << "Adding var " << aVar->name() << " to pointVars" << endl;
+ }
+ }
+ }
+ }
+
+ // TODO
+ // prompt the user to find out which fields are of interest?
+ // for now, assume all are of interest
+
+ // get points (centers of primal-mesh cells)
+
+ // TO DO check malloc return vals.
+
+ double *xCellData = (double*)malloc(nCells->size()
+ * sizeof(double));
+ CHECK_MALLOC(xCellData);
+ NcVar *xCellVar = ncFile.get_var("xCell");
+ xCellVar->get(xCellData, nCells->size());
+
+ double *yCellData = (double*)malloc(nCells->size()
+ * sizeof(double));
+ CHECK_MALLOC(yCellData);
+ NcVar *yCellVar = ncFile.get_var("yCell");
+ yCellVar->get(yCellData, nCells->size());
+
+ double *zCellData = (double*)malloc(nCells->size()
+ * sizeof(double));
+ //cout << "ptr for zCellData" << zCellData << endl;
+ CHECK_MALLOC(zCellData);
+ NcVar *zCellVar = ncFile.get_var("zCell");
+ zCellVar->get(zCellData, nCells->size());
+
+ // get dual-mesh cells
+
+ int *cellsOnVertex = (int *) malloc((nVertices->size()) * THREE->size() *
+ sizeof(int));
+ //cout << "ptr for cellsOnVertex" << cellsOnVertex << endl;
+ CHECK_MALLOC(cellsOnVertex);
+ NcVar *cellsOnVertexVar = ncFile.get_var("cellsOnVertex");
+ //cout << "getting cellsOnVertexVar</font>
<font color="gray">";
+ cellsOnVertexVar->get(cellsOnVertex, nVertices->size(), THREE->size());
+
+ // decls for variables
+ double* cellVarData[MAX_VARS];
+ double* pointVarData[MAX_VARS];
+
+ // for each variable, allocate space for variables
+
+ //cout << "cellVarIndex: " << cellVarIndex << endl;
+ for (int v = 0; v <= cellVarIndex; v++) {
+ //cout << "Allocating for cellvar: " << cellVars[v]->name() << endl;
+ cellVarData[v] = (double*)malloc((sizeof(double))
+ * nVertices->size() * nVertLevels->size());
+ //cout << "ptr for " << cellVars[v]->name() << " " << cellVarData[v] << endl;
+ //cout << "nVertices: " << nVertices->size() << endl;
+ //cout << "Time: " << Time->size() << endl;
+ //cout << "nVertLevels: " << nVertLevels->size() << endl;
+ CHECK_MALLOC(cellVarData[v]);
+ }
+
+ //cout << "pointVarIndex: " << pointVarIndex << endl;
+ for (int v = 0; v <= pointVarIndex; v++) {
+ //cout << "Allocating for point var: " << pointVars[v]->name() << endl;
+ pointVarData[v] = (double*)malloc((sizeof(double))
+ * nCells->size() * nVertLevels->size());
+ CHECK_MALLOC(pointVarData[v]);
+ }
+
+ // for each time step, write a file with the time prepended to filename
+ for (int i = 0; i < Time->size(); i++)
+ {
+
+ ostringstream vtkFileName;
+
+ vtkFileName << i << argv[2];
+
+ ofstream vtkFile(vtkFileName.str().c_str(), ios::out);
+
+ if (!vtkFile)
+ {
+ cerr << "vtk output file could not be opened" <<endl;
+ exit(1);
+ }
+
+
+ // Read data for that timestep for each variable
+
+ for (int v = 0; v <= cellVarIndex; v++) {
+ cellVars[v]->set_cur(i, 0, 0);
+ cellVars[v]->get(cellVarData[v], 1,
+ nVertices->size(), nVertLevels->size());
+ }
+
+ for (int v = 0; v <= pointVarIndex; v++) {
+ pointVars[v]->set_cur(i, 0, 0);
+ pointVars[v]->get(pointVarData[v], 1,
+ nCells->size(), nVertLevels->size());
+ }
+
+ // write header
+
+ vtkFile << "# vtk DataFile Version 2.0" << endl;
+ vtkFile << "Translated newpop data to dual grid for timestep " << i << " from netCDF by Christine Ahrens"
+ << endl;
+ vtkFile << "ASCII" << endl;
+ vtkFile << "DATASET UNSTRUCTURED_GRID" << endl;
+
+
+ // write points (the points are the primal-grid cell centers)
+
+ vtkFile << "POINTS " << numPoints << " double" << endl;
+
+ // first write a dummy point, because the climate code
+ // starts their cell numbering at 1 and VTK starts it at
+ // 0
+ vtkFile.precision(16);
+ vtkFile << (float)0.0 << "\t" << (float)0.0 << "\t" << (float)0.0
+ << endl;
+
+ for (int j = 0; j < nCells->size(); j++ )
+ {
+ vtkFile.precision(16);
+ if (abs(xCellData[j]) < 1e-126) xCellData[j] = 0;
+ if (abs(yCellData[j]) < 1e-126) yCellData[j] = 0;
+ if (abs(zCellData[j]) < 1e-126) zCellData[j] = 0;
+ vtkFile << (float)xCellData[j] << "\t" << (float)yCellData[j] << "\t"
+ << (float)zCellData[j] << endl;
+ }
+ vtkFile << endl;
+
+
+ // Write dual-mesh cells
+ // Dual-mesh cells are triangles with primal-mesh cell
+ // centers as the vertices.
+ // The number of dual-mesh cells is the number of vertices in the
+ // primal mesh.
+
+ vtkFile << "CELLS " << numCells << " "
+ << numCells * (THREE->size() + 1) << endl;
+
+ // for each dual-mesh cell
+ for (int j = 0; j < nVertices->size() ; j++) {
+
+ vtkFile << THREE->size() << "\t" ;
+
+ int* dualMeshCells = cellsOnVertex + (j * THREE->size());
+
+ for (int k = 0; k < THREE->size(); k++)
+ {
+ vtkFile << dualMeshCells[k] << "\t";
+ }
+
+ vtkFile << endl;
+ }
+ vtkFile << endl;
+
+
+ // write cell types -- all triangles
+
+ vtkFile << "CELL_TYPES " << numCells << endl;
+
+ for (int j = 0; j < numCells; j++)
+ {
+ vtkFile << VTK_TRIANGLE << endl;
+ }
+ vtkFile << endl;
+
+ // Write attributes of dual-mesh cell data (attributes of primal-mesh
+ // vertex data)
+
+ // for each var, figure out if it is by cell or by point
+
+
+ vtkFile.precision(16);
+
+ // If by point, write out point data
+
+ vtkFile << "POINT_DATA " << numPoints << endl;
+
+ for (int v = 0; v <= pointVarIndex; v++) {
+
+ vtkFile << "SCALARS " << pointVars[v]->name() << " float 1" << endl;
+
+ vtkFile << "LOOKUP_TABLE default" << endl;
+
+ double *var_target = pointVarData[v];
+
+ // write dummy
+ if (abs(*var_target) < 1e-126) *var_target = 0;
+ vtkFile << (float)*var_target << endl;
+
+ for (int j = 0; j < nCells->size(); j++)
+ {
+ if (abs(*var_target) < 1e-126) *var_target = 0;
+ vtkFile << (float)*var_target << endl;
+ var_target += nVertLevels->size();
+ }
+ }
+
+ // if by cell, then write out cell data
+
+ vtkFile << "CELL_DATA " << numCells << endl;
+
+ for (int v = 0; v <= cellVarIndex; v++) {
+
+ // write cell variable
+ vtkFile << "SCALARS " << cellVars[v]->name() << " float 1" << endl;
+
+ vtkFile << "LOOKUP_TABLE default" << endl;
+
+ double *var_target = cellVarData[v];
+
+ for (int j = 0; j < nVertices->size(); j++)
+ {
+ if (abs(*var_target) < 1e-126) *var_target = 0;
+ vtkFile << (float)*var_target << endl;
+ var_target += nVertLevels->size();
+ }
+ }
+
+ }
+
+ return(0);
+
+}
+
Added: branches/ocean/visualization/translator/vtkCellType.h
===================================================================
--- branches/ocean/visualization/translator/vtkCellType.h (rev 0)
+++ branches/ocean/visualization/translator/vtkCellType.h 2009-12-17 15:32:09 UTC (rev 86)
@@ -0,0 +1,98 @@
+/*=========================================================================
+
+ Program: Visualization Toolkit
+ Module: $RCSfile: vtkCellType.h,v $
+
+ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
+ All rights reserved.
+ See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notice for more information.
+
+=========================================================================*/
+// .NAME vtkCellType - define types of cells
+// .SECTION Description
+// vtkCellType defines the allowable cell types in the visualization
+// library (vtk). In vtk, datasets consist of collections of cells.
+// Different datasets consist of different cell types. The cells may be
+// explicitly represented (as in vtkPolyData), or may be implicit to the
+// data type (as in vtkStructuredPoints).
+
+#ifndef __vtkCellType_h
+#define __vtkCellType_h
+
+// To add a new cell type, define a new integer type flag here, then
+// create a subclass of vtkCell to implement the proper behavior. You
+// may have to modify the following methods: vtkDataSet (and subclasses)
+// GetCell() and vtkGenericCell::SetCellType(). Also, to do the job right,
+// you'll also have to modify some filters (vtkGeometryFilter...) and
+// regression tests (example scripts) to reflect the new cell addition.
+// Also, make sure to update vtkCellTypesStrings in vtkCellTypes.cxx.
+
+// .SECTION Caveats
+// An unstructured grid stores the types of its cells as a
+// unsigned char array. Therefore, the maximum encoding number for a cell type
+// is 255.
+
+typedef enum {
+ // Linear cells
+ VTK_EMPTY_CELL = 0,
+ VTK_VERTEX = 1,
+ VTK_POLY_VERTEX = 2,
+ VTK_LINE = 3,
+ VTK_POLY_LINE = 4,
+ VTK_TRIANGLE = 5,
+ VTK_TRIANGLE_STRIP = 6,
+ VTK_POLYGON = 7,
+ VTK_PIXEL = 8,
+ VTK_QUAD = 9,
+ VTK_TETRA = 10,
+ VTK_VOXEL = 11,
+ VTK_HEXAHEDRON = 12,
+ VTK_WEDGE = 13,
+ VTK_PYRAMID = 14,
+ VTK_PENTAGONAL_PRISM = 15,
+ VTK_HEXAGONAL_PRISM = 16,
+
+ // Quadratic, isoparametric cells
+ VTK_QUADRATIC_EDGE = 21,
+ VTK_QUADRATIC_TRIANGLE = 22,
+ VTK_QUADRATIC_QUAD = 23,
+ VTK_QUADRATIC_TETRA = 24,
+ VTK_QUADRATIC_HEXAHEDRON = 25,
+ VTK_QUADRATIC_WEDGE = 26,
+ VTK_QUADRATIC_PYRAMID = 27,
+ VTK_BIQUADRATIC_QUAD = 28,
+ VTK_TRIQUADRATIC_HEXAHEDRON = 29,
+ VTK_QUADRATIC_LINEAR_QUAD = 30,
+ VTK_QUADRATIC_LINEAR_WEDGE = 31,
+ VTK_BIQUADRATIC_QUADRATIC_WEDGE = 32,
+ VTK_BIQUADRATIC_QUADRATIC_HEXAHEDRON = 33,
+
+ // Special class of cells formed by convex group of points
+ VTK_CONVEX_POINT_SET = 41,
+
+ // Higher order cells in parametric form
+ VTK_PARAMETRIC_CURVE = 51,
+ VTK_PARAMETRIC_SURFACE = 52,
+ VTK_PARAMETRIC_TRI_SURFACE = 53,
+ VTK_PARAMETRIC_QUAD_SURFACE = 54,
+ VTK_PARAMETRIC_TETRA_REGION = 55,
+ VTK_PARAMETRIC_HEX_REGION = 56,
+
+ // Higher order cells
+ VTK_HIGHER_ORDER_EDGE = 60,
+ VTK_HIGHER_ORDER_TRIANGLE = 61,
+ VTK_HIGHER_ORDER_QUAD = 62,
+ VTK_HIGHER_ORDER_POLYGON = 63,
+ VTK_HIGHER_ORDER_TETRAHEDRON = 64,
+ VTK_HIGHER_ORDER_WEDGE = 65,
+ VTK_HIGHER_ORDER_PYRAMID = 66,
+ VTK_HIGHER_ORDER_HEXAHEDRON = 67,
+
+ VTK_NUMBER_OF_CELL_TYPES
+} VTKCellType;
+
+#endif
</font>
</pre>