Recently Dietrich and Newsam <1993> derived a fast and exact method for generating unconditional realizations of a stationary, multidimensional Gaussian random field on a rectangular sampling grid. The method is based on embedding the random field covariance matrix in a larger positive definite matrix with circulant/block circulant structure. The circulant structure of the embedding matrix means that a square root of this matrix can be efficiently computed by the fast Fourier transform; realizations are then generated by multiplying vectors of white noise by this square root. This paper extends the method to generating realizations conditioned on direct and/or indirect measurements of the field collected at an arbitrary set of scattered data points. ¿ American Geophysical Union 1996