How to get contour levels of picked automatically by ContourPlot?

By default, contours generated by ContourPlot have Tooltips with their numerical value as the label. You can use Cases to grab the values programmatically.

f[x_, y_] := Log[(1 - x)^2 + 100 (y - x^2)^2];
plot = ContourPlot[f[x, y], {x, -2, 0.5}, {y, -2, 2}, 
  ContourStyle -> {{Gray, Thick}}, PlotPoints -> 100, 
  ColorFunction -> Function[GrayLevel[1 - .45 #]], Frame -> None]

levels = Union@Cases[plot, Tooltip[_, label_] :> label, Infinity];

Plot3D[f[x, y], {x, -2, 0.5}, {y, -2, 2}, PlotPoints -> 100, 
 MeshFunctions -> {#3 &}, Mesh -> {levels}, Ticks -> None]

enter image description here


ContourPlot[f[x, y], {x, -2, .5}, {y, -2, 2}, 
   Contours -> ((contours = FindDivisions[{#, #2}, #3, Method -> {}]) &)];
Plot3D[f[x, y], {x, -2, .5}, {y, -2, 2}, PlotPoints -> 100, 
  MeshFunctions -> {#3 &}, Mesh -> {contours}, Ticks -> None]

enter image description here

Jason's approach, extracting Tooltip labels form ContourPlot output, is straightforward and convenient.

For the curious, how ContourPlot picks a nice set of curves as level curves by default remains a puzzle.

Speaking of "nice", searching for "nice numbers" in the Documentation Center gives: enter image description here

See also: Google Search: wolfram documentation "nice" numbers.

Documentation >> FindDivisions: enter image description here enter image description here Documentation >> Contours: enter image description here

After a few trials, we find that the function FindDivisions[{#, #2}, #3, Method -> {}]& (where the first two arguments are $z_{min}$ and $z_{max}$ and the last argument is the number of divisions which is fixed at 10) as the setting for Contours gives the same output as the default setting for this option.

Some examples:

ContourPlot[Cos[x] + Cos[y], {x, 0, 4 Pi}, {y, 0, 4 Pi}] == 
 ContourPlot[Cos[x] + Cos[y], {x, 0, 4 Pi}, {y, 0, 4 Pi}, 
  Contours -> (FindDivisions[{#, #2}, #3, Method -> {}] &)] ==
 ContourPlot[Cos[x] + Cos[y], {x, 0, 4 Pi}, {y, 0, 4 Pi}, 
  Contours -> (FindDivisions[{#, #2}, 10, Method -> {}] &)]

True

rr = RandomReal[5, 2];
ContourPlot[Evaluate[Sum[Sin[rr.{x, y}], {5}]], {x, 0, 5}, {y, 0, 5}] == 
 ContourPlot[Evaluate[Sum[Sin[rr.{x, y}], {5}]], {x, 0, 5}, {y, 0, 5},
   Contours -> (FindDivisions[{#, #2}, #3, Method -> {}] &)]

True

So for OP's example, we can extract the default contour levels using this function as follows:

f[x_, y_] := Log[(1 - x)^2 + 100 (y - x^2)^2];
cp1 = ContourPlot[f[x, y], {x, -2, .5}, {y, -2, 2}];
levels = Cases[cp1, Tooltip[_, t_]:>t, Infinity];

cp2 = ContourPlot[f[x, y], {x, -2, .5}, {y, -2, 2}, 
   Contours -> ((contours = FindDivisions[{#, #2}, #3, Method -> {}]) &)];

cp1 === cp2

True

levels == Reverse@contours

True