Soil surface texture affects many properties of soil including soil structure, water holding capacity, chemistry, organic carbon dynamic and mechanical properties. Digital soil mapping (DSM) was developed in the 1980s with the advancement of remote and proximal sensing technologies, global positioning systems, geostatistics and geographic information systems. The objective of the project was to develop digital soil surface texture maps using soil legacy and earth observation (EO) data. The project was carried out on an experimental sub-watershed (470 ha) located near Quebec City. Particle sizes (sand, silt, and clay) of 167 soil samples were determined by the hydrometer method. Elevation and reflectance data were used as ancillary variables to improve the accuracy of soil surface texture maps. Elevation data was acquired with a RTK GPS. Reflectance data was extracted from Ikonos and Quickbird satellite images acquired on May, 15^th and 30^th, 2008, respectively. Anisotropic and isotropic semivariograms, ordinary kriging and cokriging were computed using Gesotatistical Analyst toolbox of ArGIS v9.3. No significant anisotropy was detected. Geostatistical analyses revealed the presence, within the sub-watershed, of a high spatial structure (C/C_o+C > 0.9) for the sand, silt, and clay contents then allowing DSM. The cokriging approach using ancillary variables reduces the root-mean-square error (RMSE) and the average standard error (ASE). Cokriging of analytical data using elevation and reflectance data as co-variables improves the prediction accuracy of soil texture maps.