# How to extract a contour line (or a region) from a ContourPlot?

Given a ContourPlot with a set of contours, say, this: 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"] 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?

## 5 Answers

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] • Michael, what is mesh in your code above? – Alexei Boulbitch Feb 9 '16 at 10:34
• @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. – Michael E2 Feb 9 '16 at 13:54
• 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) – Alexei Boulbitch Apr 10 '17 at 7:50
• (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? – Alexei Boulbitch Apr 10 '17 at 7:52
• @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. – Michael E2 Apr 10 '17 at 10:41

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[]] Show[Table[
RegionPlot[regiondiffs[[n]], PlotStyle -> Hue[n/8]], {n, 8, 1, -1}]] 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.

• 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. – Alexei Boulbitch Feb 9 '16 at 10:30
• @AlexeiBoulbitch - comment became too long so I just edited the answer. – Jason B. Feb 9 '16 at 10:39

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]}] 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]}
] • 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. – Alexei Boulbitch Feb 9 '16 at 10:54

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 For a line, simply put the value

ContourPlot[x*Exp[-x^2 - y^2] == 0.3, {x, 0, 2}, {y, -1, 1}]

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[], __}, -1],Polygon[__], -1];
r = RegionUnion[poly];
lines = Cases[Normal@RegionPlot[r], Line[__], -1];
Graphics[lines] 