The analysis of spatially correlated binary data has received substantial attention in geo-statistical research but is very challenging due to the intricacy of the distributional form. Two principal objectives include examining the dependence of binary response on covariates of interest and quantifying the covariances or correlations between pairs of outcomes. While the literature has sufficiently addressed the modelling issue of the mean structure of a binary response, the characterization of the covariances between pairs of binary responses in terms of covariates is not clear. In this paper, we propose methods to explain such characterizations by using a latent Gaussian copula model with alternative hypersphere decomposition of the covariance matrix. Correctly specifying the covariance matrix is crucial not only for the high efficiency of mean parameters but also for scientific interest. The key is to model the marginal mean and pairwise covariance, simultaneously, for spatial binary data. Two generalized estimating equations are proposed to estimate the parameters, and the asymptotic properties of the resulting estimators are investigated. To evaluate the performance of the methods, we conduct simulation studies and provide real data analysis for illustration.