L = 
You can't. That is, you cannot use calculus to do this, at least not as it is, not directly. That contour is a set of piecewise linear segments. Could you solve the problem in some way with calc? Well, yes. You could do it in the case of some simple functions.
You CAN use my distance2curve utility, found on the File Exchange for free download. (In fact, if you allow it to do so, it does use calculus to solve the problem, in a sense. So if the curve is a spline, then it needs to formulate and solve a problem, where it implicitly minimizes the distance to a smooth curve.)
It seems like you want to learn how to formulate the problem in an analytical form, so I'll do it that way. If all you want to do is solve the problem for a specific contour, then download and use my function. I'll give an example for a simple function though, so you can learn how to solve the problem.
syms x y
z(x,y) = x^2 + x*y + y^2;
And that would be just a simple ellipse, in terms of its contours. We can see one such contour, for example at z(x,y) == 1
fimplicit(z - 1)
axis equal
xlabel x
ylabel y
Now can we formulate the problem in terms of calculus, as you seem to want? The curve is z(x,y)-1==0. Ugh. I can think of at least 4 different ways to solve even this simple problem.
The most direct seems to be to use Lagrange multipliers.
syms lambda x0 y0
We want to mimimize the expression
L = (x - x0)^2 + (y - y0)^2 + lambda*(z(x,y)-1)
We can do that by differentiating, setting the derivatives to zero. The result will be a nonlinear system of three equations, in three unknowns, here x,y, and lambda.
I'll pick some specific point, say (x0,y0) = (2,3).
xylamb = solve(gradient(subs(L,[x0,y0],[2,3]),[x,y,lambda]) == 0,[x,y,lambda],returnconditions = true)
vpa([xylamb.x,xylamb.y])
In there, we can ignore the complex solutions. The one that actually mimizes the distance happens to be the 4th one. (The first is probably the point of maximum distance.)
fimplicit(z - 1)
hold on
plot(xylamb.x(1),xylamb.y(1),'rx')
plot([2,xylamb.x(4)],[3,xylamb.y(4)],'-go')
axis equal
hold off
Ok, can you solve this in a simpler way, that does not force you to use Lagrange multipliers? We can certainly use fmincon.
zlevel = 1;
f = @(xy) xy(1)^2 + xy(1)*xy(2) + xy(2)^2 - zlevel;
nonlcon = @(xy) deal([],f(xy));
xy0 = [2 3];
obj = @(xy) sum((xy - xy0).^2);
xystart = [0,0];
[xysol,fval,exitflag] = fmincon(obj,xystart,[],[],[],[],[],[],nonlcon)
And that is the same as we got from using Lagrange multipliers. But you might argue that it implicitly used Lagrange multipliers anyway.
So, can we do it in a simpler way yet? Probably, since the equation of that ellipse simplifies greatly in a parametric form. Or, you could use my distance2curve utility.
However, your question was on a more general function, such as peaks, where the contour would be less easy to define. For a completely general surface, we might start with contour or contourc.
[X,Y] = meshgrid(linspace(-5,5,101));
Z = peaks(X,Y);
C = contour(X,Y,Z,[1,1])
The result here are two distinct contours, both of which have z(x,y)=1. So we would then find the nearest point on each of those distinct curves, then choose the solution with the smaller distance. As an example, I'll pick the point (x,y)=(-2,1).
xy0 = [-2,1];
[xymin1,d1] = distance2curve(C(:,2:60)',xy0)
[xymin2,d2] = distance2curve(C(:,62:end)',xy0)
So this point fell at a distance of 1.7177 units from the lower curve, but only 0.8364 units from the upper curve.
hold on
plot([xy0(1),xymin1(1)],[xy0(2),xymin1(2)],'r-x',[xy0(1),xymin2(1)],[xy0(2),xymin2(2)],'g-s')
axis equal
