Added course notebooks

class_2023/sgm_2023_01b_interpolation_and_inverse_distance_student_version.ipynb

@@ -0,0 +1,365 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Linear interpolation and inverse distance weighting (IDW) in 1-D\n",
+    "\n",
+    "Similar to the example of regression, we will briefly evaluate linear interpolation and IDW 1-D. This is a very simplified example - but will show us the basic principles that we can then extend to higher dimensions - and understand why some things work and some may not.\n",
+    "\n",
+    "This notebook is part of the class \"Structural Geological Models\".\n",
+    "\n",
+    "(c) Florian Wellmann, CG3, 2023\n",
+    "\n",
+    "## Brief review of the theory\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# basic imports\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import matplotlib.image as mpimg # needed to load the background image"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Input data\n",
+    "\n",
+    "We will use again the image from the lecture with the well markers from the virtual drillholes. You can use the same input points as before."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Set the position of the input points and plot them above the image:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# x-positions of points:\n",
+    "x_positions = []\n",
+    "\n",
+    "# y-positions of points:\n",
+    "y_positions = []"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import matplotlib.image as mpimg\n",
+    "\n",
+    "# Load your image using matplotlib's imread\n",
+    "image_path = 'WandC_boreholes_2.png'  # Replace with the path to your image file\n",
+    "img = mpimg.imread(image_path)\n",
+    "\n",
+    "# Determine the aspect ratio of the image\n",
+    "image_aspect_ratio = img.shape[1] / img.shape[0]\n",
+    "\n",
+    "# Create a figure with the aspect ratio of the image\n",
+    "fig, ax = plt.subplots(figsize=(12, 12 / image_aspect_ratio))\n",
+    "\n",
+    "# Display the image with the same aspect ratio\n",
+    "ax.imshow(img, extent=[0, 20, 0, 10], aspect='auto')\n",
+    "\n",
+    "# Set limits, labels, title, etc.\n",
+    "ax.set_xlim(0, 20)\n",
+    "ax.set_ylim(0, 10)\n",
+    "ax.set_xlabel('X')\n",
+    "ax.set_ylabel('Y')\n",
+    "ax.set_title('Markers in virtual boreholes')\n",
+    "\n",
+    "plt.grid()\n",
+    "\n",
+    "# Example data plot\n",
+    "ax.plot(x_positions, y_positions, 'o', markersize=15, color='blue', label='Well markers')\n",
+    "ax.legend();\n",
+    "\n",
+    "# Show the plot\n",
+    "plt.show()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Adjust the positions until you receive a good match to the markers."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Linear interpolation\n",
+    "\n",
+    "Viewing linear interpolation results is actually almost too easy in matplotlib: a line connecting input points is default behavior. But for a repetition of theory, first implement a function to calculate the Lagrange interpolating polynomial between a set of input points.\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from scipy.interpolate import lagrange\n",
+    "\n",
+    "def lagrange_interpolation(x_points, y_points):\n",
+    "    \"\"\"\n",
+    "    Calculate and plot the Lagrange interpolating polynomial.\n",
+    "\n",
+    "    :param x_points: List of x-coordinates of the points.\n",
+    "    :param y_points: List of y-coordinates of the points.\n",
+    "    \"\"\"\n",
+    "    # Calculate the Lagrange interpolating polynomial\n",
+    "    polynomial = lagrange(x_points, y_points)\n",
+    "\n",
+    "    # Generate a dense range of x values for plotting\n",
+    "    x_dense = np.linspace(0, 20, 500)\n",
+    "\n",
+    "    # Evaluate the polynomial at the dense x values\n",
+    "    y_dense = polynomial(x_dense)\n",
+    "\n",
+    "    return x_dense, y_dense\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note: by default, the degree of the polynomial is chosen to be one less than the number of provided input points\n",
+    "\n",
+    "**Your task**: test this function with several sets of input points and plot the results in the combined figure:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "x_dense, y_dense = lagrange_interpolation(x_positions[:3], y_positions[:3])\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(12, 12 / image_aspect_ratio))\n",
+    "\n",
+    "# Display the image with the same aspect ratio\n",
+    "ax.imshow(img, extent=[0, 20, 0, 10], aspect='auto')\n",
+    "\n",
+    "# Set limits, labels, title, etc.\n",
+    "ax.set_xlim(0, 20)\n",
+    "ax.set_ylim(0, 10)\n",
+    "ax.set_xlabel('X')\n",
+    "ax.set_ylabel('Y')\n",
+    "ax.set_title('Markers in virtual boreholes')\n",
+    "\n",
+    "plt.grid()\n",
+    "\n",
+    "# Plot the original points\n",
+    "plt.scatter(x_positions, y_positions, color='red', label='Original Points')\n",
+    "\n",
+    "# Plot the Lagrange interpolating polynomial\n",
+    "plt.plot(x_dense, y_dense, label='Lagrange Polynomial')\n",
+    "\n",
+    "# Add labels and legend\n",
+    "plt.xlabel('X')\n",
+    "plt.ylabel('Y')\n",
+    "plt.title('Lagrange Interpolating Polynomial')\n",
+    "plt.legend()\n",
+    "\n",
+    "# Show the plot\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Experiment with data and interpolation\n",
+    "\n",
+    "Test the following aspects to get a better insight into the results:\n",
+    "- Test the interpolation for different sets of input points. As in the regression example: calculate the difference between prediction and actual values at the points that are left out."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Inverse Distance Weighting (IDW)\n",
+    "\n",
+    "Next, we perform an interpolation using inverse distance weighting. Here, also in 1-D - but as we discussed in the lecture, the extension of the concept to 2-D is fairly straight-forward.\n",
+    "\n",
+    "Inverse Distance Weighting (IDW) is a method used for interpolating values at unknown points based on values from known points. The fundamental idea behind IDW is that the estimated value at an unknown point should be influenced more by closer known points than by those farther away. This influence typically decreases with distance.\n",
+    "\n",
+    "### Theory\n",
+    "\n",
+    "The IDW interpolated value $I(x)$ at a point $x$ can be computed using the formula:\n",
+    "\n",
+    "$$\n",
+    "I(x) = \\frac{\\sum_{i=1}^{n} \\frac{1}{d(x, x_i)^p} \\cdot y_i}{\\sum_{i=1}^{n} \\frac{1}{d(x, x_i)^p}}\n",
+    "$$\n",
+    "\n",
+    "where:\n",
+    "- $n$ is the number of known points,\n",
+    "- $x_i$ are the known points,\n",
+    "- $y_i$ are the values at the known points,\n",
+    "- $d(x, x_i)$ is the distance between the unknown point $x$ and a known point $x_i$,\n",
+    "- $p$ is a power parameter that determines how quickly the influence of known points decreases with distance. A larger value of $p$ gives greater influence to nearer points.\n",
+    "\n",
+    "In the case of 1-D interpolation, the distance $d(x, x_i)$ is simply the absolute difference between $x$ and $x_i$:\n",
+    "\n",
+    "$$\n",
+    "d(x, x_i) = |x - x_i|\n",
+    "$$\n",
+    "\n",
+    "It's important to handle the case where the unknown point $x$ coincides with one of the known points $x_i$. In this scenario, the interpolated value $I(x)$ should be exactly equal to the value at the known point, $y_i$.\n",
+    "\n",
+    "Shepard's method for smoothing is a specific implementation of IDW interpolation where the power parameter $p$ is typically set to 2, although it can be adjusted for different smoothing effects.\n",
+    "\n",
+    "### Implementation\n",
+    "\n",
+    "Here is a simple implementation in Python:\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Function for Inverse Distance Weighting using Shepard's method\n",
+    "def inverse_distance_weighting(x, xi, values, power=2):\n",
+    "    \"\"\"\n",
+    "    Apply Inverse Distance Weighting interpolation.\n",
+    "\n",
+    "    :param x: Array of points where interpolation is computed.\n",
+    "    :param xi: Known data points.\n",
+    "    :param values: Values at the known data points.\n",
+    "    :param power: Power parameter for weights. Higher value means more influence to nearer points.\n",
+    "    :return: Interpolated values at points x.\n",
+    "    \"\"\"\n",
+    "    # Initialize array to store interpolated values\n",
+    "    interpolated_values = np.zeros(len(x))\n",
+    "\n",
+    "    # Iterate over each point in x\n",
+    "    for i in range(len(x)):\n",
+    "        # Compute weights based on inverse distance\n",
+    "        weights = 1 / (np.abs(x[i] - xi) ** power)\n",
+    "        \n",
+    "        # Handle case where x[i] is exactly one of the xi (to avoid division by zero)\n",
+    "        if np.any(weights == np.inf):\n",
+    "            interpolated_values[i] = values[np.argmax(weights)]\n",
+    "        else:\n",
+    "            # Calculate the weighted average\n",
+    "            interpolated_values[i] = np.sum(weights * values) / np.sum(weights)\n",
+    "\n",
+    "    return interpolated_values\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# set the power value (Shepards method):\n",
+    "power = 1\n",
+    "\n",
+    "# Apply IDW interpolation\n",
+    "y_dense = inverse_distance_weighting(x_dense, x_positions, y_positions, power=power)\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(12, 12 / image_aspect_ratio))\n",
+    "\n",
+    "# Display the image with the same aspect ratio\n",
+    "ax.imshow(img, extent=[0, 20, 0, 10], aspect='auto')\n",
+    "\n",
+    "# Set limits, labels, title, etc.\n",
+    "ax.set_xlim(0, 20)\n",
+    "ax.set_ylim(0, 10)\n",
+    "ax.set_xlabel('X')\n",
+    "ax.set_ylabel('Y')\n",
+    "ax.set_title('Markers in virtual boreholes')\n",
+    "\n",
+    "plt.grid()\n",
+    "\n",
+    "# Plot the original points\n",
+    "plt.scatter(x_positions, y_positions, color='red', label='Original Points')\n",
+    "\n",
+    "# Plot the Lagrange interpolating polynomial\n",
+    "plt.plot(x_dense, y_dense, label='Lagrange Polynomial')\n",
+    "\n",
+    "# Add labels and legend\n",
+    "plt.xlabel('X')\n",
+    "plt.ylabel('Y')\n",
+    "plt.title('Lagrange Interpolating Polynomial')\n",
+    "plt.legend()\n",
+    "\n",
+    "# Show the plot\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Experiment with data and interpolation\n",
+    "\n",
+    "In order to get a feeling for the interpolation, experiment with the settings:\n",
+    "- Try to make the points more extreme (changing the position of well-markers, for example): how robust is the result, compared to regression and simple interpolation before?\n",
+    "- As before: remove points and calculate the (vertical) error of the prediction.\n",
+    "- Test the influence of the `power`-parameter on the results (i.e.: smoothing around the data points)\n",
+    "- What happens if you set the `power`-parameter to a high value? What does the result look like?"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "base",
+   "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.9.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}