A new estimation method for spatial binary probit models is presented. Both spatial auto-regressive (SAR) and spatial error (SEM) models are considered. The proposed estimator relies on the approximation of the likelihood function, that follows a multivariate normal distribution which parameters depend on the spatial structure of the observations. The approximation is inspired by the univariate conditioning procedure proposed by Mendell and Elston, with some modifications to improve accuracy and speed. Very accurate parameter estimations have been achieved in reasonable time for simulated data samples with as much as one million observations. The lessons learned in the Monte Carlo experiment have been applied to a case study on urban sprawl over more than forty thousands plots in Southern France.