This work aimed at setting up a multivariate and geostatistical methodology to map viticultural terroirs at the province scale (1:125,000). The methodology is based upon the creation of a GIS storing all the viticultural and oenological legacy data of experimental vineyards (1989–2009), long-term climate data, digital elevation model, soilscapes (land systems) and legacy data of soil profiles. Environmental parameters related to viticulture, selected by an explorative PCA, were: elevation, mean annual temperature, mean soil temperature, annual precipitation, clay, sand and gravel content of soils, soil water availability, redoximorphic features and rooting depth. The variables interpolated using geostatistical methods, were used for a k-means clustering aimed to map the Natural Terroir Units (NTU). Vineyards of the province of Siena was subdivided into 9 NTU. Both the historical DOCG (Chianti Classico, Brunello di Montalcino and Nobile di Montepulciano) and other DOCs were characterized by three or four NTU, whereas the wider Chianti and Chianti Colli senesi DOCG was mainly constituted by seven NTU.