This work investigates the computation of maximum likelihood estimators in Gaussian copula models for geostatistical count data. This is a computationally challenging task because the likelihood function is only expressible as a high dimensional multivariate normal integral. Two previously proposed Monte Carlo methods are reviewed, the Genz-Bretz and Geweke-Hajivassiliou-Keane simulators, and a new method is investigated. The new method is based on the so-called data cloning algorithm, which uses Markov chain Monte Carlo algorithms to approximate maximum likelihood estimators and their (asymptotic) variances in models with computationally challenging likelihoods. A simulation study is carried out to compare the statistical and computational efficiencies of the three methods. It is found that the three methods have similar statistical properties, but the Geweke-Hajivassiliou-Keane simulator requires the least computational effort. Hence, this is the method we recommend. A data analysis of Lansing Woods tree counts is used to illustrate the methods.