Geostatistics is a science. It is also an art.
Geostatistics is not simply the application of statistical methods to geology-driven spatial distributions, it also provides a conceptual framework for making inferences from Earth sciences data – data which are, more often than not, incomplete.
Some may say, as I would, that most of the problems addressed in geostatistics are inverse problems, in the sense that data are used to infer parameters of the prior model. However, a gap exists between inverse problems and geostatistical problems: in inverse problems modeling the observations can be computer intensive, but the a priori model is typically simple (or simplistic); in a geostatistical problem the data are directly related to the model parameters, this allows one to consider prior models that encapsulate properties of the real Earth, sometimes with breathtaking realism. The gap is narrowing and will disappear in the near future. For the time being, we should try to do the best we can in each of the two fields.
Most geostatistical solutions involve random functions, and a long way has been trod since the simple Gaussian models used in the beginning. The Stanford team has been at the origin of many non-Gaussian developments that have now become standard. They are at it again with the multiple-point geostatistical concept and related algorithms, which allow defining truly complex random functions.