Ternary Plot
With a little scripting and data reorganization, VisIt can produce different types of plots using the built-in plots as building blocks. For instance, VisIt can be used to create a ternary plot.
Converting Data
The following example program can be used to convert some comma-separated values files into data suitable for plotting as a ternary plot.
The input data file would consist of comma-separated values like this:
fat,protein,lactose,benefit/cost,profit in
0,100,39,12,93
0,90,45,13,95
0,80,51,13,97
Here is the maketernary.py Python source, which you would run in the VisIt cli:
visit -cli -s maketernary.py
import math, string, sys
# Add to the Python path so Python can locate visit_writer.
sys.path.append("/usr/local/apps/visit/2.5.0/darwin-x86_64/lib")
import visit_writer
def ReadCSV(filename, skipCount):
f = open(filename, "rt")
lines = f.readlines()
values = []
for line in lines[skipCount:]:
v = [float(x) for x in string.split(line, ",")]
if values == []:
for i in xrange(len(v)):
values = values + [[]]
for i in xrange(len(v)):
values[i] = values[i] + [v[i]]
return values
def WriteTernaryData(filename, axis1, axis2, axis3, value, varname, side, ndivisions):
def normalize(arr):
min = arr[0]
max = arr[0]
for v in arr[1:]:
if v < min:
min = v
if v > max:
max = v
norm = []
delta = max - min
if delta == 0.:
delta = 1.
for v in arr:
norm = norm + [(v - min) / delta]
return norm
def shepardglobal(pt, coords, values):
# Find the value using global shepard's method. This is a crude
# method for converting values at discrete points to a continuous
# field.
sum_di2inv = 0.
sum_fi_di2inv = 0.
for i in xrange(len(coords)):
dX = pt[0] - coords[i][0]
dY = pt[1] - coords[i][1]
try:
di2inv = 1. / (dX*dX + dY*dY);
except ZeroDivisionError:
di2inv = 1. / 1.e-6;
fi_di2inv = values[i] * di2inv;
sum_di2inv = sum_di2inv + di2inv
sum_fi_di2inv = sum_fi_di2inv + fi_di2inv
return sum_fi_di2inv / sum_di2inv
b1 = normalize(axis1)
b2 = normalize(axis2)
b3 = normalize(axis3)
# Enforce b1+b2+b3=1
sums = []
for i in xrange(len(b1)):
sum = b1[i]+b2[i]+b3[i]
b1[i] = b1[i] / sum
b2[i] = b2[i] / sum
b3[i] = b3[i] / sum
A1 = (-side, 0.)
A2 = (side, 0.)
A3 = (0., side * math.sqrt(3.))
x = []
y = []
# Convert barycentric to cartesian
for i in xrange(len(b1)):
x0 = b1[i]*A1[0] + b2[i]*A2[0] + b3[i]*A3[0]
y0 = b1[i]*A1[1] + b2[i]*A2[1] + b3[i]*A3[1]
x = x + [x0]
y = y + [y0]
# Make points for a triangle mesh.
a3a1 = []
a3a2 = []
for i in xrange(ndivisions+1):
t = float(i) / float(ndivisions)
a3a1 = a3a1 + [((1.-t)*A3[0] + t*A1[0], (1.-t)*A3[1] + t*A1[1])]
a3a2 = a3a2 + [((1.-t)*A3[0] + t*A2[0], (1.-t)*A3[1] + t*A2[1])]
points = []
id = 0
ids = []
for i in xrange(ndivisions+1):
nnodes_this_row = i+1
nodes = []
tids = []
for j in xrange(nnodes_this_row):
t = 0.
if nnodes_this_row > 1:
t = float(j) / float(nnodes_this_row-1)
xc = (1.-t)*a3a1[i][0] + t*a3a2[i][0]
yc = (1.-t)*a3a1[i][1] + t*a3a2[i][1]
nodes = nodes + [(xc,yc)]
tids = tids + [id]
id = id + 1
points = points + nodes
ids = ids + [tids]
# Connect the points into triangles
tris = []
for row in xrange(ndivisions):
ntri_this_row = row*2 + 1
tri_row = {}
r0 = ids[row]
r1 = ids[row+1]
# Even triangles
for i in xrange(len(r0)):
rid = i*2
a = r0[i]
b = r1[i]
c = r1[i+1]
tri_row[rid] = (a,b,c)
# Odd triangles
if row > 0:
nodds = len(r0)-1
for i in xrange(len(r0)-1):
rid = i*2+1
a = r0[i]
b = r1[i+1]
c = r0[i+1]
tri_row[rid] = (a,b,c)
keys = sorted(tri_row.keys())
for k in keys:
tris = tris + [tri_row[k]]
# Sample the data onto the triangle mesh.
samples = []
for p in points:
samples = samples + [shepardglobal(p, zip(x,y), value)]
# Write the data to a file
points3d = []
for p in points:
points3d = points3d + [p[0], p[1], 0.]
connectivity = []
for t in tris:
connectivity = connectivity + [(visit_writer.triangle, t[0], t[1], t[2])]
vars = ((varname, 1, 1, samples),)
visit_writer.WriteUnstructuredMesh(filename, 0, points3d, connectivity, vars)
# Write the original data as a point mesh that we can overlay to check the values
points3d = []
for i in xrange(len(x)):
points3d = points3d + [x[i],y[i],0.]
vars = ((varname, 1, 1, value),)
visit_writer.WritePointMesh("points_"+filename, 0, points3d, vars)
# Read the data
values = ReadCSV("input.csv", 1)
# Convert it and write it back out as VTK data
WriteTernaryData("benefit.vtk", values[0], values[1], values[2], values[3], "benefit", 1., 70)
WriteTernaryData("profit.vtk", values[0], values[1], values[2], values[4], "profit", 1., 70)
Plotting Data
The following example program can be used to set up a ternary plot in VisIt. It creates some annotation objects to label the axes instead of relying on VisIt's 2D axes.
import math
def PlotTernaryData(filename, axis1, axis2, axis3, value):
def MakeAxis(A, B, vec, nsteps):
for i in xrange(nsteps):
t = float(i) / float(nsteps-1)
x = (1.-t)*A[0] + t*B[0] + vec[0]
y = (1.-t)*A[1] + t*B[1] + vec[1]
label = CreateAnnotationObject("Text3D")
label.text = str(int(t * 100)) + "%"
label.position = (x,y,0)
label.relativeHeight = 0.01
OpenDatabase(filename)
AddPlot("Pseudocolor", value)
DisableRedraw()
# Turn off annotations
a = GetAnnotationAttributes()
a.axes2D.visible = 0
a.userInfoFlag = 0
SetAnnotationAttributes(a)
DrawPlots()
defaultTitle = "3D text annotation"
char2len = 1.1 / float(len(defaultTitle))
side = 1.
A1 = (-side, 0.)
A2 = (side, 0.)
A3 = (0., side * math.sqrt(3.))
d = 0.1
a1 = CreateAnnotationObject("Text3D")
a1.text = axis1
a1.rotations = (0,0,60)
p1 = (-0.5 - d*math.cos(math.pi/6.), math.sqrt(3.)/2. + d*math.sin(math.pi/6.))
v1 = (-0.5, -math.sqrt(3.)/2.)
len1 = len(axis1)
a1.position = (p1[0] + v1[0] * len1 * char2len/2., p1[1] + v1[1] * len1 * char2len/2., 0.)
a2 = CreateAnnotationObject("Text3D")
a2.text = axis2
p2 = (0., -0.12)
v2 = (-1., 0.)
len2 = len(axis2)
a2.position = (p2[0] + v2[0] * len2 * char2len/2., p2[1] + v2[1] * len2 * char2len/2., 0.)
a3 = CreateAnnotationObject("Text3D")
a3.text = axis3
a3.rotations = (0,0,-60)
p3 = (0.5 + d*math.cos(math.pi/6.), math.sqrt(3.)/2. + d*math.sin(math.pi/6.))
v3 = (-0.5, math.sqrt(3.)/2.)
len3 = len(axis3)
a3.position = (p3[0] + v3[0] * len3 * char2len/2., p3[1] + v3[1] * len3 * char2len/2., 0.)
MakeAxis(A2, A3, (0,0), 5)
MakeAxis(A3, A1, (-0.07, 0), 5)
MakeAxis(A1, A2, (0., -0.04), 5)
v = GetView2D()
v.windowCoords = (-1.12133, 1.12133, -0.335184, 1.83518)
v.viewportCoords = (0, 1, 0, 1)
SetView2D(v)
RedrawWindow()
# Plot the data from the benefit.vtk file.
PlotTernaryData("benefit.vtk", "fat", "protein", "lactose", "benefit")
