We introduce a numerical procedure for the construction of interpolation and quadrature formulae on bounded regions in the plane. The construction is based on the behavior of spectra of certain multiplication operators and leads to nodes which are inside a prescribed region in R2 . The resulting interpolation schemes are numerically stable and the quadrature formulae have positive weights and almost (but not quite) optimal numbers of nodes. The performance of the algorithm is illustrated by several numerical examples.