The focus of this dissertation is the interpolation and approximation of multivariate bandlimited functions via sampled (function) values. The first set of results investigates polynomial interpolation in connection with multivariate bandlimited functions. To this end, the concept of a uniformly invertible Riesz basis is developed (with examples), and is used to construct Lagrangian polynomial interpolants for particular classes of sampled square-summable data. These interpolants are used to derive two asymptotic recovery and approximation formulas. The first recovery formula is theoretically straightforward, with global convergence in the appropriate metrics; however, it becomes computationally complicated in the limit. This complexity is sidestepped in the second recovery formula, at the cost of requiring a more local form of convergence. The second set of results uses oversampling of data to establish a multivariate recovery formula. Under additional restrictions on the sampling sites and the frequency band, this formula demonstrates a certain stability with respect to sampling errors. Computational simplifications of this formula are also given.