  We propose a method for nonstationary covariance function modeling, based on
the spatial deformation method of Sampson and Guttorp [1992], but using a
low-rank, scalable deformation function written as a linear combination of the
tensor product of B-spline basis. This approach addresses two important
weaknesses in current computational aspects. First, it allows one to constrain
estimated 2D deformations to be non-folding (bijective) in 2D. This requirement
of the model has, up to now,been addressed only by arbitrary levels of spatial
smoothing. Second, basis functions with compact support enable the application
to large datasets of spatial monitoring sites of environmental data. An
application to rainfall data in southeastern Brazil illustrates the method
