{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Automatic fitting\n\nIn order to demonstrate how to automatically fit normalizer and variograms,\nwe generate synthetic log-normal data, that should be interpolated with\nordinary kriging.\n\nNormalizers are fitted by minimizing the likelihood function and variograms\nare fitted by estimating the empirical variogram with automatic binning and\nfitting the theoretical model to it. Thereby the sill is constrained to match\nthe field variance.\n\n## Artificial data\n\nHere we generate log-normal data following a Gaussian covariance model.\nWe will generate the \"original\" field on a 60x60 mesh, from which we will take\nsamples in order to pretend a situation of data-scarcity.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport gstools as gs\n\n# structured field with edge length of 50\nx = y = range(51)\npos = gs.generate_grid([x, y])\nmodel = gs.Gaussian(dim=2, var=1, len_scale=10)\nsrf = gs.SRF(model, seed=20170519, normalizer=gs.normalizer.LogNormal())\n# generate the original field\nsrf(pos)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here, we sample 60 points and set the conditioning points and values.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ids = np.arange(srf.field.size)\nsamples = np.random.RandomState(20210201).choice(ids, size=60, replace=False)\n\n# sample conditioning points from generated field\ncond_pos = pos[:, samples]\ncond_val = srf.field[samples]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Fitting and Interpolation\n\nNow we want to interpolate the \"measured\" samples\nand we want to normalize the given data with the BoxCox transformation.\n\nHere we set up the kriging routine and use a :any:`Stable` model, that should\nbe fitted automatically to the given data\nand we pass the :any:`BoxCox` normalizer in order to gain normality.\n\nThe normalizer will be fitted automatically to the data,\nby setting ``fit_normalizer=True``.\n\nThe covariance/variogram model will be fitted by an automatic workflow\nby setting ``fit_variogram=True``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "krige = gs.krige.Ordinary(\n    model=gs.Stable(dim=2),\n    cond_pos=cond_pos,\n    cond_val=cond_val,\n    normalizer=gs.normalizer.BoxCox(),\n    fit_normalizer=True,\n    fit_variogram=True,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, let's have a look at the fitting results:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(krige.model)\nprint(krige.normalizer)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As we see, it went quite well. Variance is a bit underestimated, but\nlength scale and nugget are good. The shape parameter of the stable model\nis correctly estimated to be close to `2`,\nso we result in a Gaussian like model.\n\nThe BoxCox parameter `lmbda` was estimated to be almost 0, which means,\nthe log-normal distribution was correctly fitted.\n\nNow let's run the kriging interpolation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "krige(pos)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plotting\n\nFinally let's compare the original, sampled and interpolated fields.\nAs we'll see, there is a lot of information in the covariance structure\nof the measurement samples and the field is reconstructed quite accurately.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(1, 3, figsize=[8, 3])\nax[0].imshow(srf.field.reshape(len(x), len(y)).T, origin=\"lower\")\nax[1].scatter(*cond_pos, c=cond_val)\nax[2].imshow(krige.field.reshape(len(x), len(y)).T, origin=\"lower\")\n# titles\nax[0].set_title(\"original field\")\nax[1].set_title(\"sampled field\")\nax[2].set_title(\"interpolated field\")\n# set aspect ratio to equal in all plots\n[ax[i].set_aspect(\"equal\") for i in range(3)]"
      ]
    }
  ],
  "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.8.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}