Although there is an increase in the amount of seismic data acquired with wide-azimuth geometry, it is difficult to achieve regular data distributions in spatial directions owing to limitations imposed by the surface environment and economic factors. To address this issue, interpolation is an economical solution. The current state-of-the-art methods for seismic data interpolation are iterative methods. However, iterative methods tend to incur high computational costs, which restricts their application in cases of large, high-dimensional data sets. Hence, we have developed a two-step noniterative method to interpolate nonstationary seismic data based on streaming prediction filters (SPFs) with varying smoothness in the time-space domain and we extend these filters to two spatial dimensions. Streaming computation, which is the kernel of the method, directly calculates the coefficients of nonstationary SPF in the overdetermined equation with local smoothness constraints. In addition to the traditional streaming prediction-error filter, we adopt a similarity matrix to improve the constraint condition in which the smoothness characteristics of the adjacent filter coefficient change with the varying data. We also design noncausal-in-space filters for interpolation by using several neighboring traces around the target traces to predict the signal; this is performed to obtain more accurate interpolated results than those from the causal-in-space version. Compared to the Fourier projection onto a convex sets interpolation method, our method has the advantages of fast computational speed and nonstationary event reconstruction. Application of our method on synthetic and nonstationary field data finds that it can successfully interpolate high-dimensional data with low computational cost and reasonable accuracy even in the presence of aliased and conflicting events.