You can get a pretty good surface representation using Interpolation.
First combining the data into a single list:
data = Join[contour0, contour25, contour50, contour75, contour100];
Then formatting it into something Interpolation needs for higher dimensions ({{x,y},z})
formatdata = Table[{{data[[i, 3]], data[[i, 1]]}, data[[i, 2]]}, {i, 1, Length@data}];
I should note here that given the form of your data a bit of rearranging is necessary to avoid issues related to non-uniqueness. Now constructing the Interpolation function:
fun = Interpolation[formatdata, InterpolationOrder -> 1];
Again accounting for the reordering of data required for Interpolation the resulting surface can be visualized with
Show[Plot3D[fun[z, y], {z, Min@data[[;; , 3]], Max@data[[;; , 3]]},
{y, Min@data[[;; , 2]], Max@data[[;; , 2]]}],
ListPointPlot3D[data[[;; , {3, 1, 2}]], PlotStyle -> PointSize[Large]]]
Edit:
If you are really concerned with the shape of the surface of sparsely populated regions, like the contour @ 75, create more data points along each contour to drive the behavior.
contour0fun = Interpolation[contour0[[;; , {1, 2}]], InterpolationOrder -> 1];
contour25fun = Interpolation[contour25[[;; , {1, 2}]], InterpolationOrder -> 1];
contour50fun = Interpolation[contour50[[;; , {1, 2}]], InterpolationOrder -> 1];
contour75fun = Interpolation[contour75[[;; , {1, 2}]], InterpolationOrder -> 1];
contour100fun = Interpolation[contour100[[;; , {1, 2}]], InterpolationOrder -> 1];
contour0 = Table[{{0, x}, contour0fun[x]}, {x, 0, 100, .1}];
contour25 = Table[{{25, x}, contour25fun[x]}, {x, 0, 100, .1}];
contour50 = Table[{{50, x}, contour50fun[x]}, {x, 0, 100, .1}];
contour75 = Table[{{75, x}, contour75fun[x]}, {x, 0, 100, .1}];
contour100 = Table[{{100, x}, contour100fun[x]}, {x, 0, 100, .1}];
formatdata = Join[contour0, contour25, contour50, contour75, contour100];
fun = Interpolation[formatdata, InterpolationOrder -> 1];
Show[Plot3D[fun[z, y], {z, Min@data[[;; , 3]], Max@data[[;; , 3]]},
{y, Min@data[[;; , 2]], Max@data[[;; , 2]]}],
ListPointPlot3D[data[[;; , {3, 1, 2}]], PlotStyle -> PointSize[Large]]]
EDIT 2:
Finding the shortest path between two points on a surface
We have now gone from a visualization problem to one of differential geometry. I should say right off the bat that I am not super proficent in this area, but I gave it my best shot. I think you are looking for something called a Geodesic. Having known nothing about this topic prior to your question I went straight to the Google to do some digging. There has been quite a lot of work in the area, and great summary can be found here. To the best of my understanding this is still very much an open problem in mathematics, particularly for mesh based surfaces.
Now, on to this this particular problem. In this case we have a $\left( u, v \right)$ parameterized surface defined as $\left<u,v,f\left(u,v\right)\right>$. With this type of form, the geodestics of the surface may be found through minimizing the arc length functional. For our particular surface definition, with $z=f\left(u,v\right)$, there exists a system of differential equations that can be solve to determine a geodestic (the full form of one of these equations can be found here).
I started by looking an example problem: What is the shortest distance between $\left(0,0\right)$ and $\left(10,10\right)$ on $\left<u,v,u\left(v-1\right)^2\right>$?
We can start by looking at the surface:
Plot3D[u*(v - 1)^2, {u, 0, 10}, {v, 0, 10}, ImageSize -> 400, AxesLabel -> Automatic]
Now we need to setup up our geodesic calucation. Practically, this means setting up a system of coupled differential equations. For this particular example our system looks
sol2=NDSolve[{(x^\[Prime]\[Prime])[t]==-(((-1+y[t])^2 (4 (-1+y[t])
(x^\[Prime])[t] (y^\[Prime])[t]+2 x[t] (y^\[Prime])[t]^2))/(1+4 x[t]^2
(-1+y[t])^2+(-1+y[t])^4)),(y^\[Prime]\[Prime])[t]==-((2 x[t] (-1+y[t]) (4
(-1+y[t]) (x^\[Prime])[t] (y^\[Prime])[t]+2 x[t] (y^\[Prime])[t]^2))/
(1+4 x[t]^2 (-1+y[t])^2+(-1+y[t])^4)),x[0]==0,y[0]==0,x[1]==10,y[1]==10},
{x[t],y[t]},t]
Note that in the system we are looking for parametric equations x[t] and y[t] with $t \in \left[0,1\right]$. Solving this system gives the shortest path between our two points. We can extract an approximate representation of the path with
points2=Table[{sol2[[1,1,2]],sol2[[1,2,2]],sol2[[1,1,2]]*(sol2[[1,2,2]]-1)^2},
{t,0,1,.001}];
and visualize the path on the surface with
Show[Plot3D[x*(y-1)^2,{x,0,10},{y,0,10},ImageSize->400],
ListPointPlot3D[points2]/.Point[a___]:>{Thick,Line[a]}]
You will notice, in general, the shortest path between two points on a surface is not necessarily a straight line, like in this example. Now, moving onto the question of the shortest distance between $\left\{0,0,0\right\}$ and $\left\{100,100,100\right\}$ on the surface defined by the 5 contours we need to reformulate our geodesic equations such that the surface descriptions are consistent. I ran into some issues in NDSolve with respect to extrapolation on the surface (there may be a work around with an application of WhenEvent, but I didn't look too closely into it), so I expanded the domain. In particular, we have
contour0fun=Interpolation[contour0[[;;,{1,2}]],InterpolationOrder->1];
contour25fun=Interpolation[contour25[[;;,{1,2}]],InterpolationOrder->1];
contour50fun=Interpolation[contour50[[;;,{1,2}]],InterpolationOrder->1];
contour75fun=Interpolation[contour75[[;;,{1,2}]],InterpolationOrder->1];
contour100fun=Interpolation[contour100[[;;,{1,2}]],InterpolationOrder->1];
contour10=Table[{{-10,x},contour0fun[x]},{x,0,100,.01}];
contour0=Table[{{0,x},contour0fun[x]},{x,0,100,.01}];
contour25=Table[{{25,x},contour25fun[x]},{x,0,100,.01}];
contour50=Table[{{50,x},contour50fun[x]},{x,0,100,.01}];
contour75=Table[{{75,x},contour75fun[x]},{x,0,100,.01}];
contour100=Table[{{100,x},contour100fun[x]},{x,0,100,.01}];
contour110=Table[{{110,x},contour100fun[x]},{x,0,100,.01}];
bound1=Table[{{i,-10},0},{i,-10.1,110,.0033}];
bound2=Table[{{i,110},100},{i,-10.1,110.1,.0033}];
bound3=Table[{{i,0},0},{i,-10.1,110.1,.0033}];
bound4=Table[{{i,100},100},{i,-10.1,110.1,.0033}];
Combining the data and creating a new expanded surface
formatdata2 = Join[contour10, contour0, contour25, contour50, contour75,
contour100, contour110, bound1, bound2, bound3, bound4];
fun2 = Interpolation[formatdata2, InterpolationOrder -> 1];
Because our surface is an interpolation object we have to take care in constructing the derivatives in the geodesic equations found in NDSolve. Making each derivative its own function
dx[x1_, y1_] := Evaluate[D[fun2[x, y], x]] /. x -> x1 /. y -> y1
dy[x1_, y1_] := Evaluate[D[fun2[x, y], y]] /. x -> x1 /. y -> y1
dxdx[x1_, y1_] := Evaluate[D[fun2[x, y], {x, 2}]] /. x -> x1 /. y -> y1
dydy[x1_, y1_] := Evaluate[D[fun2[x, y], {y, 2}]] /. x -> x1 /. y -> y1
dxdy[x1_, y1_] := Evaluate[D[D[fun2[x, y], x], y]] /. x -> x1 /. y -> y1
Now we can solve our system of equations directly using numeric derivatives of our surface
shortpath=NDSolve[{x''[t]==-(dx[x[t],y[t]]*(y'[t]^2*dydy[x[t],y[t]]
+2*x'[t]*y'[t]*dxdy[x[t],y[t]]+x'[t]^2*dxdx[x[t],y[t]]))
/(1+dx[x[t],y[t]]^2+dy[x[t],y[t]]^2),
y''[t]==-(dy[x[t],y[t]]*(x'[t]^2*dxdx[x[t],y[t]]
+2*y'[t]*x'[t]*dxdy[x[t],y[t]]+y'[t]^2*dydy[x[t],y[t]]))
/(1+dy[x[t],y[t]]^2+dx[x[t],y[t]]^2),x[0]==0,y[0]==0,x[1]==100,y[1]==100},
{x[t],y[t]},{t,0,1}];
We reconstuct our path with
shortpathpoints=Table[{shortpath[[1,1,2]],shortpath[[1,2,2]],
fun2[shortpath[[1,1,2]],shortpath[[1,2,2]]]},{t,0,1,.0001}];
It is possible to calculate the distance of the path by summing the individual line lengths.
Total@(EuclideanDistance @@@Partition[shortpathpoints,2,1])
(*192.702*)
Note that this is exactly the same distance as the straight line path between the two points
Total@(EuclideanDistance @@@Partition[Table[{y,y,fun2[y,y]},{y,0,100,.01}],2,1])
(*192.702*)
Visualizing
Show[Plot3D[fun2[z,y],{z,0,100},{y,0,100}],
ListPointPlot3D[shortpathpoints]/.Point[a___]:>{Thick,Line[a]},
ListPointPlot3D[Table[{y,y,fun2[y,y]},{y,0,100,.01}]]/.Point[a___]:>{Thick,Red,Line[a]}]
This shortest path is entirely dependent on the description of the surface though. If the shapes of the contours or the interpolation between then change, then the shortest path is also likely to change.
Polygons over them. $\endgroup$