Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

How can I get a ternary density plot just like the plots from OriginLab? enter image description here

ContourPlot and DensityPlot seemingly can solve the [f,{x},{y}]-style,but cannot solve the [f,{x},{y},{z}]-style

share|improve this question
3  
Here is an example on how you can reduce three coordinates into two to make a ternary list plot (and contour plot in comments): mathgis.blogspot.com/2008/12/how-to-make-tenary-plot.html –  shrx Jan 3 at 10:46
add comment

4 Answers

A basic approach

You can start with a regular density plot, restricted to the domain of x and y using RegionFunction. Then you can transform the plot to an equilateral triangle.

f[p_, q_, r_] := r Sin[10 p]^2 + (1 - r) r Cos[20 q]^2;

dp = DensityPlot[f[x, y, 1 - x - y], {x, 0, 1}, {y, 0, 1}, 
  RegionFunction -> (#1 <= 1 - #2 &), ColorFunction -> (Hue[0.85 #] &),
   Frame -> None, PlotRange -> All, AspectRatio -> Automatic]

Mathematica graphics

It's easy enough to construct the transformation by hand, but it's also easy to use FindGeometricTransform.

{error, xf} = 
  FindGeometricTransform[
    {{0, 0}, {1, 0}, {1, Tan[Pi/3]}/2},
    {{0, 0}, {1, 0}, {0, 1}}];

Graphics[
  GeometricTransformation[
   First @ dp,
   xf
   ]]

Mathematica graphics

We can apply ticks modifying this answer to create a triangleTicks function (see below).

triangleTicks[Graphics[
  GeometricTransformation[
   First@dp,
   xf
   ]]
 ]

Mathematica graphics


Update 3 - A better looking (at low-res) approach

Here's another parametrization that goes along with Mathematica's native subdivision of the plot domain. It shows that the right triangles in the subdivision of the base image are mapped to equilateral triangles in the transformed image. So it looks less distorted, although from a mathematical point of view, the denser the same points, the more faithful the representation. The plot above appears to have a roughly ENE bias due to the mesh subdivision getting compressed.

cartestianToBarycentric2 = 
 Compile[{{x, _Real}, {y, _Real}}, (x - y) {1, 0} + y {1, Sqrt[3]}/2, 
  RuntimeAttributes -> {Listable}, RuntimeOptions -> "Speed"]; base = 
 DensityPlot[f[x - y, y, 1 - x], {x, 0, 1}, {y, 0, 1},
  Mesh -> All, MaxRecursion -> 0, RegionFunction -> (#1 >= #2 &), 
  ColorFunction -> (Hue[0.85 #] &), Frame -> None, PlotRange -> All, 
  AspectRatio -> Automatic];
transformed = MapAt[
   cartestianToBarycentric2 @@ Transpose[#] &,
   base,
   {1, 1}];
Graphics[{
  Inset[Show[base, AlignmentPoint -> {1.2, 0}], {0, 0}, Automatic, 
   1],
  Thick, Arrow[{{-0.1, 0.4}, {0.15, 0.4}}],
  Inset[Show[transformed, AlignmentPoint -> {-0.1, 0}], {0, 0}, 
   Automatic, 1]},
 PlotRange -> {{-1.2, 1.15}, {-0.05, 1.0}}
 ]

Mathematica graphics

Here is a plot with a coordinate grid:

dp = DensityPlot[f[x - y, y, 1 - x], {x, 0, 1}, {y, 0, 1},
   MeshFunctions -> {#1 - #2 &, #2 &, 1 - #1 &}, Mesh -> 19,
   RegionFunction -> (#1 >= #2 &), ColorFunction -> (Hue[0.85 #] &), 
   Frame -> None, PlotRange -> All, AspectRatio -> Automatic];
triangleTicks[
 MapAt[
  cartestianToBarycentric2 @@ Transpose[#] &,
  dp,
  {1, 1}]
 ]

Mathematica graphics


Update 1 -- Getting rid of GeometricTransform from the plot

Alexey Popkov pointed out that there is a problem with GeometricTransform and exporting. It was easy to fix the density plot, and the ticks needed rewriting (see below).

cartestianToBarycentric = 
  Compile[{{x, _Real}, {y, _Real}}, x {1, 0} + y {1, Sqrt[3]}/2, 
   RuntimeAttributes -> {Listable}, RuntimeOptions -> "Speed"];

MapAt[
 cartestianToBarycentric @@ Transpose[#] &,
 dp,
 {1, 1}]

Update 2 -- Getting rid of GeometricTransform from the ticks

tickGraphics now avoids using GeometricTransform

Basically tickGraphics creates an axis on one side of the equilateral triangle. It is rotated about the sides (counter-rotating the text). The argument total represents the sum of the ternary variables (which should be constant).

ClearAll[tickGraphics, triangleTicks];
Module[{ticklen = 0.025, (*length of ticks (const parameter)*)
        textoffset = 0.04},

 tickGraphics[tickrange : {0., total_}, angle_] := 
  With[{rotSideXF = RotationTransform[angle, total {1., Tan[Pi/6]}/2.]},
   {Line[rotSideXF /@ Thread[{tickrange, 0}]],
    With[{rotTickXF = Composition[
          RotationTransform[π/3., rotSideXF @ {#, 0.}],
          rotSideXF]},
       {Text[NumberForm[N @ #, {3, 2}], 
         rotTickXF[total {-ticklen - textoffset, 0} + {#, 0.}], {0, 0.}],
        {Line[rotTickXF /@ {total {-ticklen, 0} + {#, 0.}, {#, 0.}}]}}] & /@
     Select[N @ FindDivisions[tickrange, 4], 0. <= # <= total &]
    }
   ];

 triangleTicks[g_Graphics, total_: 1] :=
  Show[
   Graphics[First@g],
   Graphics[{
     AbsoluteThickness[0.5],
     Table[
      tickGraphics[{0., total}, N@angle],
      {angle, 0, 4 Pi/3, 2 Pi/3}]}],
   Axes -> False, Frame -> None, 
   PlotRange -> total {{0, 1 + ticklen}, {0, Sqrt[3]/2 + ticklen/2}}, 
   PlotRangePadding -> 0.07 total, PlotRangeClipping -> False, 
   ImagePadding -> {{20, 5}, {20, 5}}]
 ]
share|improve this answer
    
What is the definition for function f[x, y, 1 - y] you have plotted? –  Alexey Popkov Jan 3 at 17:54
    
Ugh, it took me forever to find the transform, +1. –  bobthechemist Jan 3 at 18:09
    
Is the range you've used for f correct for a ternary plot? I thought the sum of the arguments needs to be <= 1. –  bobthechemist Jan 3 at 18:19
    
According to Wikipedia, the sum of three variables must be a constant for ternary plot. –  Alexey Popkov Jan 3 at 18:25
1  
The answer is great. It is worth to mention that if someone wants ticks to be perpendicular to the corresponding side of the equilateral triangle it is sufficient to replace π/3. with π/2. in the tickGraphics function. (+1) –  Alexey Popkov Jan 4 at 9:55
show 6 more comments

Update: now it is real ternary plot.

You can start with the 2D-adaptation of the surface plotting:

texturize[f_, n_, colf_] := # /. Polygon[{v1_, v2_, v3_}] :> {EdgeForm[], 
      Texture@ImageData@Colorize[Image@f[#1, #2, 1 - #1 - #2] &[#, Transpose[#]] &@
          ConstantArray[Range[-1./n, 1 + 1/n, 1./n], n + 3], 
         ColorFunction -> colf, ColorFunctionScaling -> False], 
      Polygon[{v1, v2, v3}, VertexTextureCoordinates -> {{1 - 1.5/(n+3), 1 - 1.5/(n+3)}, 
         {1.5/(n+3), 1.5/(n+3)}, {1.5/(n+3), 1 - 1.5/(n+3)}}]} &;

f = #3 Sin[10 #1]^2 + (1 - #3) #3 Cos[20 #2]^2 &;

triangleTicks@Graphics@N@Polygon[{{0, 0}, {1, 0}, {1/2, Sqrt[3]/2}}] 
     // texturize[f, 300, Hue]

enter image description here

Here f assumed to be Listable. For ticks I used Michael's triangleTicks function.

Note, that this approach is ~100 times faster than corresponding DensityPlot!

share|improve this answer
4  
Is this really a ternary plot? Or just the plot simply mapped/clipped to a triangle? In a ternary plot, three variables are plotted and at each point in the triangle, the sum of the three variables is a constant. –  rm -rf Jan 3 at 14:38
    
@rm-rf I have found time to update my post. Now it is really a ternary plot. –  ybeltukov Jan 4 at 14:31
    
(+1) Your code is fast but very difficult to read. Please add an explanation what it basically does. –  Alexey Popkov Jan 4 at 15:44
    
@AlexeyPopkov It replaces the triangle to the triangle with a texture. I prepare a rectangular grid with ConstantArray and calculate texture on it. Values 1./n and 1.5/n are necessary to produce small margins in the texture. –  ybeltukov Jan 4 at 17:41
    
Why not do all DensityPlots this way (with an appropriate mapping of a rectangle, etc.)? –  Michael E2 Jan 4 at 19:27
show 1 more comment

I chose to start my solution from ternary points as opposed to a plot in cartesian points. My concern with the latter route is that 1 of the axes (the hypotenuse of the original triangle) is larger than the x and y axes and therefore is transformed differently from these two axes.

First, the meat of the transformation; I find the transformation function to convert a ternary point to a point on the {{1,0},{0,1},{0.5,Sqrt[3]/2}} triangle and create a grid ternary points {a,b,c} that obey a + b + c = 1.

(* find the transformation function *)
tf = FindGeometricTransform[{{0, 0, 0}, {1, 0, 0}, {0.5, Sqrt[3]/2,0}}, 
  {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}][[2]];
(* function to make cartesian coordinates out of ternary coordinates *)
ccoords[pt_] := tf[pt][[{1, 2}]];
(* Create a set of coordinates that total 1 *)
tcoords = 
  DeleteCases[Tuples[Range[0, 100, 1], 3], x_ /; Total[x] != 100]/100;

Next I make several functions to create tickmarks. These are a bit cumbersome because I first tried to make them general (any number of tickmarks) but backed away from that idea once I realized it was over my head:

Clear[tickPoints, tickMarks, tickLabels, axesLabels, region]
(* Now assuming 10 tick marks, not debugged with num =!= 10 *)
tickPoints[poly_, num_: 10] := Flatten[Map[
    (* Segment a line *) 
    Table[
      poly[[#[[1]]]] + 1/num  i (poly[[#[[2]]]] - poly[[#[[1]]]]), {i,
        1, num}] &,
    (* Map over all faces, uses same number of ticks per face, 
    not equidistant segments *)
    Table[{i, i + 1}, {i, 1, 2}]~Join~{{3, 1}}], 1];
(* Make tick marks angled such that the point along the ternary axes *)


tickMarks[poly_, num_: 10] := (GeometricTransformation[
      Line[#], ScalingTransform[0.25, #[[2]] - #[[1]], #[[1]]]] & /@ 
    Quiet@Drop[
      With[{list = tickPoints[region, 10]},
       MapIndexed[{#1,

          Flatten[RotationTransform[-60 Degree, #1][
            list[[#2 + 1]]]]} &, list]
       ], {10, -1, 10}]);

(* Create tick labels for the ticks - note I do not use tick marks on \
the vertices *)
tickLabels[num_: 10] := Module[{t},
   t = Quiet@Drop[
      With[{list = tickPoints[region, num]}, 
       MapIndexed[
         Flatten[RotationTransform[-90 Degree, #1][
            tickPoints[region][[#2 + 1]]]] &, tickPoints[region]][[
        1 ;; -2]]], {num, -1, num}];
   MapThread[
    Text[N@#1, #2] &, {Flatten@
      ConstantArray[Range[0 + 1/num, 1 - 1/num, 1/num], 3], t}]];
(* Instead of vertex tick marks, I use a separate function and call \
these labels *)
axesLabels[a_, b_, c_] := {Text[a, {-0.01, 0}, {1, 0}], 
   Text[b, {1.01, 0}, {-1, 0}], Text[c, {0.5, 1.02 Sqrt[3]/2}]};
(* These functions require that a polygon is defined *)
region = {{0, 0}, {1, 0}, {0.5, Sqrt[3]/2}, {0, 0}};

Some trivial functions to show that the plotting makes sense:

tfunc[a_, b_, c_] := a;

ListDensityPlot[
  Partition[Flatten@Transpose[{ccoords /@ tcoords, tfunc @@@ tcoords}], 3],
    ColorFunction -> (Hue[0.85 #] &),
    Epilog -> {Line@region, Black, Quiet@tickMarks[region, 10], tickLabels[], 
      axesLabels["A", "B", "C"]}, Frame -> None, 
    PlotRange -> {{-0.15, 1.15}, {-0.15, 1.15}}, 
    PlotLegends -> Automatic]

Mathematica graphics

tfunc[a_, b_, c_] := Sin[Pi a] + Sin[2 Pi b] + Sin[3 Pi c];

Mathematica graphics

I haven't done any thorough testing at this point; however I can Export these images without problem.

share|improve this answer
add comment

First transform cartesian coordinates to simplex coordinates and then apply the function f to get the $z$-values:

f[{a_, b_, c_}] := Which[ (* simple z-values for testing *)
       a >= 1/2, 0,
       b >= 1/2, 1,
       c >= 1/2, 2,
       True, 3]; 

(* transform simplex coordinates to cartesian ones *)
(* using the following simplex vertices: {{0, 0}, {1, 0}, {1/2, Sqrt[3]/2}} *)
transform[{x_, y_}] := {1, 0, 0} + x*{-1, 1, 0} + y*{-1/Sqrt@3, -1/Sqrt@3, 2/Sqrt@3};

(* test whether a point is inside the simplex *)
insideQ[pt_List] := (Total@pt==1 && And@@NonNegative@pt);

ContourPlot[f@transform@{x, y}, {x, 0, 1}, {y, 0, 1}, 
   RegionFunction -> (insideQ@transform@{#1, #2} &), 
   BoundaryStyle -> Black,
   MaxRecursion -> 3]

Mathematica graphics

Using the same f function as Michael E2:

f[{p_, q_, r_}] := r Sin[10 p]^2 + (1 - r) r Cos[20 q]^2;
{
 DensityPlot[f@{x, y, 1 - x - y}, {x, 0, 1}, {y, 0, 1}, 
  ColorFunction -> (Hue[0.85 #] &), 
  RegionFunction -> (#1 <= 1 - #2 &)],
 DensityPlot[f@transform@{x, y}, {x, 0, 1}, {y, 0, 1}, 
  ColorFunction -> (Hue[0.85 #] &), 
  RegionFunction -> (insideQ@transform@{#1, #2} &), 
  BoundaryStyle -> Black]
 }

Mathematica graphics

share|improve this answer
    
(+1) For comparison with Michael's answer one should define f as f[{p_, q_, r_}] := r Sin[10 p]^2 + (1 - r) r Cos[20 q]^2;. –  Alexey Popkov Jan 4 at 10:07
    
Thanks @Alexey, I've included it. –  István Zachar Jan 4 at 10:44
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.