<jupyter_start><jupyter_text>3 - Semi-variance Estimators This tutorial focuses on experimental variograms. It will guide you through the main semi-variance estimators available in `scikit-gstat`. Additionally, most of the parameters available for building an experimental variogram will be discussed.**In this tutorial you will learn:*** what estimators are available* how they differ<jupyter_code>from skgstat import Variogram, OrdinaryKriging
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')
%env SKG_SUPPRESS=true<jupyter_output>env: SKG_SUPPRESS=true<jupyter_text>3.1 Load data There are three prepared data sets in the `./data` folder. Each of them is a generated random field with different underlying spatial properties. We will use only the thrid one, but you can re-run all the examples with any of the other fields.<jupyter_code>data1 = pd.read_csv('data/sample_matern_15.csv')
data2 = pd.read_csv('data/sample_matern_40.csv')
data3 = pd.read_csv('data/sample_spherical_noise.csv')
def plot_scatter(data, ax):
    art = ax.scatter(data.x, data.y, 50, c=data.z, cmap='plasma')
    plt.colorbar(art, ax=ax)
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
for data, ax in zip((data1, data2, data3), axes.flatten()):
    plot_scatter(data, ax)<jupyter_output><empty_output><jupyter_text>3.2 Comparing estimators The default estimator configured in `Variogram` is the Mathéron estimator (Mathéron, 1963). It is defined like:$$ \gamma (h) = \frac{1}{2N(h)} * \sum_{i=1}^{N(h)}(Z(x_i) - Z(x_{i+h}))^2 $$where:* $h$ is the distance lag* $N(h)$ is the number of observation pairs in $h$-lag class* $Z(x_i)$ is the observation at the $i$-th location $x$<jupyter_code>V1 = Variogram(data3[['x', 'y']].values, data3.z.values, normalize=False, n_lags=20, use_nugget=True)
V1.plot(show=False);<jupyter_output><empty_output><jupyter_text>Following the histogram, we should set a `maxlag`. This property accepts a number $0  1 will use this at an absolute limit. You can also pass `'mean'` or `'median'`. This will calculate and set the mean or median of all distances in the distance matrix as `maxlag`.<jupyter_code>V1.maxlag = 'median'
V1.plot(show=False);<jupyter_output><empty_output><jupyter_text>3.3 Alternative estimators `scikit-gstat` implements more than only the Mathéron estimator. Setting `estimator='cressie'` will set the Cressie-Hawkins estimator. It is implemented as follows (Cressie and Hawkins, 1980):$$ 2\gamma (h) = \frac{\left(\frac{1}{N(h)} \sum_{i=1}^{N(h)} |Z(x_i) - Z(x_{i+h})|^{0.5}\right)^4}        {0.457 + \frac{0.494}{N(h)} + \frac{0.045}{N^2(h)}} $$By setting `estimator='dowd'`, the Dowd estimator (Dowd, 1984) will be used:$$ 2\gamma (h) = 2.198 * {median(Z(x_i) - Z(x_{i+h}))}^2 $$Finally, `estimator='genton'` will set the Genton estimator (Genton, 1998):$$ \gamma (h) = 2.2191\{|V_i(h) - V_j(h)|; i < j\}_{(k)} $$with: $$ k = \binom{[N_h / 2] + 1}{2} $$and:$$ q = \binom{N_h}{2} $$<jupyter_code>fig, _a = plt.subplots(2, 2, figsize=(18, 12))
axes = _a.flatten()
for ax, est in zip(axes, ('matheron', 'cressie', 'dowd', 'genton')):
    V1.estimator = est
    V1.plot(axes=ax, hist=False, show=False)
    ax.set_title(est.capitalize())<jupyter_output><empty_output><jupyter_text>For Kriging, the difference on the first few lag classes is important, as no points will be used for estimation, that lies outside the range.<jupyter_code>xx, yy = np.mgrid[0:99:100j, 0:99:100j]
fig, _a = plt.subplots(1, 3, figsize=(18, 6))
axes = _a.flatten()

for ax, est in zip(axes, ('matheron', 'cressie', 'dowd')):
    V1.estimator = est
    ok = OrdinaryKriging(V1, min_points=5, max_points=15, mode='exact')
    field = ok.transform(xx.flatten(), yy.flatten()).reshape(xx.shape)
    art = ax.matshow(field, origin='lower', cmap='plasma')
    plt.colorbar(art, ax=ax)
    ax.set_title(est.capitalize())<jupyter_output>Warning: for 6 locations, not enough neighbors were found within the range.
Warning: for 27 locations, not enough neighbors were found within the range.
Warning: for 286 locations, not enough neighbors were found within the range.<jupyter_text>You can see from these results that the Cressie and the Dowd estimator pronounce the extreme values more. The original field is also in the `./data` folder. We can load it for comparison.<jupyter_code>rf = np.loadtxt('data/rf_spherical_noise.txt')
fig, ax = plt.subplots(1, 1, figsize=(8,8))
art = ax.matshow(rf, origin='lower', cmap='plasma')
plt.colorbar(art, ax=ax)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14633068081343945, 'fraction_numerical': 0.028293545534924844, 'mean_word_length': 5.529581529581529, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '72394', 'license': 'MIT', 'n_tokens_mistral': 1715, 'n_tokens_neox': 1609, 'n_words': 483, 'path': 'docs/tutorials/03_estimators.ipynb', 'repo_name': 'rhugonnet/scikit-gstat', 'size': '4524'}