The method to calculate the safety factor of slopes was implemented into a three-dimensional (3D) limit equilibrium model in the open source geographical information system (GIS) GRASS GIS. The Hovland method was exploited to evaluate the slope stability of specified or randomly selected ellipsoidal slip surfaces. The spatial plane projection theory was applied to calculate cell areas of the sliding surface, which increased the computational efficiency. By considering the relationship between the main slope sliding direction and the inclinations of cells, the expression of safety factor was proposed for cells in the sliding surface, which was used to distinguish between positive values and negative ones. Furthermore, the sliding surface was expressed by a triangulated irregular network (TIN) model which was obtained by utilizing constrained Delaunay triangulation method and mesh optimization algorithm. Thus, the shortcoming of the Grid model for accurately representing the sliding border was substantially improved. The safety factor calculated by the TIN sliding surface was closer to the exact solution than the Grid model at the same resolution condition, which demonstrated that the computational precision was obviously enhanced. Finally, the proposed method was verified by the modules programmed in Python language, which was further applied to some typical case studies. It was shown that the safety factor of cells lower than 1.0 represented the unstable region and that greater than 1.0 indicated the resistance slide area in the slip surface. The results can provide an essential basis for the development of landslide treatment measures.