The objective of our study was to expand the R factor of the RUSLE model, erosivity index by its estimation from more readily available rainfall erosivity indexes and parameters in stations without rainfall intensity data, and to determine the most accurate interpolation method for preparing an erosivity index map. Among different erosivity indexes and parameters based on rainfall amounts, only the modified fournier Index (FI_mod) was highly correlated with EI₃₀ in 20 synoptic stations. A local model was used for estimating EI₃₀ from FI_mod in the other 66 stations without rainfall intensity data. The spatial variability of the calculated EI30 in all of the stations was different at an azimuth of 32º when compared to the other directions. Moreover, the nuggetto- sill ratio of the semivariogram (0.27) confirmed a strong spatial correlation of EI₃₀. The inverse distance weighting (IDW), spline, kriging, and cokriging methods with elevation as a covariable were compared by a cross-validation technique. The root mean square error (RMSE) value of the cokriging method when compared to that of the IDW, kriging, and spline methods in the study area declined by 11%, 3%, and 4%, respectively. The output maps for all of the interpolation methods followed similar decreasing trends from west to east, with the highest erosivity index (1,450 MJ·mm·ha^-1·h^-1·y^-1) found in the west. This pattern corresponds with the pattern of climatic change from subhumid to semiarid.