Geostatistics and geographical information system (GIS) procedures are novel techniques helpful for the identification of environmental correlates sustaining contact zones among subspecies or closely related species. In this paper, we tried to infer evolutionary scenarios for Vipera ammodytes across the European part of its distribution area using geostatistics and ecological niche-based models, hence trying to solve several biogeographical questions that remained unclear after the application of classical morphological tools and genetic analyses. Eleven morphological traits from 871 vipers were analysed with geostatistics and ecological niche-based modelling. Interpolation by kriging was used to generate surfaces of morphological variation, which were combined with spatial principal components analysis (SPCA). SPCA maps were used to test putative morphological differentiated groups with discriminant function analysis (DFA). Maximum entropy modelling and seven environmental variables were used to identify factors limiting the distribution of groups and areas for the potential occurrence of such groups. Three patterns of morphological variation were observed: a north-west/south-west cline, transition zones with steep clines of variation in a west–east arc, and particular character traits that disturbed the general cline. SPCA identified between three and nine putative population groups, of which three were supported by DFA. Areas of potential occurrence of these groups were coherent with the range of the three subspecies of V. ammodytes currently recognized. The distribution of all subspecies was mostly related to precipitation in the driest month. Areas of probable sympatry between subspecies are generally small and restricted. The main patterns of geographic variation of morphological characters for V. ammodytes were similar to the patterns obtained for Vipera latastei and Vipera monticola; the same environmental factors limit the distribution of differentiated groups of vipers in the Balkans and the Iberian Peninsula. The influence of humidity on the variation of morphological traits in spatially separated viper taxa from the two European peninsulas coincides with their phylogenetic relatedness. Geostatistics and GIS procedures were successful in the identification of environmental correlates sustaining contact zones among V. ammodytes subspecies in the Balkans. The same techniques should be applied for studying other parapatric forms and refugia regions. © 2010 The Linnean Society of London, Biological Journal of the Linnean Society, 2010, 101, 651–666.