If we want to interpolate a set of data with a curve, we have several choices. We can use polynomial interpolation, spline interpolation or interpolatory subdivision. Polynomial interpolation gives high degrees for large numbers of data points, and hence large computational costs, and can give wiggles in the resulting interpolant. Spline interpolation is often better, but usually are global methods, or require us to specify derivatives. Interpolatory subdivision are local methods, and often better suited for numerical calculations. If we want to interpolate with a surface, we have almost the same methods, but some extra disadvantages for splines and subdivision. To interpolate a surface with splines or subdivision, we generally use a tensor product approach. If the data is nonuniformly spaced, these tensor product methods might give unexpected and unwanted wiggles in the interpolating surface. In this thesis we propose two new methods for surface interpolation based on interpolatory subdivision. Numerical results suggests that the new methods are smooth, and better behaved than the tensor product methods.