How, in practice, are spatial covariances determined? How, in practice, are spatial covariances determined? If one has a single realisation of some observed field, how can the spatial covariance ever be determined? Unless one has access to the many other realisations associated with that particular observation then its surely impossible to ensemble average over all of them? Is the choice of spatial covariance just an educated guess based upon the assumed underlying phenomenon that has given rise to the single observation?
What selection criteria is used to determine the spatial covariance?  
Expansion: I understand that, as some of the comments have pointed out, there are many covariance functions to choose from. Some of these covariance make the assumption that the field is stationary, i.e. an isotropic field as oppose to anisotropic one. These covariance functions are derived on the hyper-sphere and result in the Mat{'}ern family, which includes the exponential cardinal-sine and Gaussian covariance functions. Many of these functions intrinsically contain the rule that on the whole objects have more in common with their nearest neighbours than they do with their distant neighbours.
However given a single observation of a random field, I don't see how one can possibly select one covariance function over another? There is simply not enough information to make an informed choice. How then are researchers and academics supporting the justification in the choice of covariance which will later be used to model the phenomenon? It appears, at least from a naive physicists point of view, that one makes assumptions about the underlying phenomenon to later built a model to make future predictions regarding the underlying phenomenon. There's nothing wrong with this ad hoc approach however I don't see how the solutions to the model can be taken as anything more than hypothetical.
 A: To expand on my comment, commonly "spatial covariance" is associated with Gaussian Processes, which are typically assumed to be stationary. Furthermore, in practice the spatial covariance function (kernel) is assumed to be non-negative and decay with distance. The standard setup is then to assume that your spatial (sample) domain is much larger than the correlation length, so that under the ergodic hypothesis, ensemble averages are well approximated by spatial averages. In this framework, empirical (auto-)correlation functions can be estimated using e.g. spectral methods. So a parametric fit of some particular model is not always required in practice.
Now, in my work I encounter a lot of geostatistics, where all of these assumptions are commonly taken for granted, even though they blatantly do not apply (simple parametric fits are also common). As this is a pet peeve of mine, I would like to take a moment here to contrast the standard geostatistics approach with a technique from computer vision, to demonstrate how one might relax the stationary assumption.
Consider the common geostatistics paradigm of a zero-mean stationary random field $z[\mathbf{x}]$ with variance $\sigma^2$. Typically geostatistics does not use the covariance function directly, but rather uses the variogram, defined for some lag $\mathbf{h}$ as
$$\gamma[\mathbf{h}]\equiv\langle (z[\mathbf{x}+\mathbf{h}]-z[\mathbf{x}])^2\rangle=2(\sigma^2-\kappa[\mathbf{h}])$$
where $\kappa[\mathbf{h}]$ is the (auto-)covariance function, and the last equality comes from expanding the terms in the average and invoking the stationary assumption.
A quite similar construction is commonly used in computer vision, but is derived somewhat differently. If we assume that $z$ is differentiable, then we can approximate the variogram as
$$\gamma[\mathbf{h}]\equiv\langle (z[\mathbf{x}+\mathbf{h}]-z[\mathbf{x}])^2\rangle\approx\langle (\mathbf{S}[\mathbf{x}]^T\mathbf{h})^2\rangle=\mathbf{h}^T\langle \mathbf{S}\mathbf{S}^T\rangle\mathbf{h}\equiv\mathbf{h}^T\mathbf{T}\mathbf{h}$$
where $\mathbf{S}$ is the gradient of $z$ and $\mathbf{T}$ is known as the structure tensor.
A key difference between the geostatistical variogram and the structure tensor is that the latter is not assumed to be stationary. Rather, $\mathbf{T}$ is used to describe the local texture in an image (or in 3D, e.g. seismic or medical imaging data). A simple Google image search will show numerous examples of non-stationary structure, which is perhaps the norm (i.e. the classical geostatistics assumptions are quite commonly invalid).
EDIT: This structure-tensor demo on "Lena" (from a Matlab toolbox by Gabriel Peyre) gives a nice illustration of the texture information the tensor field captures.
 
