22
$\begingroup$

Given a ContourPlot with a set of contours, say, this:

enter image description here

is it possible to get the contours separating domains with the different colors in the form of lists?

For example, how to extract the boundaries of the blue domain in the image above? Or just for the sake of trial, from such a simple example:

 ContourPlot[x*Exp[-x^2 - y^2], {x, 0, 3}, {y, -3, 3}, 
 PlotRange -> {0, 0.5}, ColorFunction -> "Rainbow"]

enter image description here

The same task, let us find the lists corresponding to the blue domain boundaries.

To make it clear, I am not asking of how to get the lines from the function behind. This I understand. I ask of how to extract the contour lines that are generated by Mma.

Let us put this question another way around. Is it possible to define the areas with the same color as separate geometric regions in the sense of the computation geometry, and then work with these domains separately?

Improve this question
$\endgroup$

5 Answers 5

Reset to default
13
$\begingroup$

To answer the last question, the contour domains (since V8) are enclosed separately in GraphicsGroup, each which you can cull and turn into a region:

plot = ContourPlot[x*Exp[-x^2 - y^2], {x, 0, 3}, {y, -3, 3}, 
  PlotRange -> {0, 0.5}, ColorFunction -> "Rainbow"];

regs = With[{coords = First@Cases[plot, GraphicsComplex[p_, ___] :> p, Infinity]},
  BoundaryDiscretizeGraphics@
    GraphicsComplex[coords, #] & /@ Cases[plot, _GraphicsGroup, Infinity]
  ];

Multicolumn[regs, 5]

Mathematica graphics

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ Michael, what is mesh in your code above? $\endgroup$
    – Alexei Boulbitch
    Feb 9, 2016 at 10:34
  • $\begingroup$ @AlexeiBoulbitch Oops, I changed the name of the plot from mesh to plot, but forgot one reference I guess. Fixed now. Thanks for pointing it out. $\endgroup$
    – Michael E2
    Feb 9, 2016 at 13:54
  • $\begingroup$ Michael, please have a look. I am trying your above approach on the following case: trapezoid = Polygon[{{0, 0}, {0.05, 0.15}, {0.1, 0.15}, {0.15, 0}}]; plot = ContourPlot[(1 - 3.38*y + 265.*y^2)*(1 + (150 - 2000*y)*(x - 0.075)^2), {x, y} \[Element] trapezoid] and then I try to apply your function regs to it. The calculation takes infinity and I never can wait until its end. If, in contrast, I make first a list out of the function: (see the continuation) $\endgroup$
    – Alexei Boulbitch
    Apr 10, 2017 at 7:50
  • $\begingroup$ (Continuation) d = 0.005; plot2 = ListContourPlot[Flatten[Table[If[RegionMember[trapezoid, {i*d, j*d}], {i*d, j*d, (1 - 3.38*j*d + 265.*(j*d)^2)*(1 + (150 - 2000*j*d)*(i*d - 0.075)^2)}, Nothing], {i, 1, 0.15/d}, {j, 1, 0.15/d}], 1]]; I can apply regs, and its evaluation takes a few seconds. I use Mma 11.1, Win7. Have you an idea, why this happens? $\endgroup$
    – Alexei Boulbitch
    Apr 10, 2017 at 7:52
  • $\begingroup$ @AlexeiBoulbitch Examining them individually BoundaryDiscretizeGraphics quickly returns unevaluated on the first (little ovoid), works on parts 5 & 7, and stalls on the rest. There are some tiny components in the first (ovoid, lower left, on the edge). Maybe they overlap? -- Slightly more successful: With[{coords = First@Cases[plot, GraphicsComplex[p_, ___] :> p, Infinity]}, MeshRegion[coords, Polygon /@ Join @@ Cases[#, _Polygon, Infinity][[All, 1]]]] & /@ Cases[plot, _GraphicsGroup, Infinity] -- Not really obvious what the problem is. Some digging necessary, I guess. $\endgroup$
    – Michael E2
    Apr 10, 2017 at 10:41
17
$\begingroup$

I don't know how to do this in an automated way, but here is something at least:

Make your plot, extract the lines, convert them to regions, and then take the RegionDifference between them

plot = ContourPlot[x*Exp[-x^2 - y^2], {x, 0, 3}, {y, -3, 3}, 
  PlotRange -> {0, 0.5}, ColorFunction -> "Rainbow"]
points = Cases[Normal@plot, Line[pts__] -> pts, Infinity];
regions = BoundaryMeshRegion[#, Line[Range[Length@#]]] & /@ points;
regiondiffs = RegionDifference[#2, #1] & @@@ Partition[regions, 2, 1];
PrependTo[regiondiffs, regions[[1]]]

enter image description here

Show[Table[
  RegionPlot[regiondiffs[[n]], PlotStyle -> Hue[n/8]], {n, 8, 1, -1}]]

enter image description here

Edit @AlexeiBoulbitch - I had tried to extract the points without Normal at first, and the result was a list of integers rather than {x,y} coordinates. And Szabolcs had just shown me that Normal is useful when extracting regions from contour plots so I tried it.

What Normal is doing is removing the GraphicsComplex head:

GraphicsComplex[{$pt_1$, $pt_2$ ,...},$data$] represents a graphics complex in which coordinates given as integers i in graphics primitives in data are taken to be $pt_i$.  >>

So after using Normal, you get the actual coordinates for the lines.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you. Could you please kindly comment, what does Normal when acting on ContourPlot. Intuitively it seems to create a mesh and form the lists corresponding to the contours. Am I right? I could not find any explanation in the documentation. $\endgroup$
    – Alexei Boulbitch
    Feb 9, 2016 at 10:30
  • $\begingroup$ @AlexeiBoulbitch - comment became too long so I just edited the answer. $\endgroup$
    – Jason B.
    Feb 9, 2016 at 10:39
8
$\begingroup$

You ask whether it is "possible to define the areas with the same color as separate geometric regions in the sense of the computation geometry, and then work with these domains separately". This seems to me to be a great task for ImplicitRegion:

regions = Table[
 ImplicitRegion[i < x*Exp[-x^2 - y^2] <= i + 0.05, {{x, 0, 3}, {y, -3, 3}}],
 {i, 0, 0.4, 0.1}
]

(*Out:
{ImplicitRegion[0. < E^(-x^2 - y^2) x <= 0.05 && 0 <= x <= 3 && -3 <= y <= 3, {x, y}],
 ImplicitRegion[0.05 < E^(-x^2 - y^2) x <= 0.1 && 0 <= x <= 3 && -3 <= y <= 3, {x, y}],
 ImplicitRegion[0.1 < E^(-x^2 - y^2) x <= 0.15 && 0 <= x <= 3 && -3 <= y <= 3, {x, y}],
 ImplicitRegion[0.15 < E^(-x^2 - y^2) x <= 0.2 && 0 <= x <= 3 && -3 <= y <= 3, {x, y}],
 ImplicitRegion[0.2 < E^(-x^2 - y^2) x <= 0.25 && 0 <= x <= 3 && -3 <= y <= 3, {x, y}],
 ImplicitRegion[0.25 < E^(-x^2 - y^2) x <= 0.3 && 0 <= x <= 3 && -3 <= y <= 3, {x, y}],
 ImplicitRegion[0.3 < E^(-x^2 - y^2) x <= 0.35 && 0 <= x <= 3 && -3 <= y <= 3, {x, y}],
 ImplicitRegion[0.35 < E^(-x^2 - y^2) x <= 0.4 && 0 <= x <= 3 && -3 <= y <= 3, {x, y}],
 ImplicitRegion[0.4 < E^(-x^2 - y^2) x <= 0.45 && 0 <= x <= 3 && -3 <= y <= 3, {x, y}]}
*)

These regions can be used for further calculation or plotting. For instance we can select a random point in one of those regions:

region=ImplicitRegion[0.2 < E^(-x^2 - y^2) x <= 0.25 && 0 <= x <= 3 && -3 <= y <= 3, {x, y}];
pt = RandomPoint[region];
RegionPlot[region, Epilog -> {PointSize[0.02], Red, Point[pt]}]

Mathematica graphics

Or we can plot each one of those regions to reproduce the ContourPlot:

Show@Table[
    RegionPlot[regions[[i]], PlotStyle -> ColorData["Rainbow"][(i - 1)/Length[regions]]], 
    {i, Length[regions]}
  ]

Mathematica graphics

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Thank you, Marco, this is, however, not what I had in mind. The function x*Exp[-x^2 - y^2] above was only given for the illustrative purposes. I would like to act, as if this function is unknown and cannot be used. I will then apply it to cases, in which such a function is not known, but the ContourPlot is available. $\endgroup$
    – Alexei Boulbitch
    Feb 9, 2016 at 10:54
5
$\begingroup$

A simple way could be using the PlotRange.

Table[
  ContourPlot[x*Exp[-x^2 - y^2], {x, 0, 2}, {y, -1, 1}, 
   PlotRange -> {i, j}, ColorFunction -> "Rainbow", 
   PlotLegends -> True, PlotLabel -> {i, j}]
   , {i, 0.1, 0.3, 0.1}, {j, i + 0.1, 0.4, 0.1}] // TableForm

enter image description here

For a line, simply put the value

ContourPlot[x*Exp[-x^2 - y^2] == 0.3, {x, 0, 2}, {y, -1, 1}]
Improve this answer
$\endgroup$
5
$\begingroup$ 
p = ContourPlot[x*Exp[-x^2 - y^2], {x, 0, 3}, {y, -3, 3}, 
  PlotRange -> {0, 0.5}, ColorFunction -> "Rainbow"];

colors = Cases[p, _RGBColor, -1];
poly = Cases[Cases[Normal@p, {__, colors[[2]], __}, -1],Polygon[__], -1];
r = RegionUnion[poly];
lines = Cases[Normal@RegionPlot[r], Line[__], -1];
Graphics[lines]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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