Geostatistics and Geographic Information Systems (GIS) are basic tools used for the georeferential analysis of spatial information, and investigation of spatial variability at various scales (Chang et al., 1999). Geostatistics introduce new tools to many scientific disciplines, permitting interpretation of data in space, and obtaining information necessary in the decision making process. Geostatistical modelling involves the application of numerical methods for particular features of spatial attributes for the purpose of generating the probability model (Olea, 2009). Geostatistical methods are optimal estimation methods, if two conditions are met, i.e., data show normal distribution, and are stationary. Considerable deviations from normality or stationarity may cause problems with interpretation of results. Due to this, the first stage in the analysis of spatial data should be the preparation of a data histogram (Bohling, 2005). The histogram is the most popular way of presentation of the empirical distribution of a parameter (Fig. 1). It is primarily applied when analysing high amounts of data differing to a low degree. If a given value is a total or average value of many inconsiderable random factors, irrespective of the distribution of each of the factors, its distribution will be approximate to normal. The term of stationarity of a variable means that its values do not change with time, and the mean and variance of a given variable should not differ considerably in space. The assumption of stationarity of variables in a model is necessary in the case of introduction of distributions of typical test statistics applied in testing hypotheses. Results of many studies show that when a model involves non-stationary variables, the asymptotic distributions of test statistics are non-standard. This may lead to inaccurate results of statistical inference. During research on air quality and degree of soil pollution, the variability of natural conditions and anthropogenic activity may cause disturbances in measurements, leading to weak stationarity. This should be considered at the further stages of geostatistical modelling among others in the preparation of a semivariogram and interpolation of spatial data by means of the kriging or cokriging method (Brenning, 2001). The determination of the spatial distribution of regionalised variable Z(x) or its moments such as among others semivariance requires obtaining many realisations of z1(x), z2(x), … zn(x). In research on the natural environment, however, only a limited number of realisations of a phenomenon with location x is often available. Such observed distributions are frequently unique, e.g., distribution of soil pollution. In order to obtain the relevant number of realisations of the analysed regionalised variable, it is assumed that the studied phenomenon in a certain area is repeated in space (Zawadzki, 2011).