In this article we address two important issues common to the analysis of large spatial datasets. One is the modeling of nonstationarity, and the other is the computational challenges in doing likelihood-based estimation and kriging prediction. We model the spatial process as a convolution of independent Gaussian processes, with the spatially varying kernel function given by the modified Bessel functions. This is a generalization of the process-convolution approach of Higdon, Swall, and Kern (1999), who used the Gaussian kernel to obtain a closed-form nonstationary covariance function. Our model can produce processes with richer local behavior similar to the processes with the Matérn class of covariance functions. Because the covariance function of our model does not have a closed-form expression, direct estimation and spatial prediction using kriging is infeasible for large datasets. Efficient algorithms for parameter estimation and spatial prediction are proposed and implemented. We compare our method with methods based on stationary model and moving window kriging. Simulation results and application to a rainfall dataset show that our method has better prediction performance. Supplemental materials for the article are available online.