Subsurface stratigraphy of multi-layered slopes is essential and crucial for slope stability analysis. It is usual practice for engineers to interpret stratigraphic boundaries separating different soil layers using both site investigation data and prior knowledge of local geology, but such practice might encounter significant challenge when the site data are very limited. In addition, uncertainty in stratigraphic boundaries has not been explicitly or quantitatively considered in planning of site investigation (e.g., determination of borehole number and locations). There lacks a quantitative and objective tool to determine the optimal locations and number of boreholes for slope stability analysis while accounting for stratigraphic uncertainty. In this study, a smart sampling strategy based on multiple point statistics and information entropy is proposed for delineation of slope subsurface stratigraphy and planning of geotechnical boreholes. It is a data-driven approach that enables an ensemble of prior knowledge within a training image using multiple point statistics. The proposed method not only provides evolution of the most probable interpolation from sparse measurements and the associated interpolation uncertainties, but also adaptively determines the optimal locations of boreholes. Effectiveness of the proposed method is illustrated and validated through both a simulation example and a real case. It is found that the data-driven framework can automatically identify locations of largest interpolation uncertainty within a multi-layered slope conditional on its outcrops, and that the associated stratigraphic uncertainty gradually reduces as borehole number increases. More importantly, the optimal number and locations of boreholes required for slope stability analysis are adaptively determined by the proposed method.