In this paper (paper 1), we quantitatively show that the heterogeneous pattern of b values (of the Gutenberg-Richter relation) in the Parkfield segment of the San Andreas fault is to a high degree stationary for the past 35 years. This prepares the grounds for paper 2, where we test the hypothesis that our model of spatially varying b values forecasts future seismicity more accurately than the approach in which one assumes a constant b value equal to the average regional value. The method we develop to measure stationarity in the presence of spatial heterogeneity consists of the following steps: (1) Determine the optimal dimensions of the sampling volume by mapping b values with a wide range of radii and selecting the largest radius that gives the most detailed resolution of the b value heterogeneity. Along the selected fault segment, the high data density permits the definition of the dominant dimensions of the seismotectonic fabric, which is about 8--10 km. (2) Map the difference in b value between two periods, selecting numerous possible catalog divisions. (3) Identify significant changes of b values by the Utsu test (Utsu, 1992). Along the studied fault segment of 110 km length, only one patch of radius 5 km showed a significant increase in b, from below average to above, as a function of time. This change in b initiates around 1993 and thus correlates in space and time with a well-documented episode of creep at depth. Using the derived spatial variable b value distributions, we find that the highest probability for earthquakes with magnitude M ≥ 6 is in the Middle Mountain asperity, where the 1966 Parkfield earthquake nucleated and where all M ≥ 4.5 events in the data set occurred. In contrast, if only the regional average b value of 0.92 is used to predict future seismicity, the creeping segment north of Parkfield should produce major earthquakes most frequently, a conclusion that contradicts the observations.