{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# 4. Plotting\nAt the core of SciKit-GStat is a set of classes, that can be used interactively\nto perform variogram analysis. One important aspect of this analysis is a rich\ncollection of plotting functions. These are directly available as class methods\nof the :class:`Variogram `, \n:class:`DirectionalVariogram ` and \n:class:`SpaceTimeVariogram ` method.\nWith version ``0.3.3``, SciKit-GStat implements two different plotting backend: \n`matplotlib `_ and `plotly `_.\nGenerally speaking, matplotlib is great for creating publication ready figures\nin a variety of formats, including vector-graphic PDF files. Plotly, on the other\nhand, will translate the figure into their Javascript library and open a webbrowser\nwith an interactive plot. This way you can obtain the same figure either for\nyour publication as PDF, or as a HTML object that can be injected into a\nproject report website.\n\nWith the newly introduced :mod:`skgstat.plotting` backend, you can easily read\nand change the backend with a single convenient function. The default backend\nis matplotlib. Please be aware, that `plotly` is only a soft dependency,\nmeaning you need to take care of the installation yourself, to keep\nSciKit-GStat's dependency list shorter.\n\nThe data used to create the :class:`Variogram ` and \n:class:`DirectionalVariogram ` is from \nM\u00e4licke (2021). Here, pancake dataset is used.\nThe spatio-temporal data is derived from Fersch et al. (2020). From that data\npublication, the wireless sensor network data is used. The originaly published\n15 minutes intervals soil temperature data at 20 cm depth was taken for all 55\nstations and aggregated to mean hourly values. To further decrease the data size,\nonly every 6th data point is used here. Estimating the full data set will take\napprox. 120GB RAM and processing took about 30 minutes. The results for the\nthinned data sample are very comparable.\n\nBoth data samples can either be obtained by the orignial publications,\nor from the SciKit-GStat documentation. Both samples are published under\nCreative Commons BY 4.0 license. Please cite the original publications if you\nuse the data, and **not** SciKit-GStat.\n\n**References**\n\nFersch, Benjamin, et al. \"A dense network of cosmic-ray neutron sensors for soil moisture observation in a pre-alpine headwater catchment in Germany.\" Earth System Science Data Discussions 2020 (2020): 1-35.\n\nM\u00e4licke, M.: SciKit-GStat 1.0: A SciPy flavoured geostatistical variogram estimation toolbox written in Python, Geosci. Model Dev. Discuss. [preprint], https://doi.org/10.5194/gmd-2021-174, in review, 2021.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import skgstat as skg\nfrom skgstat.plotting import backend\nimport numpy as np\nimport json\nimport warnings\nimport matplotlib.pyplot as plt\nfrom plotly.subplots import make_subplots\nwarnings.filterwarnings('ignore')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4.1 Load data\nLoad a pancake sample from the data directory.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "c, v = skg.data.pancake(N=300, seed=42).get('sample')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load a artificial random field, generated from a Gaussian covariance function,\nwith a 2x larger range in x-axis direction:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ac, av = skg.data.aniso(N=300, seed=42).get('sample')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load the TERENO soil temperature data from Fersch et al. (2020):\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "with open('./data/tereno_fendt/tereno.json', 'r') as js:\n data_obj = json.load(js)\n\ncoords = np.array(data_obj['coordinates'])\nvals = np.array(data_obj['values'])\nprint(data_obj['description'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Estimate a variogram, with a few more lag classes, as there are enough observation points available.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "V = skg.Variogram(c,v, n_lags=25)\nprint(V)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Estimate the directional variogram with a few more lag classes and an azimuth\nof 90\u00b0. The tolerance is set rather low to illustrate the graphs better\n(fewer point connections.):\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "DV = skg.DirectionalVariogram(ac, av, n_lags=20, azimuth=40., tolerance=15.0)\nprint(DV)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Estimate the spatio-temporal variogram with a product-sum model.\nOnly every 6th hour is taken into account to decrease the memory footprint.\nIf you use the full dataset, you need ^120 GiB RAM. \nThe marginal variograms are kept as they are.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "STV = skg.SpaceTimeVariogram(coords, vals[:,::6], x_lags=20, t_lags=20, model='product-sum')\nprint(STV)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4.2 Backend\n\nYou can switch to `plotly` as a plotting backend by calling the \n:mod:`plotting.backend` function and passing the name of the backend.\nNote that plotly is only a soft dependency and will not automatically be\ninstalled along with SciKit-GStat. You can install it like:\n\n.. code-block:: bash\n\n pip install plotly\n\nNote that in a Jupyter environment you might want to use the plotly.offline\nenvironment to embed the needed Javascript into the notebook. In these cases\nyou have to catch the Figure object and use the iplot function from the\noffline submodule.\n\n## 4.3 Variogram\n\n### 4.3.1 :func:`Variogram.plot `\nThe :func:`Variogram.plot ` is the main plotting\nfunction in SciKit-GStat.\nBefore you use the variogram for further geostatistical methods, like kriging,\nor further analysis, make sure, that a suitable model was found and fitted\nto the experimental data. Further, you have to make sure that the statistical\nfoundation of this estimation is sound, the lag classes are well designed and\nbacked by a suiatable amount of data. \nOtherwise, any other geostatistical analysis or method will have to fail,\nno matter how nice the results might look like.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# from skgstat.plotting import backend\nbackend('plotly')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Plotly\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = V.plot(show=False)\nfig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A useful argument for ``plot`` is the ``ax``, this takes a \n``matplotlib.AxesSubplot`` for the ``'matplotlib'`` backend and a \n``plotly.Figure`` for the ``'plotly'`` backend.\nYou need to supply the correct amount of subplots (two). For convenience,\nthe histogram in the upper subplot can be disabled.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = make_subplots(rows=1, cols=1)\nfig.update_layout(\n width=800,\n height=200,\n template='seaborn',\n showlegend=False, \n margin=dict(l=0, r=0, b=0, t=0)\n)\n\nV.plot(axes=fig, hist=False, show=False)\nfig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The :func:`Variogram.plot ` functions is customizable\nand takes a lot of arguments. However, the same interface is used as for the\n``matplotlib`` version of that function. Many matplotlib arguments are mapped\nto the corresponding plotly arguments. Beyond that, you can either try common\nplotly arguments, or update the figure afterwards:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = V.plot(show=False)\n\nfig.update_layout(\n legend=dict(x=0.05, y=1.1, xanchor='left', yanchor='top', orientation='h'),\n template='plotly_dark',\n annotations=[dict(\n text=\"AWESOME\",\n xref=\"paper\",\n yref=\"paper\",\n x=0.5,\n y=0.5,\n font=dict(color=\"white\", size=100),\n textangle=-30,\n opacity=.3\n )]\n)\nfig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Matplotlib\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "backend('matplotlib')\n\nfig = V.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With matplotlib, you can set any ``matplotlib.AxesSubplot`` as ``axes`` to\nplot on other figures. You can send two axes, for the variogram and the\nhistogram, or only one and disable the histogram plotting.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, ax = plt.subplots(2,2, figsize=(8, 4))\n\nfig = V.plot(axes=ax.flatten()[1], hist=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.3.2 `Variogram.scattergram`\nYou can plot a scattergram of all point pairs formed by the class.\nThe pairs can be grouped by the lag classes, they were formed in. This way you\ncan analyze how the two values of the point pait (head and tail) scatter and\nif this follows a pattern (i.e. anisotropy). It is recommended to use the\n``'plotly'`` backend, as you can click on the legend entries to hide a\nspecific class, or double-click to show only the selected lag class.\nThis makes it much easier to inspect the classes.\n\n#### Plotly\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "backend('plotly')\nfig = V.scattergram(show=False)\nfig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is, however possible to re-create the plot that was used up to SciKit-GStat\nversion ``0.3.0`` with only one color. This is still the default for the\n``'matplotlib'`` backend. \n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = V.scattergram(single_color=True, show=False)\nfig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Matplotlib\nbackend('matplotlib')\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = V.scattergram()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.3.3 `Variogram.location_trend`\nAnother useful helper plot is the\n:func:`location_trend `. This will plot the\nobservation values related to their coordinate position, for each coordinate\ndimension separatedly. With the ``'plotly'`` backend, each dimension will appear\nas a coloured group in a single plot. By double-clicking the legend, you can\ninspect each group separately.\n\nThe ``'plotly'`` backend will automatically switch the used plot type from a\nordinary scatter-plot to a WebGL backed scatter-plot, if there are more than\n5000 observations. This will add some startup-overhead for the plot to appear,\nbut the interactivity actions (like pan, zoom) are speed up by magnitudes.\n\n### Plotly\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "backend('plotly')\nfig = V.location_trend(show=False)\nfig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since version ``0.3.5`` the :func:`location_trend `#\nfunction accepts a ``add_trend_line`` parameter, that defaults to ``False``.\nIf set to true, the class will fit linear models to each of the point clouds\nand output a trend line. It will also calculate the R\u00b2, which you can use to\neither accept the input data as trend free or not (a high R\u00b2 indicates a\n**linear** trend and hence you should decline using the input data).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = V.location_trend(add_trend_line=True, show=False)\nfig\n\n# Matplotlib\n# \"\"\"\"\"\"\"\"\"\"\n# \n# There is a difference between the ``'matplotlib'`` and ``'plotly'`` backend in\n# this plotting function. As Plotly utilizes the legend by default to show and\n# hide traces on the plot, the user can conveniently switch between the\n# coordinate dimensions. \n# In Matplotlib, the figures are not interactive by default and therefore\n# SciKit-GStat will create one subplot for each coordinate dimension.\nbackend('matplotlib')\nfig = V.location_trend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.3.4 :func:`distance_difference plot `\n\nThe final utility plot presented here is a scatter-plot that relates all\npairwise-differences in value to the spatial distance of the respective point\npairs. This can already be considered to be a variogram. For convenience, the\nplotting method will mark all upper lag class edges in the plot. This can already\ngive you an idea, if the number of lag classes is chosen wisely, or if you need to\nadjust. To estimate valid, expressive variograms, this is maybe the most important\npreparation step. If your lag classes do not represent your data well, you will\nnever find a useful variogram.\n\n#### Plotly\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "backend('plotly')\nfig = V.distance_difference_plot(show=False)\nfig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You might also consider to adapt the maximum lag distance using this plot, to\nexclude distances that are not well backed by data. Alternatively,\nthe binning method can be changed. Or both\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "Vcopy = V.clone()\nVcopy.bin_func = 'uniform'\n\nfig = Vcopy.distance_difference_plot(show=False)\nfig\n\n\n# Matplotlib\n# \"\"\"\"\"\"\"\"\"\"\nbackend('matplotlib')\nfig = V.distance_difference_plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4.4 Directional Variogram\n\n### 4.4.1 :func:`pair_field `\n\nThe :class:`DirectionalVariogram ` class is\ninheriting from :class:`Variogram `. Therefore all plotting\nmethod shown above are available for directional variograms, as well.\nAdditionally, there is one more plotting method,\n:func:`DirectionalVariogram.pair_field `.\nThis function will plot all coordinate locations and draw a line between all\npoint pairs, that were not masked by the directional mask array and will, thus,\nbe used for variogram estimation. By default, the method will draw all lines for\nall point pairs and you will see nothing on the plot. But there is also the\npossibility to draw these lines only for a subset of the coordinate locations.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Matplotlib\n# \"\"\"\"\"\"\"\"\"\"\nbackend('matplotlib')\nfig = DV.pair_field()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obviously, one can see the ``azimuth`` (40\u00b0) and narrow ``tolerance`` (15\u00b0)\nsettings in the cone-like shapes of the connection lines, but the whole plot\nis not really instructive or helpful like this. \nUsing the ``points`` keyword, you can show the lines only for a given set of\ncoordinate locations. You have to pass a list of coordinate indices. With\n``add_points=True``, the seleceted points will be highlighted in red.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = DV.pair_field(points=[0, 42, 104, 242], add_points=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Plotly\n\n**Note:** It is not recommended to plot the full\n:func:`pair_field ` with all points\nusing plotly. Due to the implementation, that makes the plot really,\nreally slow for rendering.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "backend('plotly')\nfig = DV.pair_field(points=[0,42, 104, 242], add_points=True, show=False)\nfig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4.5 ST Variogram\nThe :class:`SpaceTimeVariogram ` does not\ninherit from the :class:`Variogram `class and thus,\nits plotting methods are not available for space time variograms.\nHowever, the :class:`SpaceTimeVariogram `\nhas two properties, ``SpaceTimeVariogram.XMarginal`` and ``SpaceTimeVariogram.TMarginal``,\nwhich are both instances of :class:`Variogram `\nfor the spatial and temporal marginal variogram. These instances in turn,\nhave all plotting methods available, in addition to the plotting methods of\n:class:`SpaceTimeVariogram ` itself.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# 4.5.1 `plot(kind='scatter') `\n# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n# \n# The scatter plot can be used to inspect the experimental variogram data on a\n# spatial and temporal axis, with the fitted spatio-temporal model fitted to the data.\n# \n# Plotly\n# \"\"\"\"\"\"\nbackend('plotly')\nfig = STV.plot(kind='scatter', show=False)\nfig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The method can also remove the model from the plot. This can be helpful in\ncase the experimental data should be analyzed. Then, the model plane might be disturbing.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = STV.plot(kind='scatter', no_model=True, show=False)\nfig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And finally, the experimental point data can be connected to a surface grid,\nto emphasize an apparent structure more easily in a 3D plot. This can be done by switching to ``kind='surf'``.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = STV.plot(kind='surf', show=False)\nfig\n\n\n# Matplotlib\n# \"\"\"\"\"\"\"\"\"\"\nbackend('matplotlib')\nfig = STV.plot(kind='surf')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.5.2 :func:`contour `\n3D plots are great for data exploration, especially if they are interactive.\nFor publications, 3D plots are not that helpful. Additionally, it can be quite\ntricky sometimes to find a good angle to focus on the main message of a 3D plot. \nHence, there are more plotting modes. They can either be used by\nsetting ``kind='contour'`` or ``kind='contourf'``. Alternatively, these two\nplotting types also have their own method.\nIn both cases, the experimental semi-variance is plotted on a two dimensional\nplane. The spatial dimension can be found on the x-axis and the temporal\ndimension on the y-axis. The semi-variance itself is shown as a contour plot,\nthat can either only plot the lines (``'contour'``) or filled areas for each\ncontour (``'contourf'``).\n\n#### Plotly\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "backend('plotly')\nfig = STV.contour(show=False)\nfig" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = STV.contourf(show=False)\nfig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Matplotlib\n\nThe matplotlib versions of the contour plots are not that sophisticated,\nbut the returned figure can be adjusted to your needs.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "backend('matplotlib')\nfig = STV.plot(kind='contour')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = STV.plot(kind='contourf')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.5.3 :func`marginals `\nA very important step for the estimation of spatio-temporal variogram models,\nis the estimation of marginal models. While the marginal models are\nimplemented as :class:`Variogram ` instances and can be\nchanged and plotted like any other :class:`Variogram ` instance,\nit can come very handy to plot the marginal models side-by-side.\n\nThis can be done with the :func`marginals ` method.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "backend('plotly')\nfig = STV.marginals(show=False)\nfig" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "backend('matplotlib')\nfig = STV.marginals()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Additionally, the separated spatial and temporal models can be plotted into each sub-plot:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = STV.marginals(include_model=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "backend('plotly')\nfig = STV.marginals(include_model=True, show=False)\nfig" ] } ], "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 }