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.

TernaryPlot.png

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")