# Get the coordinates from ContourPlot and RegionPlot

How do I get the coordinates from a contour plot I've done in Mathematica? For example, I have a two-variable function f[x, y], for which I can make a contour plot:

    contour = ContourPlot[f[x, y] == 1, {x, -1, 1}, {y, -1, 1}];


I can access a nested list, containing lists of coordinates, from an exported file using

    Export["output.m", contour, "TEXT"]


but the nested list is inside the function Graphics.

I would like to export only the corresponding nested list. What is a straightforward way to do that? Is there a way of doing that by manipulating directly the variable contour?

Edit

@Jens solution was

    Cases[Normal@contour, Line[x_] :> x, Infinity]


Is there a way to generalize it? The issue is that Line won't give us the points inside a region. E.g., I would like to get the coordinates corresponding to a region plot:

    region = RegionPlot[f[x, y] > 1, {x, -1, 1}, {y, -1, 1}];

-
Your edit is an entirely different question, it seems to me. What do you expect to get? A grid of points? How dense? –  Jens Feb 28 '13 at 0:34
Sorry, I thought that it would be wrong to open another post for a very similar question. I would like to get the coordinates associated to the region shown by a region plot. As dense as the RegionPlot output. –  fcpenha Feb 28 '13 at 1:26
Now I understood what you mean (I think), so I updated the answer. –  Jens Feb 28 '13 at 2:53
@Jens, do you think I should change the title to Get the coordinates from ContourPlot and RegionPlot? –  fcpenha Feb 28 '13 at 9:17
Sure, that would be appropriate. –  Jens Feb 28 '13 at 15:19

Since the plot usually is a GraphicsComplex, the extraction is easiest if you first convert using Normal:

contour = ContourPlot[x^2 + y^2 == 1, {x, -1, 1}, {y, -1, 1}];

Cases[Normal@contour, Line[x_] :> x, Infinity]


This produces a list that shows the coordinates in the order they were drawn.

Explanation:

The contours in the plot are drawn using the Line command, which takes lists of points (or lists of lists of points). Usually, these points are all collected at the front of the ContourPlot output in the form of a GraphicsComplex, such that each point can later be addressed by using an index from within the Line commands. By applying Normal, these indexed points are moved to where they are actually used in the drawing part of the output. Normal@contour is the same as Normal[contour].

After that is done, we can look for all Line commands and find the coordinates in sequential order inside of them. This is done by using Cases, which selects parts of the expression that match a pattern. The pattern is specified here as Line[x_] where x_ is a "dummy variable" that gets defined whenever a Line was found, by replacing it with the contents of the line. The final step is to tell Cases that when it does find a value for x, to just output that without the wrapper Line.

This search is done throughout the whole plot, which is indicated by the Infinity level specification.

Update for RegionPlot

Extracting points from a plot using Cases can be generalized to situations where the points aren't inside a Line. In RegionPlot, for example, you may want to extract the points that form the mesh with which the region is filled. This filling is done by a polygonal tesselation, so we have to simply replace Line with Polygon:

region = Normal@RegionPlot[x^2 + y^2 <= 1, {x, -1, 1}, {y, -1, 1}];
pts = DeleteDuplicates@
Flatten[Cases[region, Polygon[x_] :> x, Infinity], 1];

Graphics[Point[pts]]


Here I had to add two more steps before being able to make the plot: first, Flatten is used to remove all nested levels except the ones grouping the coordinate tuples for each point. Then, I added DeleteDuplicates just to remove any shared vertices between polygons, so one and the same point isn't re-drawn redundantly.

-
Please, could you explain the syntax? –  fcpenha Feb 26 '13 at 18:29
OK, I've added an explanation and slightly simplified the expression. –  Jens Feb 26 '13 at 18:39
Thanks. What if ContourPlot has two or more functions as arguments (e.g., {f[x,y],{g[x,y]}})? How should we change the code above? –  fcpenha Feb 26 '13 at 18:54
I assume you mean two equations of the type {f==1, g==1}? It should work without modification, and the order of the lines lets you determine which equation they belong to. The output of Cases is a list of lists, so the structure of the point sets is preserved. –  Jens Feb 26 '13 at 18:59
Then just take the list as above, call it l, and access the first list as l[[1]], the second list as l[[2]] etc. If that doesn't work for you, it's always possible to generate two separate contour plots to be sure you capture the correct coordinates. –  Jens Feb 26 '13 at 19:14

Here the answer for @Luiz Roberto Meier, which was closed although it was no duplicate, but needed some deeper reflections.

First, ContourPlot doesn't show the complete solution.

This plot is complex-valued, therefore you may use Re and Im

So, the extraction of the plotpoints in ContourPlot won't be very useful, but as already suggested it might be done of course using Cases:

points = ContourPlot[{y^4 + 6 x^3*y + x^8 == 0},
{x, -5, 5}, {y, -5, 5}] // Normal // Cases[#, Line[$_] ->$, Infinity] &
ListPlot[points, Joined -> True, AspectRatio -> 1, PlotRange -> 5]


Basically the two functions can be solved by

Clear@"*"
f[x_] := x^3
g[x_] := Block[{y}, y /. Solve[{y^4 + 6 x^3*y + x^8 == 0 }, y]]

Plot[{f[x], g[x]}, {x, -5, 5}, PlotRange -> 5, AspectRatio -> 1]


f is easy, so we concentrate on g from now on. Inspecting g we find that it has 4 complex solutions in y, (order 4)

Plotting g using Plot having a look at Re and Im gives:

Clear@"*"
g[x_] := Block[{y}, y /. Solve[{y^4 + 6 x^3*y + x^8 == 0 }, y]]
Plot[{g[x][[1]], g[x][[2]], g[x][[3]], g[x][[4]]} // Re //
Evaluate, {x, -5, 5}, PlotRange -> 5, AspectRatio -> 1]
Plot[{g[x][[1]], g[x][[2]], g[x][[3]], g[x][[4]]} // Im //
Evaluate, {x, -5, 5}, PlotRange -> 5, AspectRatio -> 1]


Let's say, we are interested in the real part. In the plot we see that we are missing some points for x near -2.5 As we are interested in a complete point list, we try increasing option values like MaxRecursions or PlotPoints (without success).

Changing strategy using ListPlot, where we can increment ourself the x-Range provides a high quality solution and a complete set of {x,y} point coordinates by analytical means:

Clear@"*"
g[x_] := Block[{y}, y /. Solve[{y^4 + 6 x^3*y + x^8 == 0 }, y]]
sol1 = Re@{#, g[#][[1]]} & /@ Range[-5, 5, .005];
sol2 = Re@{#, g[#][[2]]} & /@ Range[-5, 5, .005];
sol3 = Re@{#, g[#][[3]]} & /@ Range[-5, 5, .005];
sol4 = Re@{#, g[#][[4]]} & /@ Range[-5, 5, .005];
ListPlot[#, Joined -> True, PlotRange -> 5,
AspectRatio -> 1] & /@ {sol1, sol2, sol3, sol4}
ListPlot[#, Joined -> True, PlotRange -> 5, AspectRatio -> 1] &@{sol1,
sol2, sol3, sol4}


Of course it would be possible to map over the 4 solutions at once but for the sake of explanation I didn't.

the first image shows the solutions separated, the next image combines them.

all points for the complete real solution you were looking for g[x]`, are now in:

{ sol1, sol2, sol3, sol4 }

-