Digital Elevation Models (DEM) are an essential component of topographic analysis with applications ranging across many fields within the Geosciences. There are several concerns regarding the accuracy and reliability of interpolated DEMs. In particular, the choice of spatial interpolation method is one of the most important factors, especially given the emergence of Light Detection and Ranging (LiDAR) technology and aerial topographic surveys that provide new sources of elevation data. There is still a matter of uncertainty regarding interpolation algorithm selection and application, especially with LiDAR point cloud data. The case study is divided into two sections. The first section compares the global accuracy of the interpolation algorithms according to the three most relevant variables, including the LiDAR point cloud density and distribution, DEM spatial resolution, and study area topography. The second section sought to use spatial autocorrelation analysis and machine learning to evaluate and model the relative interpolation error. The goal was to gain a better understanding of the influence of the sampling distribution and environmental variables on interpolation accuracy. The results for the first section show that Kriging interpolation was superior to other methods in terms of global interpolation accuracy. Due to its stability, Ordinary Kriging was selected to model the relative interpolation error. The results for the second section showed a clear relationship between terrain complexity and the residual distribution of the interpolated DEM. Overall, complex terrains such as drainage basins were associated with both a greater predictive capability for machine learning models and a higher residual error.