Variograms_Dresden.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "db261183",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47a839b6",
   "metadata": {},
   "source": [
    "## Variograms Example - Ground Elevation Data for Dresden"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "414d668f",
   "metadata": {},
   "source": [
    "### Import Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f3b6cabd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import math\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import least_squares\n",
    "from ipywidgets import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baf8597a",
   "metadata": {},
   "source": [
    "<div align=center>\n",
    "<em> <font size=4> \"Everything is related to everything else, but near things are more related than distant things!\"  </font size></em> \n",
    "<br>\n",
    "<i> Waldo Tobler (1970) </i> </div align>\n",
    "\n",
    "### Variograms\n",
    "In order to evaluate the dependency of a variable's value on its position, all samples are compared with regard to their similarity and proximity. We will use a small data set with the ground height at 35 randomly determined positions in the city of Dresden:\n",
    "### Load Data from `.csv`-file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8a808e64",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(35, 5)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Ort</th>\n",
       "      <th>Nummer</th>\n",
       "      <th>Rechtswert [m]</th>\n",
       "      <th>Hochwert [m]</th>\n",
       "      <th>Hoehe [m ue. NN]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Altmarkt</td>\n",
       "      <td>1</td>\n",
       "      <td>5411525</td>\n",
       "      <td>5656080</td>\n",
       "      <td>114</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Zwinger</td>\n",
       "      <td>2</td>\n",
       "      <td>5411250</td>\n",
       "      <td>5656485</td>\n",
       "      <td>110</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Hbf</td>\n",
       "      <td>3</td>\n",
       "      <td>5411185</td>\n",
       "      <td>5655012</td>\n",
       "      <td>115</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Chemiebau</td>\n",
       "      <td>4</td>\n",
       "      <td>5410975</td>\n",
       "      <td>5653670</td>\n",
       "      <td>143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Wasaplatz</td>\n",
       "      <td>5</td>\n",
       "      <td>5413006</td>\n",
       "      <td>5653677</td>\n",
       "      <td>123</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Ort  Nummer  Rechtswert [m]  Hochwert [m]  Hoehe [m ue. NN]\n",
       "0   Altmarkt       1         5411525       5656080               114\n",
       "1    Zwinger       2         5411250       5656485               110\n",
       "2        Hbf       3         5411185       5655012               115\n",
       "3  Chemiebau       4         5410975       5653670               143\n",
       "4  Wasaplatz       5         5413006       5653677               123"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# load the data from the excel-file, resulting in a pandas.DataFrame\n",
    "#     (a table-like data structure)\n",
    "# full path, e.g., \"C:/User/.../HGHCM/Data_Dresden.csv\" can be used as well\n",
    "df = pd.read_csv(\"Data_Dresden.csv\", sep=\";\", encoding='latin-1')\n",
    "\n",
    "# print shape of data (n_rows, n_cols)\n",
    "print(df.shape)\n",
    "# print first 5 rows of data set\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b20ab07e",
   "metadata": {},
   "source": [
    "### Cumulative Histogram\n",
    "In order to get an overview on the distribution of ground heights all over the city, we will take a look on a cumulative histogram of the values, assuming that the sum curve of the real ground heights would look similar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "558d53ca",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e18e2d69eb164015889d4571ed2ea704",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# calculate mean and median\n",
    "# with DataFrame.iloc, we can access the values in the table indexed by\n",
    "#     [row, col] --> (a colon \":\" can be used to say something like \"from here\n",
    "#     until the end\" or \"from the beginning until here\")\n",
    "mean = np.mean(df.iloc[:,-1])\n",
    "median = np.median(df.iloc[:, -1])\n",
    "\n",
    "# plot histogram\n",
    "ax = df.iloc[:, -1].hist(cumulative=True, density=True, histtype=\"step\", lw=2,\n",
    "                         bins=df.shape[0], label=\"Ground Height\", figsize=(10, 6))\n",
    "# plot mean and median\n",
    "ax.axvline(x=median, ls='--', color='red', label=\"Median\")\n",
    "ax.axvline(x=mean, ls='--', color='green', label=\"Mean\")\n",
    "\n",
    "# set axis labels, make settings\n",
    "ax.set_xlabel(\"ground height [MASL]\")\n",
    "ax.set_ylabel(\"density\")\n",
    "ax.set_title(\"Histogram of Ground Height\")\n",
    "ax.set_ylim(0, 1)\n",
    "ax.set_xlim(100, 300)\n",
    "ax.grid(True, alpha=0.3)\n",
    "l = ax.legend(loc=(0.7, 0.6))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67f66a29",
   "metadata": {},
   "source": [
    "The green line gives the arithmetic mean of all measured ground heights, the red line gives the median. As the median crosses the sum curve at approximately <span>&asymp;</span>140 m, whilst the arithmetic mean of the values is much bigger (<span>&asymp;</span>165 m), we can assume that the ground height is <b>not</b> normally distributed among the area. This indicates, that most of the <i>conventional</i> statistical measures (such as variance, standard deviation, arithmetic mean) are not meaningful, and a more sophisticated statistical evaluation is needed. However, as we are using the geographical position as additional information, the lack of a normal distribution is <b>not</b> problematic for us any more."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b6ea0c2",
   "metadata": {},
   "source": [
    "### Calculate Distances\n",
    "In order to assess the dependency of the ground height from its position (and thus from its proximity of one known point to another), we have to create an empirical variogram. For this, we need to forget about the absolute values and positions of our \"sampling points\" and calculate the distances in space and the height differences between all possible pairs of points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "dbd38506",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# initialize lists to store distance values in\n",
    "distances = []\n",
    "hdiffs = []\n",
    "\n",
    "# iterate through all points (outter loop)\n",
    "# the indexing of the table is as before\n",
    "# with zip(...), we can combine multiple lists; when we then iterate over the zip,\n",
    "#     we get a tuple (value0 from list0, value1 from list1, ...) at each iteration\n",
    "for point1 in zip(df.iloc[:, 2], df.iloc[:, 3], df.iloc[:, 4]):\n",
    "    # make lists to store information in relation to each point in the outter\n",
    "    #     iteration loop\n",
    "    dist_from_point1 = []\n",
    "    hdiff_from_point1 = []\n",
    "    # iterate through all points (inner loop)\n",
    "    for point2 in zip(df.iloc[:, 2], df.iloc[:, 3], df.iloc[:, 4]):\n",
    "        \n",
    "        # calculate the distance in the x-direction dx\n",
    "        dx = point2[0] - point1[0]\n",
    "        \n",
    "        # calculate the distance in the y-direction dy\n",
    "        dy = point2[1] - point1[1]\n",
    "        \n",
    "        # calculate abolute distance\n",
    "        dist = ((dx ** 2) + (dy ** 2)) ** 0.5\n",
    "        \n",
    "        # calculate SQUARED height difference\n",
    "        hdiff = (point1[2] - point2[2]) ** 2\n",
    "        \n",
    "        # append results to dist_from_point1 and hdiff_from_point1\n",
    "        dist_from_point1.append(dist)\n",
    "        hdiff_from_point1.append(hdiff)\n",
    "        \n",
    "    \n",
    "    # append all distances from one point to all\n",
    "    # other points to the distances array\n",
    "    distances.append(dist_from_point1)\n",
    "    hdiffs.append(hdiff_from_point1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08e7880c",
   "metadata": {},
   "source": [
    "### Visual Representation of Distances and Ground Height Differences\n",
    "We can get a overview on the distances between the different sampling points by creating a color bar. The same can be done for the (squared) height differences between the points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "34079a93",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "90d591ee082c41979bec286ffd506483",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# set up a figure with two subplots (one row, two columns)\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(10, 5))\n",
    "\n",
    "# plot representation of distances\n",
    "im1 = ax[0].imshow(distances, cmap=\"rainbow\")\n",
    "# make colorbar\n",
    "cbar1 = plt.colorbar(im1, shrink=0.3, ax=ax[0])\n",
    "# labels etc.\n",
    "cbar1.set_label(\"distance [m]\")\n",
    "ax[0].set_title(\"Distances Between Points\")\n",
    "ax[0].set_xlabel(\"point index\")\n",
    "ax[0].set_ylabel(\"point_index\")\n",
    "\n",
    "# plot represenation of squared height differences\n",
    "im2 = ax[1].imshow(hdiffs, cmap=\"rainbow\")\n",
    "# make colorbar\n",
    "cbar2 = plt.colorbar(im2, shrink=0.3, ax=ax[1])\n",
    "# labels etc.\n",
    "cbar2.set_label(\"squared height difference [m]\")\n",
    "ax[1].set_title(\"Squared Height Difference Between Points\")\n",
    "ax[1].set_xlabel(\"point index\")\n",
    "ax[1].set_ylabel(\"point_index\")\n",
    "\n",
    "# use tight layout to correct for overlapping parts of figure\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b850da2",
   "metadata": {},
   "source": [
    "### Sort Distances\n",
    "In order to create an empirical variogram, we want to rank all pairs of values (i.e., distances and height differences between points), starting from the pairs with the smallest distance between each other. Ultimately, we want to have a data-structure (or table) with the rows being individual pairs of points and the columns being something like \"point 1\", \"point 2\", \"distance\", and \"height difference\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "12e5aee3",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# initialize list to store distances and other information\n",
    "preprocessed_data = []\n",
    "\n",
    "# iterate over upper part of data matrices, excluding the main diagonal\n",
    "#     because matrix is square, symmetric and main diagonal only contains\n",
    "#     zeros in both cases\n",
    "for row in range(0, len(distances) - 1):\n",
    "    for col in range(row + 1, len(distances)):\n",
    "        \n",
    "        # information = [point1, point2, dist, hdiff]\n",
    "        # access distances and hdiffs by row and column indices\n",
    "        information = [row, col, distances[row][col], hdiffs[row][col]]\n",
    "        \n",
    "        # save data\n",
    "        preprocessed_data.append(information)\n",
    "        \n",
    "# put data to DataFrame\n",
    "df_preprocessed = pd.DataFrame(data=preprocessed_data, columns=[\"point1\", \"point2\",\n",
    "                                                                \"distance\", \"hdiff\"])\n",
    "\n",
    "# sort data according to distance\n",
    "df_sorted = df_preprocessed.sort_values(by=[\"distance\"])\n",
    "# reindex DataFrame to start at 0 again and increase because the sorting also\n",
    "#     influences the indices of the rows\n",
    "df_sorted.index = [i for i in range(len(df_sorted))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "00e8a531",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>point1</th>\n",
       "      <th>point2</th>\n",
       "      <th>distance</th>\n",
       "      <th>hdiff</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>489.540601</td>\n",
       "      <td>16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>29</td>\n",
       "      <td>993.227064</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>25</td>\n",
       "      <td>33</td>\n",
       "      <td>1090.875337</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1120.813990</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>29</td>\n",
       "      <td>1182.642803</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   point1  point2     distance  hdiff\n",
       "0       0       1   489.540601     16\n",
       "1       2      29   993.227064      1\n",
       "2      25      33  1090.875337      1\n",
       "3       0       2  1120.813990      1\n",
       "4       0      29  1182.642803      4"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_sorted.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2947648",
   "metadata": {},
   "source": [
    "### Function to Calculate Empirical Variogram Values Based on Number of Classes / Number of Values in each Class\n",
    "The sorted values are summarized in classes and the average distance between all values in one class is defined as the \"class center\". Averaging the height difference in this class allows to plot each class in the empirical variogram: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5514ed5f",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def variogram(data, class_col, var_col, nvals, plotting=True):\n",
    "    \"\"\"\n",
    "    data : pd.DataFrame containing relevant and sorted data, pd.DataFrame\n",
    "    class_col: column to use to calculate class centers, int\n",
    "    var_col : column to use for the variogram data, int\n",
    "    nvals : number of values to consider in each class, int\n",
    "    \"\"\"\n",
    "    \n",
    "    # calculate the number of classes\n",
    "    # math.ceil gets the next larger integer number\n",
    "    nclasses = math.ceil(len(data) / nvals)\n",
    "    \n",
    "    # split array into n different classes\n",
    "    data_split = np.array_split(data, nclasses)\n",
    "    \n",
    "    # get center of class and variogram values\n",
    "    ccenters = []\n",
    "    var_values = []\n",
    "    \n",
    "    # iterate over each class\n",
    "    for c in data_split:\n",
    "        # calculate class center ((min - max) / 2)\n",
    "        ccenter = (c.iloc[-1, class_col] + c.iloc[0, class_col]) / 2\n",
    "        # calculate variogram value\n",
    "        # Note: here, the difference is NOT SQUARED because we already work with\n",
    "        #     squared height differences\n",
    "        var_value = (c.iloc[:, var_col].sum() / (2 * len(c)))\n",
    "        # append the calculated values to the lists\n",
    "        ccenters.append(ccenter)\n",
    "        var_values.append(var_value)\n",
    "    \n",
    "    if plotting:\n",
    "        #Plot the empirical variogram\n",
    "        fig = plt.figure(figsize=(8,6))\n",
    "        ax = fig.add_subplot(1,1,1)\n",
    "        ax.scatter(ccenters, var_values, marker=\"x\", c=\"black\", label=\"empirical variogram\")\n",
    "        # labels etc.\n",
    "        ax.set_xlabel(\"distance [m]\")\n",
    "        ax.set_ylabel(\"variogram $[m^2]$\")\n",
    "        ax.set_title(\"Empirical Variogram\")\n",
    "        ax.grid(True, alpha=0.4, zorder=0)\n",
    "        \n",
    "    return ccenters, var_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "0a1fa9e4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2c9b1986049a4bcf8eb946ecc275a116",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_values_in_class = 30\n",
    "ccenters, var_values = variogram(data=df_sorted, class_col=2, var_col=3, nvals=n_values_in_class)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2761c3da",
   "metadata": {},
   "source": [
    "Visually interpreting the empirical variogram, we can see that the variance among the values that are closer than 2000 m to each other is <1000 m², which is much smaller than the absolute variance of the data set (<span>&asymp;</span>3600 m²). This indicates a strong relationship between values, which are close to each other and means that we can estimate the ground height at a certain point very well, if we know the height at other point(s) in close vicinity. After 6000 m distance, however, the variance within the classes exceeds the absolute variance, sometimes even by far. <br>\n",
    "It might seem strange that some classes (=set of values with a certain distance to each other) show a much higher variance than the overall data set, but you should always keep in mind that our data is <b>not</b> normally distributed and that the choice of sampling points or spatial patterns within the data set can cause higher variations in ground height than the \"average\" variance.<br>\n",
    "However, with this empirical variogram, we can roughly say that the data set allows us to draw good conclusions on all points that are closer than 6000 m to our sampling points. <br>\n",
    "<br>\n",
    "But how can we apply this knowledge for points that are not at the exact distance of our class centers? For this, we will have to find an analytical function, describing our empirical variogram best:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4d05bb2",
   "metadata": {},
   "source": [
    "### Defining Functions to Calculate Residuals Between Empirical and Theoretical Variogram Values\n",
    "For this, we will choose a certain set of theoretical variograms and define, how the height ($h$) depends on the distance ($d$). $\\sigma^2$ is the absolute variance of the ground height, and $\\lambda$, $\\alpha$ and $n$ are plotting parameters.<br>\n",
    "We will use four functions:\n",
    "- an exponential variogram <br>\n",
    "$\\large h(d)=\\sigma^2 \\cdot (1-e^{-d/{\\lambda}}) $\n",
    "- a Gaussian variogram <br>\n",
    "$\\large h(d)=\\sigma^2 \\cdot (1-e^{-(d/\\lambda)^2}) $\n",
    "- a spherical variogram<br>\n",
    "$\\large h(d)=\\sigma^2 \\cdot \\frac{h}{2\\lambda} \\cdot[3-(\\frac{d}{\\lambda})^2] $\n",
    "- a power function variogram<br>\n",
    "$\\large h(d)=n \\cdot d^{\\alpha} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "44cab3a5",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def exp_var(x, class_centers, y_true=None, fit=True):\n",
    "    \"\"\"\n",
    "    x : array [nugget, variance, correlation_length] of parameters, array\n",
    "    class_centers : class centers / x-values to calculate variogram at, array\n",
    "    y_true : observed variogram values (empirical variogram), array\n",
    "    \n",
    "    NOTES: this function takes class centers and observed variogram values\n",
    "    and returns an array of residuals between the theoretical function and\n",
    "    observed values\n",
    "    \"\"\"\n",
    "    \n",
    "    # calculate the theoretical variogram value\n",
    "    exp = [(x[0] + (x[1] - x[0]) * (1 - np.exp(- i / x[2]))) for i in class_centers]\n",
    "    \n",
    "    if fit:\n",
    "        # if fit is True, return an array of residuals (i.e., the deviations\n",
    "        #     of the theoretical variogram values from the empirical variogram\n",
    "        #     values)\n",
    "        return np.array(np.array(exp) - np.array(y_true))\n",
    "    else:\n",
    "        # if fit is False, return only the theoretical variogram values at\n",
    "        #     the given x-values\n",
    "        return exp\n",
    "        \n",
    "\n",
    "def gauss_var(x, class_centers, y_true=None, fit=True):\n",
    "    \"\"\"\n",
    "    x : array [nugget, variance, correlation_length] of parameters, array\n",
    "    class_centers : class centers to calculate variogram at, array\n",
    "    y_true : observed variogram values (empirical variogram), array\n",
    "    \n",
    "    NOTES: this function takes class centers and observed variogram values\n",
    "    and returns an array of residuals between the theoretical function and\n",
    "    observed values\n",
    "    \"\"\"\n",
    "    \n",
    "    # calculate the theoretical variogram value\n",
    "    gauss = [(x[0] + (x[1] - x[0]) * (1 - np.exp(-1 * ((i / x[2]) ** 2)))) for i in class_centers]\n",
    "    \n",
    "    if fit:\n",
    "        # if fit is True, return an array of residuals (i.e., the deviations\n",
    "        #     of the theoretical variogram values from the empirical variogram\n",
    "        #     values)\n",
    "        return np.array(np.array(gauss) - np.array(y_true))\n",
    "    else:\n",
    "        # if fit is False, return only the theoretical variogram values at\n",
    "        #     the given x-values\n",
    "        return gauss\n",
    "\n",
    "def spherical_var(x, class_centers, y_true=None, fit=True):\n",
    "    \"\"\"\n",
    "    x : array [nugget, variance, correlation_length] of parameters, array\n",
    "    class_centers : class centers to calculate variogram at, array\n",
    "    y_true : observed variogram values (empirical variogram), array\n",
    "    \n",
    "    NOTES: this function takes class centers and observed variogram values\n",
    "    and returns an array of residuals between the theoretical function and\n",
    "    observed values\n",
    "    \"\"\"\n",
    "    \n",
    "    # calculate the theoretical variogram value\n",
    "    spherical = [(x[0] + (x[1] - x[0]) * 0.5 * i / x[2] * (3 - (i / x[2]) ** 2)) if i < x[2] else x[1] for i in class_centers]\n",
    "    \n",
    "    if fit:\n",
    "        # if fit is True, return an array of residuals (i.e., the deviations\n",
    "        #     of the theoretical variogram values from the empirical variogram\n",
    "        #     values)\n",
    "        return np.array(np.array(spherical) - np.array(y_true))\n",
    "    else:\n",
    "        # if fit is False, return only the theoretical variogram values at\n",
    "        #     the given x-values\n",
    "        return spherical\n",
    "\n",
    "def power_var(x, class_centers, y_true=None, fit=True):\n",
    "    \"\"\"\n",
    "    x : array [nugget, slope, power] of parameters, array\n",
    "    class_centers : class centers to calculate variogram at, array\n",
    "    y_true : observed variogram values (empirical variogram), array\n",
    "    \n",
    "    NOTES: this function takes class centers and observed variogram values\n",
    "    and returns an array of residuals between the theoretical function and\n",
    "    observed values\n",
    "    \"\"\"\n",
    "    \n",
    "    # calculate the theoretical variogram value\n",
    "    power = [(x[1] * i ** x[2] + x[0]) for i in class_centers]\n",
    "    \n",
    "    if fit:\n",
    "        # if fit is True, return an array of residuals (i.e., the deviations\n",
    "        #     of the theoretical variogram values from the empirical variogram\n",
    "        #     values)\n",
    "        return np.array(np.array(power) - np.array(y_true))\n",
    "    else:\n",
    "        # if fit is False, return only the theoretical variogram values at\n",
    "        #     the given x-values\n",
    "        return power"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2da7a28f-cce5-4cc1-bc60-38989caa2b29",
   "metadata": {},
   "source": [
    "# Manual Fitting of Theoretical Variograms\n",
    "Now, we will fit the functions for our theoretical variograms to our empirical variograms manually (change the number of values in each class and change the sliders!):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e992e365-4f49-4b2f-a275-b6d9e2ab95ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_values_in_class = 30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "c3d6a401-b403-4a9c-8819-f9a21167314d",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d26c5d1cff17442b95507297aa8c136b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8aa7d7e7511f411284acab1e2c490e36",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='Nugget (Exponential)', layout=Layout(width='500px'),…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# get empirical variogram\n",
    "ccenters, var_values = variogram(data=df_sorted, class_col=2, var_col=3, nvals=n_values_in_class, plotting=False)\n",
    "xs = np.arange(0, max(ccenters), 100)\n",
    "\n",
    "# make the base-plot\n",
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "ax.scatter(ccenters, var_values, marker=\"x\", c=\"black\", label=\"empirical variogram\")\n",
    "\n",
    "line_exp, = ax.plot(xs, exp_var(x=[0., 1000., 100.], class_centers=xs, fit=False), label=\"Exponential\", c=\"red\")\n",
    "line_gauss, = ax.plot(xs, gauss_var(x=[0., 1000., 100.], class_centers=xs, fit=False), label=\"Gaussian\", c=\"deepskyblue\")\n",
    "line_spherical, = ax.plot(xs, spherical_var(x=[0., 1000., 100.], class_centers=xs, fit=False), label=\"Spherical\", c=\"limegreen\")\n",
    "line_power, = ax.plot(xs, power_var(x=[0., 1., 1.], class_centers=xs, fit=False), label=\"Power\", c=\"purple\")\n",
    "\n",
    "ax.set_ylim(0, max(var_values) * 1.2)\n",
    "ax.set_xlim(0, max(xs) * 1.1)\n",
    "ax.grid(which=\"both\")\n",
    "ax.set_xlabel(\"Distance\")\n",
    "ax.set_ylabel(\"Variogram\")\n",
    "ax.legend(loc=\"best\")\n",
    "\n",
    "def update(exp_ng=0.,\n",
    "           exp_vr=1000.,\n",
    "           exp_cl=1000.,\n",
    "           gau_ng=0.,\n",
    "           gau_vr=1000.,\n",
    "           gau_cl=1000.,\n",
    "           sph_ng=0.,\n",
    "           sph_vr=1000.,\n",
    "           sph_cl=1000.,\n",
    "           pwr_ng=0.,\n",
    "           pwr_sp=1.,\n",
    "           pwr_pw=1.):\n",
    "    line_exp.set_ydata(exp_var(x=[exp_ng, exp_vr, exp_cl], class_centers=xs, fit=False))\n",
    "    line_gauss.set_ydata(gauss_var(x=[gau_ng, gau_vr, gau_cl], class_centers=xs, fit=False))\n",
    "    line_spherical.set_ydata(spherical_var(x=[sph_ng, sph_vr, sph_cl], class_centers=xs, fit=False))\n",
    "    line_power.set_ydata(power_var(x=[pwr_ng, pwr_sp, pwr_pw], class_centers=xs, fit=False))\n",
    "    fig.canvas.draw_idle()\n",
    "\n",
    "style1 = {'description_width': '200px', \"handle_color\": \"red\"}\n",
    "style2 = {'description_width': '200px', \"handle_color\": \"deepskyblue\"}\n",
    "style3 = {'description_width': '200px', \"handle_color\": \"limegreen\"}\n",
    "style4 = {'description_width': '200px', \"handle_color\": \"purple\"}\n",
    "interact(update,\n",
    "         exp_ng=widgets.FloatSlider(value=0., min=0., max=15000., step=1., description=\"Nugget (Exponential)\", layout=widgets.Layout(width=\"500px\"), style=style1),\n",
    "         exp_vr=widgets.FloatSlider(value=1000., min=0., max=15000., step=1., description=\"Variance (Exponential)\", layout=widgets.Layout(width=\"500px\"), style=style1),\n",
    "         exp_cl=widgets.FloatSlider(value=1000., min=0., max=15000., step=1., description=\"Corr. Len. (Exponential)\", layout=widgets.Layout(width=\"500px\"), style=style1),\n",
    "         gau_ng=widgets.FloatSlider(value=0., min=0., max=15000., step=1., description=\"Nugget (Gaussian)\", layout=widgets.Layout(width=\"500px\"), style=style2),\n",
    "         gau_vr=widgets.FloatSlider(value=1000., min=0., max=15000., step=1., description=\"Variance (Gaussian)\", layout=widgets.Layout(width=\"500px\"), style=style2),\n",
    "         gau_cl=widgets.FloatSlider(value=1000., min=0., max=15000., step=1., description=\"Corr. Len. (Gaussian)\", layout=widgets.Layout(width=\"500px\"), style=style2),\n",
    "         sph_ng=widgets.FloatSlider(value=0., min=0., max=15000., step=1., description=\"Nugget (Spherical)\", layout=widgets.Layout(width=\"500px\"), style=style3),\n",
    "         sph_vr=widgets.FloatSlider(value=1000., min=0., max=15000., step=1., description=\"Variance (Spherical)\", layout=widgets.Layout(width=\"500px\"), style=style3),\n",
    "         sph_cl=widgets.FloatSlider(value=1000., min=0., max=15000., step=1., description=\"Corr. Len. (Spherical)\", layout=widgets.Layout(width=\"500px\"), style=style3),\n",
    "         pwr_ng=widgets.FloatSlider(value=0., min=0., max=15000., step=1., description=\"Nugget (Power)\", layout=widgets.Layout(width=\"500px\"), style=style4),\n",
    "         pwr_sp=widgets.FloatSlider(value=0.1, min=0., max=10., step=.01, description=\"Slope (Power)\", layout=widgets.Layout(width=\"500px\"), style=style4),\n",
    "         pwr_pw=widgets.FloatSlider(value=1., min=0., max=5., step=.01, description=\"Power (Power)\", layout=widgets.Layout(width=\"500px\"), style=style4));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1056410a",
   "metadata": {},
   "source": [
    "### Fitting and Plotting Empirical and Fitted Theoretical Variograms\n",
    "Now, we will fit the functions for our theoretical variograms to our empirical variograms. In the cell below, we define a function to do this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "6c9373ef",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plot_variograms(class_centers, variogram_values, manual_fitting=False,\n",
    "                    exp_ng=1000.,\n",
    "                    exp_vr=1000.,\n",
    "                    exp_cl=1000.,\n",
    "                    gau_ng=1000.,\n",
    "                    gau_vr=1000.,\n",
    "                    gau_cl=1000.,\n",
    "                    sph_ng=1000.,\n",
    "                    sph_vr=1000.,\n",
    "                    sph_cl=1000.,\n",
    "                    pwr_ng=0.,\n",
    "                    pwr_sp=1.,\n",
    "                    pwr_pw=1.):\n",
    "    \"\"\"\n",
    "    class_centers : class centers, array\n",
    "    variogram_values : variogram values, array\n",
    "    manual_fitting : whether to manually fit parameters in the theoretical variograms\n",
    "        or to fit with scipy.least_squares, bool\n",
    "    (exp / gau / sph / pwr)_ng : nugget effect for exponential / gaussian / spherical / power variogram; floats\n",
    "    (exp / gau / sph)_vr : variance for exponential / gaussian / spherical variogram; floats\n",
    "    (exp / gau / sph)_cl : corr. length for exponential / gaussian / spherical variogram; floats\n",
    "    pwr_sp : slope for power variogram; float\n",
    "    pwr_pw : power for power variogram; float\n",
    "    \"\"\"\n",
    "    \n",
    "    # x-values for theoretical variograms\n",
    "    xs = np.arange(0, max(class_centers), 100)\n",
    "    \n",
    "    \"\"\" exponential variogram \"\"\"\n",
    "    if manual_fitting:\n",
    "        # change these parameters for manual fitting\n",
    "        # [nugget, variance, correlation_length]\n",
    "        params_exp = [exp_ng, exp_vr, exp_cl]\n",
    "    else:\n",
    "        # scipy least squares optimization for exp\n",
    "        exp_x0 = np.array([0, 3000, 3000])\n",
    "        result_exp = least_squares(\n",
    "            exp_var,\n",
    "            exp_x0,\n",
    "            kwargs={\"class_centers\": class_centers, \"y_true\": variogram_values},\n",
    "            bounds=([0, 0, 0], [1000, 10000, 10000])\n",
    "        )\n",
    "        params_exp = result_exp.x\n",
    "        resids = result_exp.fun ** 2\n",
    "        ssres = resids.sum()\n",
    "        sstot = np.array([(i - np.array(var_values).mean()) ** 2 for i in var_values]).sum()\n",
    "        r2_exp = 1 - (ssres / sstot)\n",
    "    \n",
    "    # fitted exponential variogram    \n",
    "    exp_fit = [(params_exp[0] + (params_exp[1] - params_exp[0]) * (1 - np.exp(- i / params_exp[2]))) for i in xs]\n",
    "    \n",
    "    \"\"\" gaussian variogram \"\"\"\n",
    "    if manual_fitting:\n",
    "        # change these parameters for manual fitting\n",
    "        # [nugget, variance, correlation_length]\n",
    "        params_gauss = [gau_ng, gau_vr, gau_cl]\n",
    "    else:\n",
    "        # scipy least squares optimization for gauss\n",
    "        gauss_x0 = np.array([0, 3000, 3000])\n",
    "        result_gauss = least_squares(\n",
    "            gauss_var,\n",
    "            gauss_x0,\n",
    "            kwargs={\"class_centers\": class_centers, \"y_true\": variogram_values},\n",
    "            bounds=([0, 0, 0], [1000, 10000, 10000])\n",
    "        )\n",
    "        params_gauss = result_gauss.x\n",
    "        resids = result_gauss.fun ** 2\n",
    "        ssres = resids.sum()\n",
    "        sstot = np.array([(i - np.array(var_values).mean()) ** 2 for i in var_values]).sum()\n",
    "        r2_gauss = 1 - (ssres / sstot)\n",
    "    \n",
    "    # fitted gauss variogram    \n",
    "    gauss_fit = [(params_gauss[0] + (params_gauss[1] - params_gauss[0]) * (1 - np.exp(-1 * ((i / params_gauss[2]) ** 2)))) for i in xs]\n",
    "    \n",
    "    \"\"\" shperical variogram \"\"\"\n",
    "    if manual_fitting:\n",
    "        # change these parameters for manual fitting\n",
    "        # [nugget, variance, correlation_length]\n",
    "        params_spherical = [sph_ng, sph_vr, sph_cl]\n",
    "    else:\n",
    "        # scipy least squares optimization for spherical\n",
    "        spherical_x0 = np.array([0, 3000, 3000])\n",
    "        result_spherical = least_squares(\n",
    "            spherical_var,\n",
    "            spherical_x0,\n",
    "            kwargs={\"class_centers\": class_centers, \"y_true\": variogram_values},\n",
    "            bounds=([0, 0, 0], [1000, 10000, 10000])\n",
    "        )\n",
    "        params_spherical = result_spherical.x\n",
    "        resids = result_spherical.fun ** 2\n",
    "        ssres = resids.sum()\n",
    "        sstot = np.array([(i - np.array(var_values).mean()) ** 2 for i in var_values]).sum()\n",
    "        r2_spherical = 1 - (ssres / sstot)\n",
    "    \n",
    "    # fitted spherical variogram    \n",
    "    spherical_fit = [(params_spherical[0] + (params_spherical[1] - params_spherical[0]) * 0.5 * i / params_spherical[2] * \\\n",
    "                      (3 - (i / params_spherical[2]) ** 2)) if i < params_spherical[2] else params_spherical[1] for i in xs]\n",
    "    \n",
    "    \"\"\" power variogram \"\"\"\n",
    "    if manual_fitting:\n",
    "        # change these parameters for manual fitting\n",
    "        # [nugget, slope, power]\n",
    "        params_power = [pwr_ng, pwr_sp, pwr_pw]\n",
    "    else:\n",
    "        # scipy least squares optimization for power\n",
    "        power_x0 = np.array([0, 1, 1])\n",
    "        result_power = least_squares(\n",
    "            power_var,\n",
    "            power_x0,\n",
    "            kwargs={\"class_centers\": class_centers, \"y_true\": variogram_values},\n",
    "            bounds=([0, 0, 0], [1000, 10000, 10000])\n",
    "        )\n",
    "        params_power = result_power.x\n",
    "        resids = result_power.fun ** 2\n",
    "        ssres = resids.sum()\n",
    "        sstot = np.array([(i - np.array(var_values).mean()) ** 2 for i in var_values]).sum()\n",
    "        r2_power = 1 - (ssres / sstot)\n",
    "    \n",
    "    # fitted power variogram    \n",
    "    power_fit = [(params_power[0] + (params_power[1] - params_power[0]) * i ** params_power[2]) for i in xs]\n",
    "    \n",
    "    \n",
    "    fig, ax = plt.subplots(figsize=(12, 6))\n",
    "    ax.scatter(class_centers, variogram_values, marker=\"x\", c=\"black\", label=\"empirical variogram\")\n",
    "    \n",
    "    ax.plot(xs, exp_fit, label=\"Exponential\", c=\"red\")\n",
    "    ax.plot(xs, gauss_fit, label=\"Gaussian\", c=\"deepskyblue\")\n",
    "    ax.plot(xs, spherical_fit, label=\"Spherical\", c=\"limegreen\")\n",
    "    ax.plot(xs, power_fit, label=\"Power\", c=\"purple\")\n",
    "    \n",
    "    # print out the fitted parameters\n",
    "    ax.text(1.05, 0.8, \"Exponential: \\n  nugget = %1.2f \\n  variance = %1.2f \\n  corr. length = %1.2f \\n  $R^2$ = %1.2f\" \\\n",
    "            %(params_exp[0], params_exp[1], params_exp[2], r2_exp), transform=ax.transAxes, fontsize=12)\n",
    "    \n",
    "    ax.text(1.05, 0.55, \"Gaussian: \\n  nugget = %1.2f \\n  variance = %1.2f \\n  corr. length = %1.2f \\n  $R^2$ = %1.2f\" \\\n",
    "            %(params_gauss[0], params_gauss[1], params_gauss[2], r2_gauss), transform=ax.transAxes, fontsize=12)\n",
    "    \n",
    "    ax.text(1.05, 0.3, \"Spherical: \\n  nugget = %1.2f \\n  variance = %1.2f \\n  corr. length = %1.2f  \\n  $R^2$ = %1.2f\" \\\n",
    "            %(params_spherical[0], params_spherical[1], params_spherical[2], r2_spherical), transform=ax.transAxes, fontsize=12)\n",
    "    \n",
    "    ax.text(1.05, 0.05, \"Power: \\n  nugget = %1.2f \\n  slope = %1.2f \\n  power = %1.2f \\n  $R^2$ = %1.2f\" \\\n",
    "            %(params_power[0], params_power[1], params_power[2], r2_power), transform=ax.transAxes, fontsize=12)\n",
    "    \n",
    "    ax.legend(loc=\"best\", fontsize=14)\n",
    "    ax.set_xlabel(\"distance [m]\", fontsize=14)\n",
    "    ax.set_ylabel(\"variogram $[m^2]$\", fontsize=14)\n",
    "    ax.set_title(\"Empirical and Fitted Theoretical Variograms\", fontsize=16)\n",
    "    ax.set_ylim(0, 8500)\n",
    "    ax.set_xlim(0, 20000)\n",
    "    ax.grid(True, alpha=0.4, zorder=0)\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    \n",
    "    return result_exp.x, result_gauss.x, result_spherical.x, result_power.x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "872b53ea",
   "metadata": {},
   "source": [
    "### Creating Empirical Variogram Values and Fit Theoretical Variograms\n",
    "\n",
    "- try different number of values in the classes\n",
    "- try different initial conditions for theoretical variograms\n",
    "- try different parameter bounds for theoretical variograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7e4a8459",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6dd4f34963a747c8b6df395ce608aeeb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ccenters, var_values = variogram(data=df_sorted, class_col=2, var_col=3, nvals=100, plotting=False)\n",
    "p_exp, p_gauss, p_spherical, p_power = plot_variograms(class_centers=ccenters, variogram_values=var_values, manual_fitting=False)"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}