Methods for interpolating and approximating three-dimensional scattered data are presented. These methods consist of several local least squares approximations, followed by a piecewise bicubic Hermite interpolant to gridded data, and then optionally followed by a modified Shepard's method. Error bounds are derived for the interpolation and approximation methods that depend on the maximum distance from the nearest data point. The visual smoothness and the discrete errors for these methods applied to known functional data compare favorably with other methods in the literature. The storage and computational complexities of these methods are linear in the number of data points.