{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# 2 - Semi-variance Estimators\n\nThis tutorial focuses on experimental variograms.\nIt will guide you through the main semi-variance estimators available in\nSciKit-GStat. Additionally, most of the parameters available for building ###\nan experimental variogram will be discussed.\n\n**In this tutorial you will learn:**\n\n * what estimators are available\n * how they differ\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import skgstat as skg\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom scipy.ndimage import zoom\nskg.plotting.backend('matplotlib')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.1 Load data\nThis tutorial will use the same data samples as the models tutorial.\nThese are samples from the pancake image sampled at different sample sizes.\nWe will use only the second one, but you can re-run all the examples with\nany of the other fields.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "s = [30, 80, 300]\ndata1 = skg.data.pancake(N=s[0], seed=42, as_dataframe=True).get('sample')\ndata2 = skg.data.pancake(N=s[1], seed=42, as_dataframe=True).get('sample')\ndata3 = skg.data.pancake(N=s[2], seed=42, as_dataframe=True).get('sample')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def plot_scatter(data, ax):\n art = ax.scatter(data.x, data.y, 50, c=data.v, cmap='plasma')\n plt.colorbar(art, ax=ax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "make a nice overview of the samples\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 3, figsize=(18, 5))\nfor data, ax in zip((data1, data2, data3), axes.flatten()):\n plot_scatter(data, ax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.2 Comparing estimators\ncalculate all variograms, but use only the second one\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "V1 = skg.Variogram(data2[['x', 'y']].values, data2.v.values, normalize=False, n_lags=8, use_nugget=True)\nV2 = skg.Variogram(data2[['x', 'y']].values, data2.v.values, normalize=False, n_lags=15, use_nugget=True)\nV3 = skg.Variogram(data2[['x', 'y']].values, data2.v.values, normalize=False, n_lags=25, use_nugget=True)\n\nvario = V2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The default estimator configured in :class:`Variogram ` \nis the :func:`Math\u00e9ron estimator `\n(Math\u00e9ron, 1963). It is defined like:\n\n\\begin{align}\\gamma (h) = \\frac{1}{2N(h)} * \\sum_{i=1}^{N(h)}(Z(x_i) - Z(x_{i+h}))^2\\end{align}\n\nwhere:\n\n * $h$ is the distance lag\n * $h$ is the number of observation pairs in $h$-lag class\n * $Z(x_i)$ is the observation at the $i$-th location $x$\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "vario.plot(show=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Following the histogram, we should set a :func:`maxlag `.\nThis property accepts a number 0 < maxlag < 1 to set the :func:`maxlag `\nto this ratio of the maximum separating distance. A number > 1 will use this at an absolute limit.\nYou can also pass ``'mean'`` or ``'median'``.\nThis will calculate and set the mean or median of all distances in the\ndistance matrix as :func:`maxlag `.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "vario.maxlag = 'median'\nvario.plot(show=False);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.3 Alternative estimators\nSciKit-GStat implements more than only the Math\u00e9ron estimator.\nSetting :func:`estimator='cressie' `\nwill set the Cressie-Hawkins estimator.\nIt is implemented as follows (Cressie and Hawkins, 1980):\n\n\\begin{align}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)}}\\end{align}\n\nBy setting :func:`estimator='dowd' `,\nthe Dowd estimator (Dowd, 1984) will be used:\n\n.. matho::\n\n 2\\gamma (h) = 2.198 * {median(Z(x_i) - Z(x_{i+h}))}^2 \n\nFinally, :func:`estimator='genton' `\nwill set the Genton estimator (Genton, 1998):\n\n\\begin{align}\\gamma (h) = 2.2191\\{|V_i(h) - V_j(h)|; i < j\\}_{(k)}\\end{align}\n\nwith: \n\n\\begin{align}k = \\binom{[N_h / 2] + 1}{2}\\end{align}\n\nand:\n\n\\begin{align}q = \\binom{N_h}{2}\\end{align}\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, _a = plt.subplots(1, 3, figsize=(12, 3), sharey=True)\naxes = _a.flatten()\nfor ax, estimator_name in zip(axes, ('matheron', 'cressie', 'dowd')):\n vario.estimator = estimator_name\n vario.plot(axes=ax, hist=False, show=False)\n ax.set_title(estimator_name.capitalize())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The important part is here that the effective range as well as the sill is \nchangeing for the estimator. This will likely change the Kriging result.\nFor Kriging, the difference on the first few lag classes is important,\nas no points will be used for estimation, that lies outside the range.\nWe will zoom in, to actually use a higher resolution. Thus the results\nwill be quite smooth.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "xx, yy = np.mgrid[100:199:100j, 0:99:100j]\nfig, _a = plt.subplots(1, 3, figsize=(18, 6))\naxes = _a.flatten()\n\nfields = []\nfor ax, est in zip(axes, ('matheron', 'cressie', 'dowd')):\n vario.estimator = est\n ok = skg.OrdinaryKriging(vario, min_points=5, max_points=15, mode='exact')\n field = ok.transform(xx.flatten(), yy.flatten()).reshape(xx.shape)\n fields.append(field)\n art = ax.matshow(field, origin='lower', cmap='plasma', vmin=vario.values.min(), vmax=vario.values.max())\n plt.colorbar(art, ax=ax)\n ax.set_title(est.capitalize())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's hard to see this by eye, but there are actually some differences.\nLet's load the original field as well and see which estimator performed best.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pan = skg.data.pancake_field().get('sample')\n\n# rescale to 100x100\n# pan = zoom(pan, 0.2)\n# select the same extend\npan = pan[100:200, 0:100]\n\nfig, axes = plt.subplots(2, 2, figsize=(10,10), sharex=True, sharey=True)\n\n# original\nm = axes[0][0].matshow(pan, origin='lower', cmap='Reds_r')\nplt.colorbar(m, ax=axes[0][0])\naxes[0][0].set_title('Original')\n\n# matheron\nm = axes[0][1].matshow(np.abs(pan - fields[0]), origin='lower', cmap='YlGn_r', vmin=0, vmax=25)\nplt.colorbar(m, ax=axes[0][1])\naxes[0][1].set_title(f'|Original - Matheron|; Diff: {np.abs(pan - fields[0]).mean().round(1)}')\n\n# cressie\nm = axes[1][0].matshow(np.abs(pan - fields[1]), origin='lower', cmap='YlGn_r', vmin=0, vmax=25)\nplt.colorbar(m, ax=axes[1][0])\naxes[1][0].set_title(f'|Original - Cressie| Diff: {np.abs(pan - fields[1]).mean().round(1)}')\n\n# dowd\nm = axes[1][1].matshow(np.abs(pan - fields[2]), origin='lower', cmap='YlGn_r', vmin=0, vmax=25)\nplt.colorbar(m, ax=axes[1][1])\naxes[1][1].set_title(f'|Original - Dowd| Diff: {np.abs(pan - fields[2]).mean().round(1)}')\n\n\nplt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the extracted part of the field, quite some differences are visible.\nThe Matheron estimator is less robust and can't capture the higher values in\nthe upper left corner not as quite well as the other estimators. One can also\nsee, that a substantial amount of the deviations are caused by the noisy\ncharacter of the original image. Note that we loaded the field without\napplying any kind of filter to it. \n\n## 2.4 References\nCressie, N., and D. Hawkins (1980): Robust estimation of the variogram. Math. Geol., 12, 115-125.\n\nDowd, P. A., (1984): The variogram and kriging: Robust and resistant estimators, in Geostatistics for Natural Resources Characterization. Edited by G. Verly et al., pp. 91 - 106, D. Reidel, Dordrecht.\n\nGenton, M. G., (1998): Highly robust variogram estimation, Math. Geol., 30, 213 - 221.\n\nMatheron, G. (1963). Principles of geostatistics. Economic Geology, 58(8), 1246\u20131266. https://doi.org/10.2113/gsecongeo.58.8.1246\n\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" } }, "nbformat": 4, "nbformat_minor": 0 }