The analysis of spatially correlated binary data has received substantial attention in geo-statistical research but is very challenging due to the intricacy of thedistributional 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 through using a latent Gaussian copula model with alternative hypersphere decomposition of covariance matrix. Correctly specifying the covariance matrix is crucial not only for 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 estimatingequations are proposed to estimate the parameters, and 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.