A geostatistically based approach with a local regression method is used to predict the magnitude of seasonal streamflow using ocean-atmospheric signals and the hydrological condition of a basin as predictors. The model characterizes the stochastic behavior of a forecast variable by generating a conditional distribution of the predicted value for different hydroclimatic conditions. The correlation structure between dependent and independent variables is represented by the variography of the predicted values in which the distance variable in the variogram is determined by measuring the distance between the predictors. This variogram in a virtual field constructed from the predictors makes it possible to predict variables as unmeasured points while considering historic information as measurement points of the field. Different types of kriging, as well as a generalized linear model regression, are used to predict data in interpolation and extrapolation modes. The forecast skill is evaluated using a linear error in probability space score for different combinations of predictors and different kriging methods. The method is applied to a case study of the Zayandeh-rud River in Isfahan, Iran. The utility of the method is demonstrated for forecasting autumn-winter and spring streamflow using the Southern Oscillation Index, the North Atlantic Oscillation, serial correlation between seasonal streamflow series, and the snow budget. The study analyzes the application of the proposed method in comparison with a K-nearest neighbor regression method. The results of this study show that the proposed method can significantly improve the long-lead probabilistic forecast skill for a nonlinear relationship between hydroclimatic predictors and streamflow in a region.