{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# 3 - Variogram Models\n\nThis tutorial will guide you through the theoretical variogram models available for the :class:`Variogram ` class. \n\n**In this tutorial you will learn:**\n\n * how to choose an appropiate model function\n * how to judge fitting quality\n * about sample size influence\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\nimport warnings\nwarnings.filterwarnings('ignore')\nskg.plotting.backend('matplotlib')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3.1 Load data\nFor this example we will use the pancake dataset. You can use the\n:mod:``skgstat.data`` submodule to directly sample the dataset. This is the\nred-channel of an image of an actual pancake. The intersting thing about this pancake is,\nthat it shows some clear spatial structures in its browning, but of different \nshapes at different scales. This should be reflectable with different samples.\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": "markdown", "metadata": {}, "source": [ "Plotting:\n\n" ] }, { "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)\n\n# run\nfig, 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": [ "## 3.2 Comparing theoretical models\nOne of the features of :mod:`skgstat` is the fact that it is programmed object oriented.\nThat means, we can just instantiate a :class:`Variogram ` object\nand start changing arguments unitl it models spatial dependency in our observations well.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "V1 = skg.Variogram(data1[['x', 'y']].values, data1.v.values, maxlag='median', normalize=False)\nV1.plot(show=False);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the others as well\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "V2 = skg.Variogram(data2[['x', 'y']].values, data2.v.values, maxlag='median', normalize=False)\nV3 = skg.Variogram(data3[['x', 'y']].values, data3.v.values, maxlag='median', normalize=False)\n\nfig, _a = plt.subplots(1, 3, figsize=(12, 3), sharey=True)\naxes = _a.flatten()\n\nx = np.linspace(0, V1.maxlag, 100)\n\n# plot each variogram\nfor i, v in enumerate([V1, V2, V3]):\n axes[i].plot(v.bins, v.experimental, '.b')\n axes[i].plot(x, v.fitted_model(x), '-g')\n axes[i].set_title(f'N = {s[i]}')\n axes[i].set_xlabel('Lag (-)')\n if i == 0:\n axes[0].set_ylabel('semivariance (matheron)')\n axes[i].grid(which='major', axis='x')\nplt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see how the experimental variogram changes dramatically with sample size.\nDepending on the sample size, we can also choose different number of lag classes.\nAs the :class:`Variogram `` is object oriented, we can simply\nupdate the binning function. First we set the number of lags directly, then we derive\nit from the distance matrix distribution. In the code below, we build the plot\nfrom scratch, demonstrating how you can access the empirical data and how it is updated, when new parameters are supplied.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axes = plt.subplots(3, 3, figsize=(12, 9), sharey=True, sharex=True)\n\nx = np.linspace(0, V1.maxlag, 100)\nmanual_lags = (6, 12, 18)\ncol_lab = ['10 lags', 'varying lags', 'Scott rule']\n\n# plot each variogram\nfor i in range(3):\n for j, v in enumerate([V1, V2, V3]):\n # first row - use same settings\n if i == 0:\n v.bin_func = 'even'\n v.n_lags = 10\n # second row - use the manual lags\n if i == 1:\n v.n_lags = manual_lags[j]\n # last row - use scott\n if i == 2:\n v.bin_func = 'scott'\n axes[i][j].set_xlabel('Lag (-)')\n \n # plot\n axes[i][j].plot(v.bins, v.experimental, '.b')\n axes[i][j].plot(x, v.fitted_model(x), '-g')\n axes[i][j].grid(which='major', axis='x')\n \n # label first col\n if j == 0:\n axes[i][j].set_ylabel(col_lab[i])\n # title first row\n if i == 0:\n axes[i][j].set_title(f'N = {s[j]}')\nplt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That actually demonstrates how the selection of the experimental variogram can \nhave huge influence on the base data for fitting. Now consider the center column.\nIn each of the plots, the selection of model is not deterministic.\nYou can argue for at least two different models here, that might actually be supported by the empirical data.\nThe :class:`Variogram ` class has several goodness of fit\nmeasures to help you on assessing the fit. This does not replace a careful\nvisual inspection of the models, but can assist you in making an decision.\nRemember that the Kriging will be influenced by the quality of the spatial model,\nespecially on short distances.\nWe can quickly cycle all available models for the sample size of 80 to see\nif spherical fits best. The histogram plot can be turned off.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# we use the settings from before - scott rule\nV2.bin_func = 'scott'\nfig, _a = plt.subplots(2,3, figsize=(12, 6), sharex=True, sharey=True)\naxes = _a.flatten()\nfor i, model in enumerate(('spherical', 'exponential', 'gaussian', 'matern', 'stable', 'cubic')):\n V2.model = model\n V2.plot(axes=axes[i], hist=False, show=False)\n axes[i].set_title('Model: %s; RMSE: %.2f' % (model, V2.rmse))\n axes[i].set_ylim(0, 2000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is quite important. We find all 6 models to describe the experimental\nvariogram more or less equally well in terms of RMSE. Think of the\nimplications: We basically can use any model we like. \nThis is a problem as i.e. the gaussian and the spherical model describe\nfundamentally different spatial properties. Thus, our model selection\nshould be driven by interpretation of the variogram, and not the difference\nin RMSE of only 0.4%, which might very likely not be significant at all.\n\n\nBut what does this difference look like, when it comes to interpolation?\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def interpolate(V, ax):\n xx, yy = np.mgrid[0:499:100j, 0:499:100j]\n ok = skg.OrdinaryKriging(V, min_points=5, max_points=15, mode='exact')\n field = ok.transform(xx.flatten(), yy.flatten()).reshape(xx.shape)\n art = ax.matshow(field, origin='lower', cmap='plasma', vmin=V.values.min(), vmax=V.values.max())\n ax.set_title('%s model' % V.model.__name__)\n plt.colorbar(art, ax=ax)\n return field" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fields = []\nfig, _a = plt.subplots(2,3, figsize=(12, 10), sharex=True, sharey=True)\naxes = _a.flatten()\nfor i, model in enumerate(('spherical', 'exponential', 'gaussian', 'matern', 'stable', 'cubic')):\n V2.model = model\n fields.append(interpolate(V2, axes[i]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Get some basic statistics about the fields\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pd.DataFrame({'spherical': fields[0].flatten(), 'exponential': fields[1].flatten(), 'gaussian': fields[2].flatten(),\n 'matern': fields[3].flatten(), 'stable': fields[4].flatten(), 'cubic': fields[5].flatten()}).describe()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This should illustrate, how important the selection of model is, even if no observation uncertainties are propagated into the analysis.\n\n 1. Gaussian model is far off, producing estimations far outside the observed value ranges\n 2. All other models seem produce quite comparable mean values\n 3. BUT: the standard deviation is quite different\n 4. The median of the field can vary by more than 3 units, even if we took the Gaussian model out\n\nYou have to remind that we had quite some observations. The selection of model becomes even more arbitrary with smaller samples and more importantly: We have to consider more than one equally probable parameterization of each model when the experimental is more spreaded.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Finally, we can calculate the difference between the kriging fields to inspect the spread of estimations spatially:\n#\nfield_min = np.nanmin(np.stack(fields, axis=2), axis=2)\nfield_max = np.nanmax(np.stack(fields, axis=2), axis=2)\n\nfig, ax = plt.subplots(1, 1, figsize=(7,7))\nm = ax.matshow(field_max - field_min, origin='lower', cmap='Reds')\nplt.colorbar(m)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The colorbar is spanning the entire value range. Thus, given the minor differences in the fitting of the models, we would have to reject just any estimation based on an automatic fit, which is considering some uncertainties in the selection of parameters, because the RMSE values were too close.\n\nTo use the result from above, we need to justfy the selection of model first and manually fit the model based on expert knowledge.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3.3 Using other sample sizes\nLet's have a look at the sparse sample again\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "V1.plot(show=False);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is a nugget-effect variogram. Thus we have to reject any geostatistical \nanalysis based on this sample. It just does not expose any spatial pattern that \ncan be exploited.\n\nWhat about the denser sample. Increasing the sample size should reject some \nof the models. Remind, that we are sampling at more short distances and thus,\nthe variogram will be governed by the short ranged patterns of the field, while\nthe other samples are more dependent on the medium and large range patterns, as\nthere were less short separating distances sampled.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# we use the settings from before - scott rule\nV3.bin_func = 'scott'\nfig, _a = plt.subplots(2,3, figsize=(12, 6), sharex=True, sharey=True)\naxes = _a.flatten()\nfor i, model in enumerate(('spherical', 'exponential', 'gaussian', 'matern', 'stable', 'cubic')):\n V3.model = model\n V3.plot(axes=axes[i], hist=False, show=False)\n axes[i].set_title('Model: %s; RMSE: %.2f' % (model, V3.rmse))\n axes[i].set_ylim(0, 2000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now clearly reject the cubic, gaussian and exponential model. \nI personally would also reject the spherical model we used in the fist place,\nas it is systematically underestimating the semi-variance on short distances. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "d_fields = []\nfig, _a = plt.subplots(2,3, figsize=(18, 12), sharex=True, sharey=True)\naxes = _a.flatten()\nfor i, model in enumerate(('spherical', 'exponential', 'gaussian', 'matern', 'stable', 'cubic')):\n V3.model = model\n d_fields.append(interpolate(V3, axes[i]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again some statistics\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pd.DataFrame({'spherical': d_fields[0].flatten(), 'exponential': d_fields[1].flatten(), 'gaussian': d_fields[2].flatten(),\n 'matern': d_fields[3].flatten(), 'stable': d_fields[4].flatten(), 'cubic': d_fields[5].flatten()}).describe()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, if we only concentrate on the not-rejected models: matern and stable,\nwe can see hardly any difference in the field. Additionally, except for extrema,\nthe statistical properties of the two fields are largely the same.\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 }