{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# 1 - Getting Started\n\nThe main application for :mod:`scikit-gstat <skgstat>` is variogram analysis and `Kriging <https://en.wikipedia.org/wiki/Kriging>`_.\nThis tutorial will guide you through the most basic functionality of ``scikit-gstat``.\nThere are other tutorials that will explain specific methods or attributes in ``scikit-gstat`` in more detail.\n\n**What you will learn in this tutorial**\n\n    * How to instantiate :class:`Variogram <skgstat.Variogram>` and :class:`OrdinaryKriging <skgstat.OrdinaryKriging>`\n    * How to read a variogram\n    * Perform an interpolation\n    * Most basic plotting\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom pprint import pprint\nimport warnings\nwarnings.filterwarnings(\"ignore\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The module is called ``skgstat`` and by convention imported as ``skg``.\nThe :class:`Variogram <skgstat.Variogram>` and :class:`OrdinaryKriging <skgstat.OrdinaryKriging>` classes are available at top level.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import skgstat as skg"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 1.1 Load data\nSciKit-GStat includes a data submodule, that contains some sample datasets. It also offers some basic random sampling on data sources.\nHere we use the well-known Meuse dataset from the R package ``sp`` (https://cran.r-project.org/web/packages/st/index.html). \nIf not specified different, the loading function will only export the lead measurements from the data source.\n\n**Note:** The data is distributed along with the package sp under a GPL-3 license.\nIf you use the data, cite the original authors, not SciKit-GStat for the dataset.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "src = skg.data.meuse()\nprint(src.get('origin'))\n\ncoords, vals = src.get('sample')\n# make a nice table\npd.DataFrame({'x': coords[:, 0], 'y': coords[:, 1], 'lead': vals.flatten()}).head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Get a first overview of your data by plotting the `x` and `y` coordinates and visually inspect how the `z` spread out.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(1, 1, figsize=(9, 9))\nart = ax.scatter(coords[:, 0], coords[:, 1], s=50, c=vals.flatten(), cmap='plasma')\nplt.colorbar(art)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can already see a lot from here:\n\n* There are a few hot-spots, most of them aligned along the river-bank.\n* Most observations show rather low values\n\nThese are already very important insights.\n\n## 1.2 Build a Variogram\nAs a quick reminder, the variogram relates pair-wise separating distances of `coordinates` and relates them to the *semi-variance* of the corresponding `values` pairs. The default estimator used is the Matheron estimator:\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\nFor more details, please refer to the `User Guide <https://mmaelicke.github.io/scikit-gstat/userguide/variogram.html#experimental-variograms>`_.\n\nThe :class:`Variogram <skgstat.Variogram>` class takes at least two arguments. \nThe :func:`coordinates <skgstat.Variogram.coordinates>` and the :func:`values <skgstat.Variogram.values>` observed at these locations.\nIf you use older versions, <ou should also at least set the ``normalize`` parameter to explicitly, as it changed it's default value in\nversion `0.2.8` to ``False``. This attribute affects only the plotting, not the variogram values.\nAdditionally, the number of bins is set to 15, because we have fairly many observations and the default value of 10 is unnecessarily small.\nThe ``maxlag`` set the maximum distance for the last bin. As we have no other additional information about expected correlation lengths,\nwe can either set nothing and use the full distance matrix, or set it i.e. to ``'median'``, which will restrict the distance matrix to only\nuse half of all combinations.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "V = skg.Variogram(coords, vals.flatten(), maxlag='median', n_lags=15, normalize=False)\nfig = V.plot(show=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The upper subplot show the histogram for the count of point-pairs in each lag class. You can see various things here:\n\n  * As expected, there is a clear spatial dependency, because semi-variance increases with distance (blue dots)\n  * The default `spherical` variogram model is well fitted to the experimental data\n  * The shape of the dependency is **not** captured quite well, but fair enough for this example\n  * The first two bins are not well captured, suggesting either the use of a nugget, or a different model.\n\nThe sill of the variogram should correspond with the field variance. The field is unknown, but we can compare the sill to the *sample* variance:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print('Sample variance: %.2f   Variogram sill: %.2f' % (vals.flatten().var(), V.describe()['sill']))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The ``describe`` method will return the most important parameters as a dictionary. \nAnd we can simply print the variogram ob,ect to the screen, to see all parameters.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pprint(V.describe())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "An even shorter summary is printed by the :class:`Variogram <skgstat.Variogram>` object itself:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(V)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 1.3 Kriging\nThe Kriging class will now use the Variogram from above to estimate the Kriging weights for each grid cell. This is done by solving a linear equation system. For an unobserved location $s_0$, we can use the distances to 5 observation points and build the system like:\n\n\\begin{align}\\begin{pmatrix}\n  \\gamma(s_1, s_1) & \\gamma(s_1, s_2) & \\gamma(s_1, s_3) & \\gamma(s_1, s_4) & \\gamma(s_1, s_5) & 1\\\\\n  \\gamma(s_2, s_1) & \\gamma(s_2, s_2) & \\gamma(s_2, s_3) & \\gamma(s_2, s_4) & \\gamma(s_2, s_5) & 1\\\\\n  \\gamma(s_3, s_1) & \\gamma(s_3, s_2) & \\gamma(s_3, s_3) & \\gamma(s_3, s_4) & \\gamma(s_3, s_5) & 1\\\\\n  \\gamma(s_4, s_1) & \\gamma(s_4, s_2) & \\gamma(s_4, s_3) & \\gamma(s_4, s_4) & \\gamma(s_4, s_5) & 1\\\\\n  \\gamma(s_5, s_1) & \\gamma(s_5, s_2) & \\gamma(s_5, s_3) & \\gamma(s_5, s_4) & \\gamma(s_5, s_5) & 1\\\\\n  1 & 1 & 1 & 1 & 1 & 0 \\\\\n  \\end{pmatrix} *\n  \\begin{bmatrix}\n  \\lambda_1 \\\\\n  \\lambda_2 \\\\\n  \\lambda_3 \\\\\n  \\lambda_4 \\\\\n  \\lambda_5 \\\\\n  \\mu \\\\\n  \\end{bmatrix} =\n  \\begin{pmatrix}\n  \\gamma(s_0, s_1) \\\\\n  \\gamma(s_0, s_2) \\\\\n  \\gamma(s_0, s_3) \\\\\n  \\gamma(s_0, s_4) \\\\\n  \\gamma(s_0, s_5) \\\\\n  1 \\\\\n  \\end{pmatrix}\\end{align}\n\nFor more information, please refer to the `User Guide <https://mmaelicke.github.io/scikit-gstat/userguide/kriging.html#kriging-equation-system>`_.\n\nConsequently, the :class:`OrdinaryKriging <skgstat.OrdinaryKriging>` class needs a :class:`Variogram <skgstat.Variogram>` \nobject as a mandatory attribute. Two very important optional attributes are ``min_points`` and ``max_points```.\nThey will limit the size of the Kriging equation system. As we have 200 observations,\nwe can require at least 5 neighbors within the range. More than 15 will only unnecessarily slow down the computation.\nThe ``mode='exact'`` attribute will advise the class to build and solve the system above for each location.\n\n**Note:** The recommended way for kriging applications is to use the interface to :any:`gstools`.\nThere is an easy-to-use interface via :func:`Variogram.to_gstools <skgstat.Variogram.to_gstools>` \nand :func:`Variogram.to_gs_krige <skgstat.Variogram.to_gs_krige>`.\nThe getting started tutorial will use the builtin kriging class, anyway.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ok = skg.OrdinaryKriging(V, min_points=5, max_points=15, mode='exact')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The :func:`OrdinaryKriging.transform <skgstat.OrdianryKriging.transform>` method will apply the interpolation for passed arrays of coordinates.\nIt requires each dimension as a single 1D array. We can easily build a meshgrid of 100x100 coordinates and pass them to the interpolator.\nTo recieve a 2D result, we can simply reshape the result. The Kriging error will be available as the ``sigma`` attribute of the interpolator.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# build the target grid\nx = coords[:, 0]\ny = coords[:, 1]\nxx, yy = np.mgrid[x.min():x.max():100j, y.min():y.max():100j]\nfield = ok.transform(xx.flatten(), yy.flatten()).reshape(xx.shape)\ns2 = ok.sigma.reshape(xx.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And finally, we can plot the result.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axes = plt.subplots(1, 2, figsize=(16, 8))\n\n# rescale the coordinates to fit the interpolation raster\nx_ = (x - x.min()) / (x.max() - x.min()) * 100\ny_ = (y - y.min()) / (y.max() - y.min()) * 100\n\nart = axes[0].matshow(field.T, origin='lower', cmap='plasma', vmin=vals.min(), vmax=vals.max())\naxes[0].set_title('Interpolation')\naxes[0].plot(x_, y_, '+w')\naxes[0].set_xlim((0, 100))\naxes[0].set_ylim((0, 100))\nplt.colorbar(art, ax=axes[0])\nart = axes[1].matshow(s2.T, origin='lower', cmap='YlGn_r')\naxes[1].set_title('Kriging Error')\nplt.colorbar(art, ax=axes[1])\naxes[1].plot(x_, y_, '+w')\naxes[1].set_xlim((0, 100))\naxes[1].set_ylim((0, 100))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "From the Kriging error map, you can see how the interpolation is very certain close to the observation points,\nbut rather high in areas with only little coverage. You could now define a threshold of maximum accepted Kriging error\nvariance to mask the interpolation result.\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
}