{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# 6 - GSTools\nWith version ``0.5`` ``scikit-gstat`` offers an interface to the awesome `gstools `_\nlibrary. This way, you can use a :class:`Variogram ` estimated with ``scikit-gstat`` \nin `gstools `_ to perform random field generation, kriging and much, much more.\n\nFor a :class:`Variogram ` instance, there are three possibilities to export into `gstools `_ : \n\n 1. :func:`Variogram.get_empirical(bin_center=True) ` returns a pair of distance lag bins and experimental semi-variance values, like `gstools.variogram.vario_estimate `_. \n 2. :func:`Variogram.to_gstools ` returns a parameterized :any:`CovModel ` derived from the Variogram.\n 3. :func:`Variogram.to_gs_krige ` returns a :any:`GSTools Krige ` instance based on the variogram\n\n\n## 6.1 ``get_empirical``\n\n### 6.1.1 Reproducing the gstools example\n\nYou can reproduce the `Getting Started example for variogram estimation from GSTools docs `_ \nwith ``scikit-gstat``, and replace the calculation of the empirical variogram with :class:`skg.Variogram `. \n\nNote: This does only make sense if you want to use a distance metric, binning procedure or semi-variance estimator, that is not included in `gstools` or are bound to `scikit-gstat` for any other reason. :class:`Variogram ` will _always_ perform a full model fitting cycle on instantiation, which could lead to some substantial overhead here.\nThis behavior might change in a future version of `scikit-gstat`.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# import\nimport skgstat as skg\nimport gstools as gs\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport warnings\nwarnings.filterwarnings('ignore')\nskg.plotting.backend('matplotlib')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "use the example from gstools\ngenerate a synthetic field with an exponential model\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x = np.random.RandomState(19970221).rand(1000) * 100.\ny = np.random.RandomState(20011012).rand(1000) * 100.\nmodel = gs.Exponential(dim=2, var=2, len_scale=8)\nsrf = gs.SRF(model, mean=0, seed=19970221)\nfield = srf((x, y))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "combine x and y for use in skgstat\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "coords = np.column_stack((x, y))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the example, :any:`gstools.variogram.vario_estimate` is used to estimate the empirical variogram:\n\n.. code-block:: python\n\n # estimate the variogram of the field\n bin_center, gamma = gs.vario_estimate((x, y), field)\n\n\nHere, we can use :class:`skg.Variogram `. \nFrom the shown arguments, :func:`estimator ` and\n:func:`bin_func ` are using the default values:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "V = skg.Variogram(coords, field, n_lags=21, estimator='matheron', maxlag=45, bin_func='even')\nbin_center, gamma = V.get_empirical(bin_center=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And finally, the exact same code from the GSTools docs can be called:\nfit the variogram with a stable model. (no nugget fitted)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fit_model = gs.Stable(dim=2)\nfit_model.fit_variogram(bin_center, gamma, nugget=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Output the model\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ax = fit_model.plot(x_max=max(bin_center))\nax.scatter(bin_center, gamma)\nprint(fit_model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.1.2 ``bin_center=False``\n\nIt is important to understand, that ``gstools`` and ``skgstat`` are handling lag bins different.\nWhile ``skgstat`` uses the upper limit, ``gstools`` assumes the bin center.\nThis can have implications, if a model is fitted. C\nonsider the example below, in which only the ``bin_center`` setting is different.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "bin_edges, _ = V.get_empirical(bin_center=False)\n\n# fit the variogram with a stable model. (no nugget fitted)\nedge_model = gs.Stable(dim=2)\n_ = edge_model.fit_variogram(bin_edges, gamma, nugget=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make a nice plot\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n\n# plot first\nfit_model.plot(ax=axes[1], label='center=True')\n# plot second\nedge_model.plot(ax=axes[1], label='center=False')\n\n# bins\naxes[0].scatter(bin_center, gamma, label='center=True')\naxes[0].scatter(bin_edges, gamma, label='center=False')\n\naxes[0].set_title('Empirical Variogram')\naxes[1].set_title('Variogram Model')\naxes[0].legend(loc='lower right')\nprint(fit_model)\nprint(edge_model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice the considerable gap between the two model functions. This can already lead to seroius differences, i.e. in Kriging.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.1.3 Using other arguments\n\nNow, with the example `from the GSTools docs `_ working,\nwe can start chaning the arguments to create quite different empirical variograms.\n\n**Note**: This should just illustrate the available possibilities, the result is by no means producing a better\nestimate of the initially created Gaussian random field.\n\nIn this example different things will be changed:\n\n- use only 15 lag classes, but distribute the point pairs equally. Note the differing widths of the classes. (``bin_func='uniform'``)\n- The :func:`Dowd ` estimator is used. (``estimator='dowd'``)\n- The Taxicab metric (https://en.wikipedia.org/wiki/Taxicab_geometry) (aka. Manhattan metric or cityblock metric) is used over\nEuklidean for no obvious reason. (``dist_func='cityblock'``)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "V = skg.Variogram(coords, field, n_lags=15, estimator='dowd', maxlag=45, bin_func='uniform', dist_func='cityblock')\nbin_center, gamma = V.get_empirical(bin_center=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "fit the variogram with a stable model. (no nugget fitted)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fit_model = gs.Stable(dim=2)\nfit_model.fit_variogram(bin_center, gamma, nugget=True)\n\n# output\nax = fit_model.plot(x_max=max(bin_center))\nax.scatter(bin_center, gamma)\nprint(fit_model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you fit the `gs.Stable` with a nugget, it fits quite well. But keep in mind that this does not necessarily describe the original field very well and was just fitted for demonstration.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.2 ``to_gstools``\n\nThe second possible interface to ``gstools`` is the :func:`Variogram.to_gstools ` function.\nThis will return one of the classes `listed in the gstools documentation `_.\nThe variogram parameters are extracted and passed to gstools. You should be able to use it, just like any other :any:`CovModel `.\n\nHowever, there are a few things to consider:\n\n- ``skgstat`` can only export isotropic models.\n- The ``'harmonize'`` cannot be exported\n\n### 6.2.1 exporting :class:`Variogram `\n\nIn this example, the same Variogram from above is estimated, but we use the :func:`exponential ` model. \nAn exponential covariance function was used in the first place to create the field that was sampled.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "skg.plotting.backend('plotly')\nV = skg.Variogram(coords, field, n_lags=21, estimator='matheron', model='exponential', maxlag=45, bin_func='even')\nV.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now export the model to ``gstools``:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "exp_model = V.to_gstools()\nprint(exp_model)\n\n# get the empirical for the plot as well\nbins, gamma = V.get_empirical(bin_center=True)\n\nax = exp_model.plot(x_max=45)\nax.scatter(bins, gamma)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note:** It is important to understand, that ``skgstat`` and ``gstools`` handle coordinates slightly different.\nIf you export the :class:`Variogram ` to a :any:`CovModel `\nand you want to use the :class:`Variogram.coordinates `, you **must** transpose them.\n\n.. code-block:: python\n\n # variogram is a skgstat.Variogram instance\n model = variogram.to_gstools()\n cond_pos = variogram.coordinates.T\n\n # use i.e. in Kriging\n krige = gs.krige.Ordinary(model, cond_pos, variogram.values)\n\n### 6.2.2 Spatial Random Field Generation\n\nWith a :any:`CovModel `, we can use any of the great tools implemented in ``gstools``.\nFirst, let's create another random field with the exponential model that we exported in the last section:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x = y = range(100)\nnew_field = gs.SRF(exp_model, seed=13062018)\nnew_field.structured([x, y])\nnew_field.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Keep in mind, that we did not call a Kriging procedure, but created **another** field.\nOf course, we can do the same thing with the more customized model, created in 6.1.3:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "malformed = gs.SRF(fit_model, seed=24092013)\nmalformed.structured([x, y])\nmalformed.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice how the spatial properties as well as the value range has changed. \nThat's why it is important to estimate :class:`Variogram ` or :any:`CovModel ` \ncarefully and not let the GIS do that for you somewhere hidden in the dark.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.3 ``to_gs_krige``\n\nFinally, after carefully esitmating and fitting a variogram using SciKit-GStat, \nyou can also export it directly into a :any:`GSTools Krige ` instance. \nWe use the variogram as in the other sections:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# export\nkrige = V.to_gs_krige(unbiased=True) # will result in ordinary kriging\nprint(krige)\n\n# create a regular grid\nx = y = range(100)\n\n# interpolate\nresult, sigma = krige.structured((x, y))\n\nfig, axes = plt.subplots(1, 2, figsize=(8, 4))\n\n# plot\naxes[0].imshow(result, origin='lower')\naxes[1].imshow(sigma, origin='lower', cmap='RdYlGn_r')\n\n# label\naxes[0].set_title('Kriging')\naxes[1].set_title('Error Variance')\n\nplt.tight_layout()" ] } ], "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 }