{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### X lines of Python\n", "\n", "# Wedge model\n", "\n", "This is part of [an Agile blog series](http://ageo.co/xlines00) called **x lines of Python**.\n", "\n", "We start with the usual preliminaries." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Make an earth model\n", "\n", "We'll start off with an earth model --- an array of 'cells', each of which has some rock properties.\n", "\n", "Line 1 sets up some basic variables, then in line 2 I've used a little matrix-forming trick, `np.tri(m, n, k)`, which creates an *m* × *n* matrix with ones below the *k*th diagonal, and zeros above it. The `dtype` specification just makes sure we end up with integers, which we need later for the indexing trick.\n", "\n", "Then line 3 just sets every row above `depth//3` (the `//` is integer division, because NumPy prefers integers for indexing arrays), to 0." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "length, depth = 40, 100\n", "model = 1 + np.tri(depth, length, -depth//3, dtype=int)\n", "model[:depth//3,:] = 0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll have a quick look with some very basic plotting commands." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAADKCAYAAAC11LviAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADNNJREFUeJzt3V+spHV9x/H3p/uHlT+GPy6E7mIBQ1qJrUhOkZaGGLAN\nWtOlCaTY1mwakq2NWqg2ZeUG2pQEmlb0SrMN6F6gQFZauDBaipC2NysLruXPiqyosOyWXS1UegOC\n317MQ3qynNmZs2fOeeb8zvuVnJyZZ55hPvmx57Nffs/MIVWFJGn5+4W+A0iSJsNCl6RGWOiS1AgL\nXZIaYaFLUiMsdElqhIUuSY1YUKEnuSzJU0n2Jtk6qVCSpPnL0X6wKMkq4HvAbwP7gIeBD1fVk5OL\nJ0ka1+oFPPcCYG9VPQOQ5E5gEzC00NfmmFrHcQt4SUlaeV7mxR9X1fpR5y2k0DcAz826vw947+En\nJdkCbAFYx7G8N5cu4CUlaeX519rxo3HOW8geeuY49qb9m6raVlUzVTWzhmMW8HKSpCNZSKHvA86Y\ndX8jsH9hcSRJR2shhf4wcE6Ss5KsBa4C7ptMLEnSfB31HnpVvZbk48A3gFXA7VX1xMSSSZLmZSEX\nRamqrwFfm1AWSdIC+ElRSWqEhS5JjbDQJakRFrokNcJCl6RGWOiS1AgLXZIaYaFLUiMsdElqhIUu\nSY2w0CWpERa6JDXCQpekRljoktQIC12SGmGhS1IjLHRJaoSFLkmNsNAlqREWuiQ1wkKXpEZY6JLU\nCAtdkhphoUtSI0YWepIzkjyYZE+SJ5Jc0x0/Ocn9SZ7uvp+0+HElScOMM6G/Bnyqqt4JXAh8LMm5\nwFbggao6B3iguy9J6snIQq+qA1X1aHf7ZWAPsAHYBGzvTtsOXL5YISVJo81rDz3JmcB7gJ3AaVV1\nAAalD5w66XCSpPGNXehJjge+ClxbVT+dx/O2JNmVZNfPeOVoMkqSxjBWoSdZw6DM76iqe7rDLyQ5\nvXv8dODgXM+tqm1VNVNVM2s4ZhKZJUlzGOddLgFuA/ZU1WdmPXQfsLm7vRm4d/LxJEnjWj3GORcB\nHwEeS7K7O3Y9cDNwd5KrgWeBKxcnoiRpHCMLvar+A8iQhy+dbBxJ0tHyk6KS1AgLXZIaYaFLUiMs\ndElqhIUuSY2w0CWpERa6JDXCQpekRljoktQIC12SGmGhS1IjLHRJaoSFLkmNsNAlqREWuiQ1Ypz/\nwcXEvHLGcez91IVL+ZKStPxdu2Os05zQJakRSzqh/+pJh/jWH3xh6OPvuOujS5hGktrihC5JjVjS\nCX2U7x9hen+DU7wkzW2qCn0co0rfwpe0UrnlIkmNWHYT+ihO8JJWKid0SWrE2BN6klXALuD5qvpQ\nkrOAO4GTgUeBj1TVq4sTc3K88CqpVfOZ0K8B9sy6fwtwa1WdA7wIXD3JYJKk+RlrQk+yEfhd4Cbg\nk0kCXAL8YXfKduBG4POLkHHJuQ8vaTkad8vls8BfASd0908BXqqq17r7+4ANE842tdy2kTSNRm65\nJPkQcLCqHpl9eI5Ta8jztyTZlWTXoZ+8fpQxJUmjjDOhXwT8XpIPAuuAtzKY2E9Msrqb0jcC++d6\nclVtA7YBzLx73Zyl3yK3bSQttZETelV9uqo2VtWZwFXAN6vqj4AHgSu60zYD9y5aSknSSAv5YNF1\nwJ1J/hb4NnDbZCKtDE7wkiZtXoVeVQ8BD3W3nwEumHwkgRdeJc2fnxSVpEY097tcVhK3bSTN5oQu\nSY1wQm+YE7y0sljoK5gXXqW2uOUiSY1wQtcRuW0jLR9O6JLUCCd0LYj78NL0sNC16Ny2kZaGWy6S\n1AgndPXOCV6aDCd0SWqEE7qmnhdepfFY6GqC2zaSWy6S1AwndK0ITvBaCZzQJakRTugSXnhVGyx0\naUxu22jaueUiSY1wQpcmxG0b9c0JXZIa4YQuLSH34bWYLHRpilj4WoixtlySnJhkR5LvJtmT5DeS\nnJzk/iRPd99PWuywkqThxp3QPwd8vaquSLIWOBa4Hnigqm5OshXYCly3SDkl4YVXHdnICT3JW4GL\ngdsAqurVqnoJ2ARs707bDly+WCElSaONM6GfDRwCvpjk3cAjwDXAaVV1AKCqDiQ5dfFiShqX+/Ar\n1ziFvho4H/hEVe1M8jkG2ytjSbIF2ALw9g1eg5X6ZuG3a5yLovuAfVW1s7u/g0HBv5DkdIDu+8G5\nnlxV26pqpqpm1p+yahKZJUlzGDkyV9V/JXkuyS9X1VPApcCT3ddm4Obu+72LmlTSkvDC6/I17h7I\nJ4A7une4PAP8CYPp/u4kVwPPAlcuTkRJ0jjGKvSq2g3MzPHQpZONI2k5cB9+Ovm7XCSpEb7tRNLE\nuQ/fDwtdUi/ctpk8t1wkqRFO6JKmkhP8/DmhS1IjnNAlLUteeH0zC11Ss1bato1bLpLUCCd0SStW\naxO8E7okNcIJXZKGWG4XXi10SVqAadq2cctFkhrhhC5Ji2gpt22c0CWpEU7oktSzUVP8qmvH++c4\noUtSIyx0SWqEhS5JjbDQJakRFrokNcJCl6RGWOiS1IixCj3JXyR5IsnjSb6SZF2Ss5LsTPJ0kruS\nrF3ssJKk4UYWepINwJ8DM1X1LmAVcBVwC3BrVZ0DvAhcvZhBJUlHNu6Wy2rgLUlWA8cCB4BLgB3d\n49uByycfT5I0rpGFXlXPA38PPMugyP8HeAR4qape607bB2xYrJCSpNHG2XI5CdgEnAX8InAc8IE5\nTq0hz9+SZFeSXYd+8vpCskqSjmCcLZf3Az+oqkNV9TPgHuA3gRO7LRiAjcD+uZ5cVduqaqaqZtaf\nsmoioSVJbzZOoT8LXJjk2CQBLgWeBB4ErujO2QzcuzgRJUnjGGcPfSeDi5+PAo91z9kGXAd8Msle\n4BTgtkXMKUkaYazfh15VNwA3HHb4GeCCiSeSJB0VPykqSY2w0CWpERa6JDXCQpekRljoktQIC12S\nGmGhS1IjLHRJaoSFLkmNsNAlqREWuiQ1wkKXpEZY6JLUCAtdkhphoUtSIyx0SWqEhS5JjbDQJakR\nFrokNcJCl6RGWOiS1AgLXZIaYaFLUiMsdElqhIUuSY2w0CWpERa6JDXCQpekRqSqlu7FkkPAj2Yd\nehvw4yULcPTMOVnLIedyyAjmnLRpzflLVbV+1ElLWuhvevFkV1XN9BZgTOacrOWQczlkBHNO2nLJ\nOYxbLpLUCAtdkhrRd6Fv6/n1x2XOyVoOOZdDRjDnpC2XnHPqdQ9dkjQ5fU/okqQJ6a3Qk1yW5Kkk\ne5Ns7SvHkST5YZLHkuxOsqvvPLMluT3JwSSPzzp2cpL7kzzdfT9pCjPemOT5bk13J/lgnxm7TGck\neTDJniRPJLmmOz5t6zks51StaZJ1Sb6V5Dtdzr/ujp+VZGe3nnclWTuFGb+U5Aez1vK8vjIelapa\n8i9gFfB94GxgLfAd4Nw+sozI+UPgbX3nGJLtYuB84PFZx/4O2Nrd3grcMoUZbwT+su/1Oyzn6cD5\n3e0TgO8B507heg7LOVVrCgQ4vru9BtgJXAjcDVzVHf8C8GdTmPFLwBV9r+HRfvU1oV8A7K2qZ6rq\nVeBOYFNPWZalqvo34L8PO7wJ2N7d3g5cvqShDjMk49SpqgNV9Wh3+2VgD7CB6VvPYTmnSg38b3d3\nTfdVwCXAju54r+t5hIzLWl+FvgF4btb9fUzhH0wG/4L/JckjSbb0HWYMp1XVARj88AOn9pxnmI8n\n+c9uS6bXbYzDJTkTeA+DiW1q1/OwnDBla5pkVZLdwEHgfgb/Rf5SVb3WndL7z/zhGavqjbW8qVvL\nW5Mc02PEeeur0DPHsWn82/Giqjof+ADwsSQX9x2oAZ8H3gGcBxwA/qHfOP8vyfHAV4Frq+qnfecZ\nZo6cU7emVfV6VZ0HbGTwX+TvnOu0pU112IsfljHJu4BPA78C/DpwMnBdjxHnra9C3wecMev+RmB/\nT1mGqqr93feDwD8x+IM5zV5IcjpA9/1gz3nepKpe6H6Qfg78I1OypknWMCjJO6rqnu7w1K3nXDmn\ndU0Bquol4CEG+9MnJlndPTQ1P/OzMl7WbWtVVb0CfJEpWstx9FXoDwPndFe91wJXAff1lGVOSY5L\ncsIbt4HfAR4/8rN6dx+wubu9Gbi3xyxzeqMgO7/PFKxpkgC3AXuq6jOzHpqq9RyWc9rWNMn6JCd2\nt98CvJ/Bfv+DwBXdab2u55CM3531F3gY7PH3/udzPnr7YFH31qrPMnjHy+1VdVMvQYZIcjaDqRxg\nNfDlacqY5CvA+xj8drgXgBuAf2bwToK3A88CV1ZVbxclh2R8H4OtgWLwLqI/fWOfui9Jfgv4d+Ax\n4Ofd4esZ7E9P03oOy/lhpmhNk/wag4ueqxgMjXdX1d90P1N3MtjK+Dbwx90kPE0ZvwmsZ7AtvBv4\n6KyLp1PPT4pKUiP8pKgkNcJCl6RGWOiS1AgLXZIaYaFLUiMsdElqhIUuSY2w0CWpEf8Hxjh98qev\nSIcAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.imshow(model, cmap='viridis', aspect=0.2)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,\n", " 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model[60]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can make some Vp-rho pairs (rock 0, rock 1, rock 2) and select from those with `np.take`. This works like `vlookup` in Excel --- it says \"read this array, `model` in this case, in which the values *i* are like 0, 1, ... n, and give me the *i*th element from this other array, `rocks` in this case." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "rocks = np.array([[2700, 2750], # Vp, rho\n", " [2400, 2450],\n", " [2800, 3000]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Edit:** I was using `np.take` here, but ['fancy indexing'](http://docs.scipy.org/doc/numpy/user/basics.indexing.html) is shorter and more intuitive. We are just going to index `rocks` using the integers in `model`. That is, if `model` has a `1`, we take the second element, `[2400, 2450]`, from `rocks`. We'll end up with an array containing the rocks corresponding to each element of `earth`." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "earth = rocks[model]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now apply `np.product` to those Vp-rho pairs to get impedance at every sample.\n", "\n", "This might look a bit magical, but we're just telling Python to apply the function `product()` to every set of numbers it encounters on the last axis (index `-1`). The array `earth` has shape (100, 40, 2), so you can think of it as a 100 row x 40 column 'section' in which each 'sample' is occupied by a Vp-rho pair. That pair is in the last axis. So product, which just takes a bunch of numbers and multiplies them, will return the impedance (the product of Vp and rho) at each sample location. We'll end up with a new 100 x 40 'section' with impedance at every sample." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "imp = np.apply_along_axis(np.product, -1, earth)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could have saved a step by taking from `np.product(rocks, axis=1)` but I like the elegance of having an earth model with a set of rock properties at each sample location. That's how I think about the earth --- and it's similar to the concept of a geocellular model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model seismic reflections\n", "\n", "Now we have an earth model — giving us acoustic impedance everywhere in this 2D grid — we define a function to compute reflection coefficients for every trace.\n", "\n", "I love this indexing trick though I admit it looks weird the first time you see it. It's easier to appreciate for a 1D array. Let's look at the differences:\n", "\n", " >>> a = np.array([1,1,1,2,2,2,3,3,3])\n", " >>> a[1:] - a[:-1]\n", " array([0, 0, 1, 0, 0, 1, 0, 0])\n", "\n", "This is equivalent to:\n", "\n", " >>> np.diff(a, axis=0)\n", " \n", "But I prefer to spell it out so it's analogous to the sum on the denominator." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAADGCAYAAAAZrdG7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGZxJREFUeJzt3W1sY9l93/Hvn6REiqQoUhQl7ezsbjYN6qZpVddFDLRO\nYRdOCncCZ40W3TqOATsL1H3hxG5jFLv2GzsvAmQNtIHf9EUbJ9hkbGS3BlIbwbR2Fo5RNGjideut\nNrbHG+/szmpnRk+kREl8EkWevhDvMSmRQ83oiVfz+wAH99w7nOHRtec3d8859xxzziEiIuEQOe8G\niIjI0Sm0RURCRKEtIhIiCm0RkRBRaIuIhIhCW0QkRI4V2mb2PjO7bmavmtnTJ9UoERHpz+53nraZ\nRYBXgfcCt4GXgA86566fXPNERKTbcZ603wn8tXPupnOuCfwR8MTJNEtERPo5Tmg/DCx1nb/VuSYi\nIqckdtpfYGZ6T15E5D445+zgteOE9i3g0a7zy51rh3zoQx/y9ZWVFT71qU8d42vPxtWrV/nwhz98\n3s0YSu08WWFoZxjaCGrnvVpcXGRxcdGff/nLX+77ueOE9kvAT5nZY8Ad4IPAL/f7YPcNuXr16jG+\nUkTkYlpYWGBhYcGfn3hoO+daZvZrwDfY7xv/onPuB/f754mIyHDH6tN2zv0P4G338nu6/yUZZWrn\nyVI7T04Y2ghq52m573naR/4CM3ft2rVT/Q4RkYvmypUrfQci9Rq7iEiIKLRFREJEoS0iEiIKbRGR\nEFFoi4iEiEJbRCREFNoiIiGi0BYRCRGFtohIiCi0RURCRKEtIhIiCm0RkRBRaIuIhIhCW0QkRBTa\nIiIhotAWEQkRhbaISIgotEVEQkShLSISIgptEZEQUWiLiISIQltEJEQU2iIiIaLQFhEJkaGhbWaX\nzeybZvY9M3vFzD7RuZ4zs2+Y2Q/N7OtmNnX6zRURebAd5Ul7D/gN59zPAP8Q+LiZ/S3gGeBF59zb\ngG8Cnz69ZoqICBwhtJ1zy865lzv1HeAHwGXgCeC5zseeAz5wWo0UEZF999SnbWY/Abwd+Atgzjm3\nAvvBDsyedONERKRX7KgfNLM08BXgk865HTNzBz5y8Ny7evWqry8sLLCwsHCv7RQRudAWFxdZXFwc\n+jlzbmDW/vhDZjHgT4D/7pz7QufaD4D3OOdWzGwe+DPn3E/3+b3u2rVr99p+EZEH2pUrV3DO2cHr\nR+0e+T3g+0Fgd3wN+Gin/hHgq8dqoYiIDDW0e8TM3gX8CvCKmX2X/W6QzwDPAi+Y2VPATeDJ02yo\niIgcIbSdc38ORAf88s+fbHNERORu9EakiEiIKLRFREJEoS0iEiIKbRGREFFoi4iEiEJbRCREFNoi\nIiGi0BYRCRGFtohIiCi0RURCRKEtIhIiCm0RkRBRaIuIhMiRd645jqkpbdQuInISjrRzzbG+wMxt\nbm6e6neIiFw02Wy27841Z/Kk/fzzz/t6NBolkUgwMTHBxMRE37qIiPR3Jk/a0WiU4HsSiQSFQsGX\n2dnZQ0ezQ/+4iIg8UAbtEXkmT9of+9jHfN05R6vVYm9vj1arRbPZZGlpiddff91fS6VSTE5Okk6n\nmZyc7Kmn02ni8fhZNFtEZOScSWi///3v9/Xd3V1KpVJP2djYoFwus7GxQalUIpVKMTU1RSaT6Vsm\nJiaIRqPEYrGe0n1NROQiOvN0i8Vi5HI50uk0Dz30EM1m81CpVCrs7Oyws7NDpVKhWCxy8+ZNfz0a\njZLL5cjlcmSz2Z5jLpfTbBURubDOPLQjkQjxePyuXRzb29tsbW35cvC81WoxPj5Oq9Via2uLer1O\nsVhkfHzcl0Qi0bfE43ESiQTR6KC9ikVERtdI9iOk02lSqRTz8/M453DO0W63fb1arbK+vs76+jrF\nYpH19XVu3brlr21tbTEzM0OhUPDHoATnCm0RCaORDG0zu+sMkqDvOpPJMD8/T7Va9aVWq1Gr1Wi3\n2z3lzp073Lp1y58nEgn/j0M6ne4pqVRKUw9FZCSNZGgPE3SBDOq7brValEolisUixWKxb318fJxs\nNks2m2VqaurQMZ1OE41G71o0NVFEztqRQ9vMIsB3gLecc79kZjngeeAx4A3gSedc+VRaeY8ikQiZ\nTIZEIsHs7Cy7u7u+NJtNdnd3qdVqVCoV/4S+tbXF8vKyP2+32z7Uu8M9GOjMZrMKbRE5c/fypP1J\n4PtApnP+DPCic+7zZvY08OnOtXNnZkMHOyuVCuVyma2tLcrlck89FotRr9cBqNVqvh+9VCpx584d\nP5gZfMegoqmHInLSjpQqZnYZuAL8FvAbnctPAO/u1J8DvsWIhPZRTExMEI/HKRQKPX3fwaBno9Hw\nXSkbGxsUi0WWlpZ65pfncjlmZmZ8yefzFAoF8vk8MzMzCm0ROXFHTZXfAf490N2JPOecWwFwzi2b\n2exJN+40RSIRIpHBK9Pu7e0Rj8fJZDLMzs5SrVZ9d0qlUqFSqfhX84OgX11dZXV11V8fGxsjlUr5\nkk6nSSaTfrAzmUyeyc8qIhfH0NA2s18EVpxzL5vZe+7y0dNdxOSMBbNTMpnMwM+USqWeaYfdx2Kx\nCOBf+OlXstms/8cjEokQjUYP1dVvLiLdjvKk/S7gl8zsCjABTJrZHwLLZjbnnFsxs3lgddAfcPXq\nVV9fWFhgYWHhmM0eDalUirGxMaanp3nsscdoNBo0Gg12d3d9vVarUa/X/fHGjRv+vNlsMjU1NbBk\ns1l1sYg8IBYXF1lcXBz6uXta5c/M3g18qjN75PNA0Tn3bGcgMuecO9SnbWbu2rVr99D0i6NWq/lB\nzs3NzUP17e1tJiYmSCaTfnnagyUY1Ox+27P7fGxs7Lx/TBE5Baexyt9vAy+Y2VPATeDJY/xZF1I8\nHmdmZoZcLke73abVavljsNLh5uYmGxsbbGxssLm5yfLysq9vbGyQTCbJ5/OHyszMDNPT0wptkQfM\nmayn/aA+aQ/Tbrf9IljB4ObBervd9m+I9iuxWIxkMukHNvvV1S8uEj7nup629BeJRPx64YOUy2XW\n1tb8uipBWVtbo1gsUq/XmZ6e7nkKn56e9k/i7Xbbv70ZDGwGx+66iISDQnvETUxMMDc3x9TUFJcv\nX6bRaFCv1/1AZ71eZ3d31x8bjQZvvfUWr732mv9MOp3265P3W6c8kUic948pIkek0B5xwYDjoKfx\n3d1dNjc3fQk2lNjb22N3d5dyuexDfWdnh3K57LtNgm6URCLhBzUHHUVkNCi0Qy7YVCKTyfDwww/7\nLdu6t3TrflV/c3OTUqnEjRs3/Kv7kUjEd7FMT0/7Epzncrnz/jFFpEOhHXLByzh3exoOdgHqLtvb\n276+t7fnX+hptVoUi0U2Nze5efMmkUiEWCzWMzWx31Hrk4ucDYX2AyBYJ3yQnZ0d1tfXWV1d9cc7\nd+74883NzZ5NJYISnJsZ4+PjPeugD6qLyPEotIV4PE4+nyeVSvHQQw9Rr9d73uIM3vIMlrVtNpus\nrKywtLTkr01MTDA5OUkmk+l71DorIidDoS2MjY0xNjY2cLCz1Wr5vvDul4Hq9TpbW1tsbGwQi8WY\nnJw8tAtQUJLJpP+eWCzm693nehoXGU6hLUMFm0okk0nm5+dpNpvs7e31lEql4jdhDjZivnXrlj/f\n29sjl8v5gc3uY1BXaIsMp9CWoczMPxUPUq1WfUAHZWdnh62tLXZ2dqjX6/6JutlsUiwWKZfL3Lp1\ni1gsRiwWI5FI9Ky70n2eSCQ09VAEhbackCBcZ2f7L6veaDRYW1tjdXXVH5eXl/0a5Gtra0xPTzM7\nO0uhUOg5BnWFtohCW07IsK6NsbExcrkciUSCQqHA448/Tq1W6ynd3S3B5syrq6v+2tjYGOl02r/6\n312GzZARuSgU2nImotHo0GDd2NigVCr5Y6lUYmdnxw98tttt/+r9wRkqQT3oaolGo4yNjRGNRv21\nWCx2192KRMJAoS0jI51O++Vsm82mL3t7ezSbTWq1mn8hKFgFcX193V+r1+tks9mBuwTlcjmFtoSe\nQltGxrDBznq97memBKX7vFKpMD4+TjQapVar+Vf479y549dRicfjfnOJRCJxqIyPj5/hTyxy7xTa\nEhrBjj3T09PA/obK3aXVah1awnZ5ebnnPJlM+jc5C4WCr8/OzjIzM6PQlpGn0JbQGPY6fLvdJp/P\nMzExQT6f59FHH6VarVKtVqnValSrVT/I2W63abfb/qWh69ev02q1iEQifrAzlUr1HIM+ec0nl/Ok\n0JYLIxKJ+CVnB9na2qJYLPqBzqAeHOv1ut9UufvYXQ8GOoNBzqAeFPWby2lSaMsDJdhUIpfL8cgj\nj/i1U4I1Ver1un86r1arVCoVlpaWuH79ur+WyWTIZrMDizaVkNOk0JYHSjDYOWgBq2az6dcf73cM\nphLu7e2xvb3tP7+yskI8HvcDnEG933k8Hj/jn1ouEoW2SJdYLOa7QdrtNs453/8d1Dc2Nnq6VYIS\ndLlEIhHy+TyFQsEfu5e0VWjLcSi0RbqY2dANHWKxGMlkknw+z6VLl3w3StB9sru7i3MO2J/hEkxJ\nfO2113DOYWa+7z2VSpFOpw+da1MJGUShLXKPksnkXdcHr1QqPU/g6+vr/lgqlSiXy2SzWaanpwe+\nDBTMNw92Jgrq3dfkwaTQFjlhwVudk5OTXLp0yW8i0X0MNpkIyvLyMq+//rq/PjEx4WetHJzFMjU1\nddcZMnKxKbRFTliwzsmgp/F2u+03WQ42XA7OYX9FxGazSaVSod1uU6vVKJfLh5asHR8f92959jtq\nPvnFdKTQNrMp4HeBvwO0gaeAV4HngceAN4AnnXPl02mmyMVhZmQyGb+9W/CyT/dLP9vb236hrGC3\noDfffNOf7+7uks/nmZmZIZ/P+9J9rtC+mI76pP0F4Jpz7l+aWQxIAZ8BXnTOfd7MngY+DTxzSu0U\nuTCCwc67DTYGb2Hm83m/OFalUvGlXq/7N0TNzC+a9eabb/prQd97MpkklUr1HJPJpF7ZD6mhoW1m\nGeAfO+c+CuCc2wPKZvYE8O7Ox54DvoVCW+REBF0hMzMzfX+90WgcWmdlfX2dtbU1isUia2trZLNZ\n8vk809PTh47T09Mkk0nMzA9sBvXuazJ6jvKk/Tiwbma/D/w94DvAvwXmnHMrAM65ZTPrv2WJiJy4\nWCxGLpcjmUwyNzdHvV6n0Wj0lO5r9XqdYrHI7du3/bVYLObXIp+amjp0HLTRs5yvo4R2DHgH8HHn\n3HfM7HfYf6J2Bz538Ny7evWqry8sLLCwsHAfTRWRQDQaHTr1sFwus7Gx4Y+bm5tsbm76V/Xb7Ta7\nu7vUajW2t7cHdqGMj4/7pW0PHvU0fnIWFxdZXFwc+jkLXgIY+AGzOeB/O+d+snP+c+yH9t8A3uOc\nWzGzeeDPnHM/3ef3u2vXrt3HjyAix9FqtfyqhkHpPq9Wqz2zV7pnsZTLZSqVSt+ule6j9u08PVeu\nXME5d2g0eeiTdieUl8zsbzrnXgXeC3yvUz4KPAt8BPjqyTZZRI5j2GBno9Egl8v5QcyDpVqt9rzQ\ns729TbVa5fbt2/56PB73T+UTExOH6npl/+QddfbIJ4AvmdkYcAP4VSAKvGBmTwE3gSdPp4kichqC\nxavy+XzfX2+1WqytrbG2tsb6+jqrq6t+sDMowWDpoJLJZPzUw+7ZLt1ro2tq4r0Z2j1y7C9Q94hI\nKLXbbRqNBrVajXq93rd0L2vbfQzqZnZoE+buzZi7Q1163Xf3iIg8mCKRiJ96OEjwElDwAlBQ6vW6\nH/QMdvwZVILlcg+WWCzG2NiYFs86QKEtIvctmGGSz+fZ29s7VBqNxqFNmJeWlnrOp6am/Nzx6elp\ncrlcz/Fu/2g8iBTaInLfgsHOQQOOzWaTnZ0dtre3fTl4HmzbZmZsbW1RrVZZWVnxa7gEG0kET/0H\ny4O2U5BCW0ROzdjYmF9udpBgkHN1dZW1tbWe+traGgCzs7PMzs5SKBQOHY8yQ+Ui9ZsrtEXkXKVS\nKebn58lmszzyyCN+udpgmdpGo9HT5VKr1bhx4wavvvoqe3t7tNttJicne0o6ne45j8UuTtRdnJ9E\nREJp2GBntVplY2ODUqnkj6VSia2tLX8MZqIMmqkyPj7uu1uCAc7u8zANdiq0RWSkBXPJp6amuHz5\nMs1mk2azyd7enp9i2P1S0Pb2NsvLy/zoRz/y1yYmJsjlcgN3Ckqn0+f9Yx6ZQltERtqwNzvb7XbP\nbJTg2F03M8bGxvxna7Ua6+vrfg2V8fFxEonEoRIMdCYSiZHpF1doi0iomZlfnfDSpUs9myoHZWtr\ny7/Z2a80Gg1mZmaYnZ1lZmaGQqHQU+LxuEJbROQkHOV1+EwmQywWY2pqikuXLlGtVqnVav5Yr9d7\ndg/a3d1laWmJN954w19PpVKHXgwKBj3T6fSZbSqh0BaRCy9YZ2XQ1MNGo+EHOEulEsVi8dBxcnKy\nZ6PlTCbjN1zOZDIkk0nflRPMPT9YPwkKbRF54AXzydPpNPPz84fWUgk2Wq5Wq37Lt1KpxNLSkr8e\njUb9QOfU1JQf+AzK1NTUibRVoS0iD7xIJOIHHAcJBjfL5XLf+t7eHmZGvV7HOUetVqNYLPqn/KAk\nEomB9aNsKqHQFhE5gnQ6TSqVYm5uDucc7Xabdrvt69Vq1XelBOXWrVu+i2Vra6tn2dpCoUA+n6dQ\nKPhrR+kXV2iLiBxB8BQ8qG86eGLOZDLMzc1RrVZ9d0ow4Nk9o6XZbHLnzh1u377dc30YhbaIyAkI\n5ntns9m+v95qtXqewtfX1ymVSqyvr/trzWZz6PcotEVEzkAkEiGTyZBIJCgUCjQaDb9hRFBvt9v+\n85/73Of6/jkKbRGRM2BmQwc7j2L4UKWIiIwMhbaISIgotEVEQkShLSISIgptEZEQUWiLiITIkULb\nzP6dmf2VmS2a2ZfMbNzMcmb2DTP7oZl93cxOZjUUEREZaGhom9kl4NeBdzjnFtif2/3LwDPAi865\ntwHfBD59mg0VEZGjd49EgZSZxYAJ4BbwBPBc59efAz5w8s0TEZFuQ0PbOXcb+A/Am+yHddk59yIw\n55xb6XxmGZg9zYaKiMjRukey7D9VPwZcYv+J+1eAg8tRDV+eSkREjuUoa4/8PHDDOVcCMLM/Bv4R\nsGJmc865FTObB1YH/QFXr1719YWFBRYWFo7XahGRC2ZxcZHFxcWhn7Nh67ea2TuBLwI/CzSA3wde\nAh4FSs65Z83saSDnnHumz+93165du/efQETkAXblyhWcc4d2Kx76pO2c+7aZfQX4LtDsHP8zMAm8\nYGZPATeBJ0+2ySIictCRlmZ1zv0m8JsHLpfY7zoREZEzojciRURCRKEtIhIiCm0RkRBRaIuIhIhC\nW0QkRBTaIiIhotAWEQkRhbaISIgotEVEQkShLSISIgptEZEQUWiLiISIQltEJEQU2iIiIaLQFhEJ\nEYW2iEiIKLRFREJEoS0iEiIKbRGREFFoi4iEiEJbRCREFNoiIiGi0BYRCZEzD+3FxcWz/sr7onae\nLLXz5IShjaB2nhaF9gBq58lSO09OGNoIaudpUfeIiEiIKLRFRELEnHOn+wVmp/sFIiIXlHPODl47\n9dAWEZGTo+4REZEQUWiLiITImYa2mb3PzK6b2atm9vRZfve9MLM3zOz/mdl3zezb592egJl90cxW\nzGyx61rOzL5hZj80s6+b2dQItvGzZvaWmf3fTnnfebax06bLZvZNM/uemb1iZp/oXB+1+3mwnb/e\nuT5S99TM4mb2l52/M6+Y2Wc710fmft6ljSN1L4c5sz5tM4sArwLvBW4DLwEfdM5dP5MG3AMzuwH8\nA+fcxnm3pZuZ/RywA/yBc26hc+1ZoOic+3znH8Kcc+6ZEWvjZ4Ft59x/PK92HWRm88C8c+5lM0sD\n/wd4AvhVRut+Dmrnv2L07mnSOVc1syjw58AngH/BaN3Pfm38Z4zYvbybs3zSfifw1865m865JvBH\n7P+fbxQZI9h15Jz7X8DBf0ieAJ7r1J8DPnCmjTpgQBth/56ODOfcsnPu5U59B/gBcJnRu5/92vlw\n55dH7Z5WO9U4EAMco3c/+7URRuxe3s1ZBtPDwFLX+Vv8+P98o8YBf2pmL5nZvz7vxgwx65xbgf2/\n4MDsObdnkF8zs5fN7HfPu8vhIDP7CeDtwF8Ac6N6P7va+ZedSyN1T80sYmbfBZaBP3XOvcSI3c8B\nbYQRu5d3M3JPkyPiXc65dwBXgI93/pM/LEZxDud/An7SOfd29v+yjMx/hna6HL4CfLLzJHvw/o3E\n/ezTzpG7p865tnPu77P/XyzvNLOfYcTuZ582/m1G8F7ezVmG9i3g0a7zy51rI8c5d6dzXAP+mP2u\nnVG1YmZz4Ps/V8+5PYc459bcjwdP/gvws+fZnoCZxdgPwj90zn21c3nk7me/do7qPQVwzm0B3wLe\nxwjeT+ht4yjfy37OMrRfAn7KzB4zs3Hgg8DXzvD7j8TMkp2nGswsBfxT4K/Ot1U9jN7+t68BH+3U\nPwJ89eBvOAc9bez8ZQ38c0bnfv4e8H3n3Be6ro3i/TzUzlG7p2Y2E3QrmNkE8Avs97+PzP0c0Mbr\no3YvhznTNyI7U2m+wP4/Fl90zv32mX35EZnZ4+w/XTv2Byq+NCrtNLMvA+8B8sAK8FngvwH/FXgE\nuAk86ZzbHLE2/hP2+2LbwBvAvwn6Oc+Lmb0L+J/AK+z/b+2AzwDfBl5gdO7noHZ+iBG6p2b2d9kf\naIx0yvPOud8ys2lG5H7epY1/wAjdy2H0GruISIhoIFJEJEQU2iIiIaLQFhEJEYW2iEiIKLRFREJE\noS0iEiIKbRGREFFoi4iEyP8Hut1Zi7YrJlsAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "rc = (imp[1:,:] - imp[:-1,:]) / (imp[1:,:] + imp[:-1,:])\n", "\n", "plt.imshow(rc, cmap='Greys', aspect=0.2)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll use a wavelet function from [`bruges`](https://github.com/agile-geoscience/bruges). This is not cheating! Well, I don't think it is... we could use `scipy.signal.ricker` but I can't figure out how to convert frequency into the 'width' parameter that function wants. Using the Ricker from `bruges` keeps things a bit simpler." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import bruges\n", "\n", "w = bruges.filters.ricker(duration=0.100, dt=0.001, f=40)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's make sure it looks OK:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD7CAYAAACIYvgKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGtRJREFUeJzt3XmQFPXdx/H3d4FFOUQMCsqyCBJATQhBRYyPcT3igfGo\n0hiIR6QSD6IRE5JK1KTAoyo+poxHQEXFMypRTBQjiucWUYyP4oOonILsLtfqI/eh7PF7/vjNuMuy\n93RP90x/XlVT7M72dv+mme3P/M425xwiIpJsBVEXQEREoqcwEBERhYGIiCgMREQEhYGIiKAwEBER\noGPUBWjIzDTWVUSkjZxzlsnvx7Jm4JzTwzkmTZoUeRni8NB50LnQuWj+EYRYhoGIiGSXwkBERBQG\ncVZSUhJ1EWJB56GOzkUdnYtgWVDtTUExMxe3MomIxJmZ4fKxA1lERLJLYSAiIgoDERFRGIiICAoD\nERFBYSAiIigMREQEhYGIiKAwEBERFAYiIoLCQEREUBiIiAgKAxERQWEgIiIEFAZmNt3MKs1sYTPb\n3GVmy81sgZkND+K4IiISjKBqBg8Bpzb1QzM7HTjEOfdN4HLg3oCOKyIiAQgkDJxzbwIbm9nkbODR\n1LbvAD3MrHcQxxYRkcxlq8+gL1BR7/s1qedEckJNDdTWRl0KkfB0jLoAjZk8efLXX5eUlOhepxK5\niy6C7t1h2rSoSyICpaWllJaWBrrPwO6BbGb9geedc8Ma+dm9wBvOub+nvl8CHO+cq2xkW90DWWJl\n7ly48EKoqoIXXoARI6Iukcju4nYPZEs9GjMLuBjAzEYBmxoLApG4qamBCRPgz3+Gm27yX+uziuSj\noIaWPgHMAwabWbmZjTOzy83sMgDn3GzgUzP7BJgG/CKI44qE7cEHoVs3OP98GDcOtm2Dp56KulQi\nwQusmSgoaiaSuNi8GYYMgdmz65qG5s71/QeLF0OXLtGWTyQtbs1EInnlppvgzDN37yP4/vdh1Ci4\n9dboyiUSBtUMRBpRWwv77Qcffwx9GwyCXrIETjoJ1qyJpmwiDalmIBKSTz6Bnj33DALwTUdffgnr\n1mW/XCJhURiINOK99+DIIxv/mZn/2fz52S2TSJgUBiKNaC4MwP/svfeyVx6RsCkMRBqhMJCkUQey\nSAM1NbDvvlBe7vsNGlNeDiNH+n4Dy6jbTiRz6kAWCcHSpdC7d9NBANCvnx9xpBFFki8UBiINtNRE\nBHWdyGoqknyhMBBpoDVhAAoDyS8KA5EGFAaSROpAFqmnuhp69PAdw/vs0/y2a9fCsGHw+efqRJZo\nqQNZJGCLFvnO4ZaCAOCgg6CwEMrKwi+XSNgUBiL1tLaJKE1NRZIvFAYi9SgMJKkUBiL1KAwkqdSB\nLJKya5efefzZZ/7uZq1RWQlDh8KGDepEluioA1kkQB9/DAMGtD4IwM9U7tYNVq4Mr1wi2aAwEElZ\ntsx/ym+roUP974rkMoWBSMrKlXDIIW3/vUMOUc1Acp/CQCRlxYr2h8GKFcGXRySbFAYiKStWwMCB\nbf+9gQMVBpL7FAYiKWomkiTT0FIR4Kuv/BIU27dDx45t+90tW+DAA2HbNg0vlWhoaKlIQFatgqKi\ntgcB+BDp0gXWrw+8WCJZozAQof1NRGlqKpJcpzAQof0jidI0okhyncJABIWBiMJAhPYPK03T8FLJ\ndQoDEdRnIKIwkMRzzl/IM6kZqJlIcp3CQBJv3Tq/8mj37u3fR58+sHWrf4jkIoWBJF6mTUQABQV+\n+etPPw2mTCLZpjCQxMt0JFGamooklykMJPEUBiIKA5GMO4/TBg7UiCLJXQoDSTzVDEQUBiIKAxG0\nhLUk3Nat/qb227dnvvx0Jstgi2RCS1iLZCjdXxDEfQg6d/bBUlGR+b5Esk1hIIkWVBNRmpqKJFcp\nDCTRFAYinsJAEu3TT/3M4aBoeKnkKoWBJFp5ORQXB7e/fv3UZyC5SWEgiVZREWwYFBcrDCQ3KQwk\n0Soq/Kf5oKhmILlK8wwksbZvh/33D2aOQVp6rsGOHdChQzD7FGmJ5hmIZKCiAoqKggsC8HMNevaE\nysrg9imSDQoDSazy8mCbiNL69fP7FsklCgNJrKA7j9PUiSy5SGEgiRV053GaOpElFykMJLEUBiJ1\nFAaSWAoDkTqBhIGZnWZmS8xsmZn9rpGfH29mm8zs/dTjD0EcVyQT6kAWqZPxqutmVgBMAU4C1gLv\nmtlzzrklDTad65w7K9PjiQTBufBqBupAllwURM1gJLDcOVfmnKsCZgBnN7JdgKO5RTKzcSN06uQn\niAWtTx/YsAF27Qp+3yJhCSIM+gL1PwetTj3X0DFmtsDMXjCzwwI4rki7hVUrAD/zuE8fWLMmnP2L\nhCFbN+ebDxQ753aY2enAs8DgpjaePHny11+XlJRQUlISdvkkYcIMA6jrRA5yeWyRtNLSUkpLSwPd\nZ8ZrE5nZKGCyc+601Pe/B5xz7r+b+Z1PgSOccxsa+ZnWJpLQ3X03fPABTJsWzv7HjoUzzoALLwxn\n/yL1xWVtoneBQWbW38wKgTHArPobmFnvel+PxIfQHkEgki1hzT5OUyey5JqMm4mcczVmdhXwMj5c\npjvnFpvZ5f7H7j7gPDMbD1QBO4EfZ3pckUxUVMApp4S3/379YNGi8PYvErRA+gyccy8BQxo8N63e\n11OBqUEcSyQI2egzmDMnvP2LBE0zkCWRstWBLJIrFAaSOLW1fthnUVF4x9AsZMk1CgNJnMpK2Hdf\n2Guv8I7Rqxfs3OnvoiaSCxQGkjhhNxGBv3uamooklygMJHGyEQagMJDcojCQxFEYiOxJYSCJE9bS\n1Q2pE1lyicJAEifs2cdpmoUsuURhIImjZiKRPSkMJHEUBiJ7ynjV0qBp1VIJU1UVdO0KO3ZAx5AX\ncN+yBQ46CLZu9UNNRcISl1VLRXLGunVwwAHhBwH4u6gVFMDmzeEfSyRTCgNJlGw1EaWpqUhyhcJA\nEmX16nDXJGqoqMgfUyTuFAaSKKoZiDROYSCJopqBSOMUBpIoqhmINE5hIImimoFI4xQGkiiqGYg0\nTpPOJDGyOeEsTRPPJBs06UykDdauhd69sxcE4CeedegAmzZl75gi7aEwkMTIdn9BmvoNJBcoDCQx\nst1fkKZ+A8kFCgNJDNUMRJqmMJDEUM1ApGkKA0kM1QxEmqYwkMRQzUCkaQoDSQzVDESapklnkgi7\ndkG3brBzpx/3n01bt0KfPrBtmyaeSTg06Uykldat8xPOsh0EAN27+4lumngmcaYwkESIqr8gTf0G\nEncKA0mEqPoL0tRvIHGnMJBEUM1ApHkKA0mE1aujDwPVDCTOFAaSCBUV0TcTqWYgcaYwkERQzUCk\neQoDSQTVDESap0lnkveinHCWpolnEiZNOhNphbVr/YU4qiAAP/GsUyfYuDG6Mog0R2EgeS/qJqI0\nDS+VOFMYSN4rL4f+/aMuBRQX+7KIxJHCQPLeqlXxCIP+/X1ZROJIYSB5r6wMDj446lL4MpSVRV0K\nkcYpDCTvqWYg0jKFgeQ91QxEWqYwkLxWW+s7bYuLoy6JagYSbwoDyWuffebH+HftGnVJ/M11tm2D\n7dujLonInhQGktfi0l8AfuZxcbGaiiSeFAaS1+LSX5CmfgOJK4WB5LU41QxA/QYSXwoDyWuqGYi0\njsJA8ppqBiKtozCQvKaagUjrBBIGZnaamS0xs2Vm9rsmtrnLzJab2QIzGx7EcUWa45xqBiKtlXEY\nmFkBMAU4FTgcGGtmQxtsczpwiHPum8DlwL2ZHlekJV98AZ07wz77RF2SOgceCBs2wJdfRl0Skd0F\nUTMYCSx3zpU556qAGcDZDbY5G3gUwDn3DtDDzHoHcGyRJsWtVgD+BjtFRVrKWuKnYwD76AvUv2XH\nanxANLfNmtRzlY3tcO3aAEqVo5q6JWL6+fr/ph8FBf4iU1AAHTv6O2p16KDbK8atvyAt3W8weHDU\nJYmWc1BdXfeorYWaGv+vc3WP9Lb1/21sX0nTrVuwtd4gwiBwQ4ZM/vrrwsISOncuiaws2dTSG73+\nv/UftbV1f0jpP6yaGh8Ke+8NXbr4f3v0gJ49/aN3bz8btrgYBgyAYcPisWRDkOJYM4D87TfYtg0+\n+AA+/dTXfMrL/XIgGzfCpk2weTPs2OHvRb1zJ1RV+Q8t6Q8v6Q80BQW7f9hp7INQY5L04eerr0oZ\nPryU444Lbp9BhMEaoP4yYEWp5xpu06+Fbb62devkAIqVbLW1/o8t/Ye3Y4f/Y9y40T/Wr/d/rB98\nAMuXw6JFPhSOPBJGj/aP7t2jfhWZKSuDgQOjLsWe8mVE0ebN8PzzMGcOvPeefz8dfjgMGuQ/ZAwb\n5j90pD+A9OjhP5SkH506+Qu/tEdJ6uHdcMMNGe8xiDB4FxhkZv2BdcAYYGyDbWYBVwJ/N7NRwCbn\nXKNNRBKMggLfedq5M+y7b8vb79oFH38M77wDjzwCl14KJ54I48bBWWfl5qeuVavghBOiLsWe+veH\nV16JuhTtU1sLzzwDDz8Mb74Jxx/v3x+//S0ceqi/wEtuyjgMnHM1ZnYV8DK+Q3q6c26xmV3uf+zu\nc87NNrPRZvYJsB0Yl+lxJViFhfDd7/rHFVf4av2sWTBpEtx4o3+MHp1boRD3PoNc4hz885/+/bD3\n3nDNNfDkk/EaqSWZMReznhczc3ErU5LV1tZdBHr2hMcei+cFtjE9evj26/32i7oku1u1Co47Dioq\nWtw0FpYvhwsu8H1RN94IZ5yRWx8KksDMcM5l9L+iFjtpVkEBnHsuLFwI55wDRx8NL7wQdalatmmT\n/zTbs2fUJdlTUZHvWN21K+qStGzmTDj2WN9cOH8+/PCHCoJ8FcvRRBI/BQUwcSKMGgVjxsCFF8LN\nN/sRIHGUHkkUxwtXx47Qpw+sXh3PDm7wgw9+8xvfQfzii3DEEVGXSMKmmoG0ybHHwvvvw7x5MH58\nfMd3x7W/IC3O/Qa1tXDxxbB0qa8NKAiSQWEgbbb//vCvf8GCBXD99VGXpnFxnWOQFte5Bs7B1VfD\nunW+ryiOzWwSDjUTSbt07w6zZ/uO0F694Ne/jrpEu8uFmkEcw+DGG+Gtt6C01I8akuRQzUDarVcv\nePlluOsueOKJqEuzu08+iW97PPiyrVgRdSl2N20a/O1v8NJLfiSWJIvCQDLSr5/vZJwwIV6fdJcs\ngaFDW94uKkOH+jLGxeLFvslv9mw/a1iSR/MMJBC33upHnbz2WvRLDOza5SdDbd7sZ2DH0caNfsmG\nLVuiH/FUXQ3f+54fPjp+fLRlkfbRPAOJjYkT/UV4ypSoS+KbX/r1i28QgO+Y7dIlHiv03nKLX7Lk\niiuiLolESR3IEogOHfyaRsccA6eeCkOGRFeWpUujPX5rDRniy9q3b3RlWLDA9/m8/370NRSJlmoG\nEphBg+CGG+CSS/xY9ajEvb8gLep+g+pqP5/gttv8rGhJNoWBBOqKK/xFZubM6MqQazWDqDz8MHzj\nG342uYjCQAJVUAB/+pMfmVJVFU0ZVDNo2c6dMHmy7y9Q85CAwkBCcPLJflLV9OnZP7Zzqhm0xl//\nCiNH+oUHRUBDSyUk8+fDmWf65Y+zeTvNzz7zn7i/+CL+n3irq/1M7i++8COLsmXjRn//5blz/Q1p\nJPdpaKnE1hFH+KUq7roru8dN1wriHgTgVy8dONAHZjbdeiucfbaCQHanMJDQ3HQT/OUvsGFD9o6Z\nK/0FadnuN1i7Fu67z/cXiNSnMJDQDB7s7487dWr2jpkr/QVp2e43uP12uOgiDSWVPSkMJFQTJ/ow\n+PLL7BxPNYOmbdkCDz4Iv/pVdo4nuUVhIKE67DAYMQIefzw7x1PNoGnTp8MPfhDv+zxIdDSaSEL3\n2mvwy1/CRx+Fu4jdV1/5pZe3bIHCwvCOE6RNm3yTzdat4XZ6V1fDIYf4yYBHHRXecSQaGk0kOeHE\nE/3F+aWXwj3OihV+JdBcCQLwC8R16wZr1oR7nJkzfY1AQSBNURhI6Mx838Ftt4V7nFzrL0gLu9/A\nOX/uJ04M7xiS+xQGkhU//rFvG1+wILxj5Fp/QVrY/Qb//re/t8OZZ4Z3DMl9CgPJisJC329w++3h\nHUM1g8bdcQdcc030Nx2SeNPbQ7LmZz+D557zyy+EQTWDPa1e7W9uf/HF4exf8ofCQLKmVy/fVPHw\nw8Hv2zn/6ToXw2DoUFi0KJx9338/jB3rO6lFmqMwkKwaPx7uvTf4m98sW+ZH5uy/f7D7zYaDD/ZL\nSgd9C8yqKnjgAd3XWFpHYSBZdcwxfoXO114Ldr/z5vmbuueiggJ/Xt5+O9j9zprlF8L71reC3a/k\nJ4WBZJWZ/6R6zz3B7jeXwwB82efNC3af99yjWoG0nsJAsu6CC3yn5urVwe1TYbC7Zcvgww/h3HOD\n26fkN4WBZF337jBmjO/cDMLGjVBeDsOGBbO/KBx1FCxcGNyCfvfeC+PGQefOwexP8p/CQCIxfrzv\n3AziPsn/+Y+/mHbsmPm+otK1q7/ZzPz5me9r50549FG4/PLM9yXJoTCQSHz7237htGefzXxfud5E\nlBZUU9GTT8KoUTBgQOb7kuRQGEhkrroqmBvfKAzqOAdTpvhzK9IWWsJaIlNV5VfSnDPH1xTao7oa\nevaEsjLYb79gy5dt5eW+uWv9+vYvZ/3223628dKlWn4iSbSEteS0Tp3gssvg7rvbv48PP4R+/XI/\nCMC/jk6dYOXK9u9j6lTfH6MgkLbSW0YiddllMGOGX1WzPfKliQh8bSCTpqLKSnjhBT+KSKStFAYS\nqYMOglNOgUcead/v51MYQGZh8MADcN55vtlMpK0UBhK5dEdye9YrUhh41dV+bsGVVwZfJkkGhYFE\n7r/+C/bay3ckt8Xatf5+x4MHh1OuKAwf7m/fuWVL237vH//wnfHDh4dTLsl/CgOJnBlcdx1MmuSH\nRrbWK6/4IMmnztLCQhg5sm0L+dXUwA03wPXXh1cuyX959GckuexHP/JLMTz/fOu2T4+nv/TScMsV\nhZ//3L+21poxwy/ffdpp4ZVJ8p/CQGKhoABuugn++MfW9R3MmwebNsHo0eGXLdvOO8/fqGfhwpa3\nraqCyZPh5pvbPzdBBBQGEiNnneWbSWbObHnbO+6Aq6/OryaitMJC+MUv4M47W9720UehuBhOOCH8\nckl+0wxkiZU5c2DCBPjoo6YXnisrgxEjYNUqvwJqPvr8c98xvmxZ03dv++orv82TT+bXiCppO81A\nlrxzyin+4vf4401vM3UqXHJJ/gYB+HNw3nkwbVrT29x/v7+LmYJAgqCagcTOO+/AmWf6FU0bXui2\nbfP3DH733fxflfPDD+HUU30NqLBw95+9+ir85Cf+31y+j4MEQzUDyUtHHw2PPQbnnOODIc05v47R\n97+f/0EAfvG+ww6DBx/c/fk33oCxY+GZZxQEEhzVDCS20uvsPP20bzufMsXfuOWpp5Izuertt30N\noFcvP7u4b1///dNPQ0lJ1KWTuAiiZqAwkFibNctf/E46yV8MTz45P0cQNaemBl580feVvPWWn218\n8slRl0riRGEgieCcxtCn6VxIY9RnIImgi18dnQsJS0a3EDeznsDfgf7AKuB859weK9Ob2SpgM1AL\nVDnnRmZyXBERCVamNYPfA68654YArwPXNrFdLVDinPuugkBEJH4yDYOzgfRtSR4BzmliOwvgWCIi\nEpJML9AHOOcqAZxz64EDmtjOAa+Y2btmlofrTIqI5LYW+wzM7BWgd/2n8Bf3PzSyeVPDgI51zq0z\ns/3xobDYOfdmU8ecPHny11+XlJRQogHVIiJfKy0tpbS0NNB9ZjS01MwW4/sCKs2sD/CGc+7QFn5n\nErDVOfeXJn6uoaUiIm0Qh6Gls4BLUl//FHiu4QZm1sXMuqW+7gqcAnyU4XFFRCRAmdYM9gOeAvoB\nZfihpZvM7EDgfufcD81sAPBPfBNSR+Bx59wtzexTNQMRkTbQDGQREYlFM5GIiOQBhYGIiCgMRERE\nYSAiIigMYi3oSSW5Suehjs5FHZ2LYCkMYkxvdk/noY7ORR2di2ApDERERGEgIiIxnXQWdRlERHJN\n3s1AFhGR7FMzkYiIKAxERCRGYWBmp5nZEjNbZma/i7o82WRmRWb2upl9bGYfmtnVqed7mtnLZrbU\nzOaYWY+oy5otZlZgZu+b2azU94k8F2bWw8yeNrPFqffH0Qk+F78ys4/MbKGZPW5mhUk5F2Y23cwq\nzWxhveeafO1mdq2ZLU+9b05pzTFiEQZmVgBMAU4FDgfGmtnQaEuVVdXAr51zhwPHAFemXv/vgVed\nc0OA14FrIyxjtk0AFtX7Pqnn4k5gduqmUd8BlpDAc2FmBwG/BEY454bhl8MfS3LOxUP462N9jb52\nMzsMOB84FDgduNvMWuxcjkUYACOB5c65MudcFTADODviMmWNc269c25B6uttwGKgCH8OHklt9ghw\nTjQlzC4zKwJGAw/Uezpx58LM9gGOc849BOCcq3bObSaB5yKlA9DVzDoCewNrSMi5SN0meGODp5t6\n7WcBM1Lvl1XAcvw1tllxCYO+QEW971ennkscMzsYGA78B+jtnKsEHxjAAdGVLKtuB37L7vfUTuK5\nGAD8n5k9lGoyu8/MupDAc+GcWwvcBpTjQ2Czc+5VEngu6jmgidfe8Hq6hlZcT+MSBgKkbg86E5iQ\nqiE0HPeb9+OAzewMoDJVU2quapv35wLfFDICmOqcGwFsxzcNJPF9sS/+k3B/4CB8DeECEngumpHR\na49LGKwBiut9X5R6LjFSVd+ZwGPOufS9pCvNrHfq532Az6IqXxYdC5xlZiuBJ4ETzewxYH0Cz8Vq\noMI5917q+2fw4ZDE98XJwErn3AbnXA3+VrrfI5nnIq2p174GfyvitFZdT+MSBu8Cg8ysv5kVAmOA\nWRGXKdseBBY55+6s99ws4JLU1z8Fnmv4S/nGOXedc67YOTcQ/z543Tl3EfA8yTsXlUCFmQ1OPXUS\n8DEJfF/gm4dGmdleqc7Qk/ADDJJ0Lozda8tNvfZZwJjUaKsBwCDgf1rceVxmIJvZafiREwXAdOfc\nLREXKWvM7FhgLvAhvqrngOvw/4FP4VO+DDjfObcpqnJmm5kdD0x0zp1lZvuRwHNhZt/Bd6R3AlYC\n4/AdqUk8F5PwHxCqgP8Ffg50JwHnwsyeAEqAbwCVwCTgWeBpGnntZnYt8DP8uZrgnHu5xWPEJQxE\nRCQ6cWkmEhGRCCkMREREYSAiIgoDERFBYSAiIigMREQEhYGIiKAwEBER4P8Bukco0pEsDOIAAAAA\nSUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(w)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now one more application of `apply_along_axis`. We could use a loop to step over the traces, but the rule of thumb in Python is \"if you are using a loop, you're doing it wrong.\". So, we'll use `apply_along_axis`.\n", "\n", "It looks a bit more complicated this time, because we can't just pass a function like we did with `product` before. We want to pass in some more things, not just the trace that `apply_along_axis` is going to send it. So we use Python's 'unnamed function creator', `lambda` (in keeping with all things called `lambda`, it's a bad name that no-one can quite explain)." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAADHCAYAAADS8QIeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnVusbdd51//f3ic2adqGAPKxiN3TW6IecwuNalVKkQMJ\nNASpjgKculSNHaMqD20TKEJ2/BL1oVJtCarwwANtinyiotq1DnEeUJo62EVUtLFdu5jmGFfypbng\nE0PBEEVJzt578LDX2B77O99tjDnXbZ/vJy3NcZlrzjHHHus//+sbY65NpRQkSZIk28HOuhuQJEmS\nxEnRTpIk2SJStJMkSbaIFO0kSZItIkU7SZJki0jRTpIk2SImiTYRvYeIniWi54jorrkalSRJksjQ\n6DptItoB8ByAdwH4CoDHAdxWSnl2vuYlSZIkLVOc9s0A/riU8lIp5TKA3wBw6zzNSpIkSSROTXjv\nmwF8scl/CYdCfgwiykcukyRJBiilEC+bItph3ve+9x2lX3nlFXzoQx9axWknceHCBbz//e9fdzNc\nsp3zsg3t3IY2AtnOXi5evIiLFy8e5T/1qU+J+00R7S8D+K4mf8Oi7AraDrlw4cKEUyZJkpxMzp49\ni7Nnzx7lNdGeEtN+HMD3E9EZIroGwG0APj3heEmSJInDsNMupewT0c8C+CwOxf8TpZSLztuO3Uk2\nmWznvGQ752Mb2ghkO5fF8JK/8AmIyvnz55d6jiRJkpPGBz7wAXEiMp+ITJIk2SJStJMkSbaIFO0k\nSZItIkU7SZJki0jRTpIk2SJStJMkSbaIFO0kSZItIkU7SZJki0jRTpIk2SJW8it/Lct+AjNJkuQk\nk047SZJki0jRTpIk2SJStJMkSbaIFO0kSZItYiUTkXNPPuZkZpIkVyvptJMkSbaIlS/5izDVSacT\nT5LkpJJOO0mSZIvY+IdrvP2jx0v3nSTJSWDlok1EqoBawtrznhToJElOKmuJaRMd/q/KVly50Hp5\nrzxanyRJsk2sRLSrSAMxEe0Rc++Yc4h2Cn+SJJvCRq4eqWjiPVXIR84/1zGTJEmmsHKnXSmlHMW3\nqxi26ZofTfNzRVm2a5/ahiRJrm7WJtrAlYLbirdUz9Nzive6Y+Mp1EmSRHBFm4huAHAewGkABwB+\npZTyr4joTQAeAHAGwIsAzpVSXg0cTxUorc4SaUvopfwyyuYihTtJEo+I094D8POllKeJ6NsBPElE\nnwXwQQCPlFLuI6K7AHwUwN3SATSnzamiVcU7KuIRAZ+6UmWO0EuEFO4kSSxc0S6lvAzg5UX6a0R0\nEcANAG4FcMtit/sBPIaJoi2c+9jyQEnELcHuEXV+PK+spz66z0nharrWJFk1XTFtIvpuAG8D8HsA\nTpdSLgGHwk5E12nv29k5/rR8FVTr1e7L13VHHHiPiPP3S3VW+ZQ15O21WKEj6xhJklw9hEV7ERp5\nCMBHFo6bq4WqHg899NBR+qabbsJNN90UbqDlrr2HdCKxcK9srny0zqLnHEmSbBcvvPACXnjhBXe/\nkGgT0SkcCvYnSykPL4ovEdHpUsolIroewFe19992221H6YjL1lw3a5NYZ7lwL9TS7tturXTPftp5\nPJZxA0iSZLM4c+YMzpw5c5R/9NFHxf2iTvvXAHyhlPLxpuzTAO4AcC+A2wE8LLwPwFh4ZFTcvbJ2\nbbjm4r2J0HZfno6K9lQhnirWKfZJsp1Elvy9A8BPAniGiJ7CYRjkHhyK9YNEdCeAlwCc046xu7t7\nLN8jugcHB8OC3SPoWloTcF4WEfre8rnq5yKFPknWT2T1yO8C2FWq3x05iSbaWnpUmKe4cunYkbbW\ndLv1ynheis1L9dpxen/bZRStnUmSrI6VPBEpiba3jQjpMl7a8a22RLaRtOTkeTkQi9svA+0mkiTJ\n6liJaPOYtoTmMjWBnFPIR98vta9XzEdEXsp75SP7pqNOks1jbaItuTXta77lSiOi3ivMbRw98j6+\nj5W32mxtvX7h9Ih6tCxJkvWzlh+MqnltK72nJeJAI2LJBXoOhy6VeXVeu61rjjjw6Hu1sl5S8JNk\neazt97SJSH3V+nZfXqbhuXK+7RXYudw831fLR7Zzpq2ykX2mkjeAJDnOSkS7/eC1AkxE2NnZObbV\nBFwTdI71lKRWFhF1ad8R0R7Zp6e93vV6+V4XPzcp0klisxLRPjg4OEpXgT44OMDOzg5KKcfEWBNv\ny41HwyotyxT1USdvHUM7f2TrlUXfa/VhhBTkJJnOSkR7f3//KF2Xq1Vx5iLV7ic58bpt95HS7ZYf\nt8ULDfQIuiXqvUJvHXMZ22id1VfLqE+S5DgrEe29vb2jdBXeuna7fmiJCAcHB8fEt31Pfd8UJz4a\nG+f5HqHzwhtRwZ5L1KN10euW3iOR4p0k87AS0f7Wt751lCYi7O7u4uDgAAcHB9jd3T0Sn1bIa15y\nq/U43I3zPN+vzbfH4USe/LPErc1b21WKdY+IS9cUFeyePuut790vSU4iKxHtb37zm0fpKtqnTp06\n2vJld1WsgeMf0OrGeSy81k1x42169Ik/SUxGxHRKmbQdbYd2DOuaoyK+Te48bxLJJrES0f7GN75x\nlG5F+9SpU9jb2zu2PXXq1FH4pL729/eP3Pnu7u6RMNdXFbD2IZ4q4G2+TUshk+jEZo+o97jVOUV2\nGQLO0z11vWU99ct67yadI0kqaxVtvq3p9iWVceFu89WlVzGv5dWhA8fj6BzuxPl+PZOcLaMCNqcI\n9+zrCXYkpDKS18qi9Lw3xTbZRtYm2tZLEnRJwFtH3or37u7usVCJ9LLi4S2WYM8ZWgF8Aet1516Z\nVD9FtCPpSL6nTGKO0EuSbCprE20usu02KtqS+5ZEXHLkrXi3Wy8eXq+BC3Ypx3972xLv0fAKz88V\nepHKpgr0XCGVHvFOsU6uBtYyESmJJhfXSJiEl0n7SELOxVtz5d4KlXo99VUFG8Cx9BQB1+p7xXzU\nSa9iX+96InmtjNf1hrCSZNNYiWhfvnz5WN5apteu45bi1ZqoW05cEnctHq6Jt5YeXTOulUn5XqQb\nh7RP++2gLeNp/p7efSW0GxwXV+0c3vH5uUbq5tg/SeZmbaLNtzxthS8099zmrYnOVsQtJ66JOI+b\nW1vrGrW+kLY9WMIiia0mnp4r5ftaZdEbgHYdc7tyqzxaH90nSeZk5Y+xt2gfroh7lQSVu+bohKcX\nN4+EWLwJTy/E0l43T7f7aPkoo+LXEwfviY+Pxtej6ei1amU99R4p8MkcrPxX/toPovZq0cRKC1u0\nZdaEZ69wW0sPoy5dip2PhFakOq2/In8Tr7xXnJe1HSmT8lHBXoaIL1u488Zw8ln572kT0dEj7KWU\no3Sbj3x4JVfa5j0HrDnpXjeupTVhj65ikZz5aMy8l4hwa3+faFmkrue8ve/RrnUuQZ9zn6mkkJ8s\n1vLvxnics354uYBLr/ZnXut7tXPysIQ22SkJaU94Jbo0UQu9aCGeiKBLq1giIZZeQfdCESPb3ro5\nRH9Z4ZWeby499VP3T04eK/9v7PzDwYVFE/CIoLcvS8CkUIrmzq3Qh+W8vYlQL3buxc3bOs2B1zrp\n+iMTnfzv0uIJm+d2PQGW8tE66fzLcuPezYAzh4incF/drES0r7nmmqN0KQX7+/tHgrO/v38kInt7\ne8c+BLVcEm7+4vXt+VpqXhO4Nt8TG+91472xc+01R9zcCq30OPIecfNctifYXpn2XumcVnt69x91\n8pG6SH1y8gmLNhHtAHgCwJdKKT9GRG8C8ACAMwBeBHCulPKq9F5JtPf397G3t3eFsLRCXled1A/c\nzs7OFYO/inQ9Zk23jr39ULcxc0ukInFxa1WJJdyeeHv78TLt5tKmtW8VPSGW3klPTcQ98fMEPCri\nVt46Z0/dlG2kz0bLl0neONZLj9P+CIAvAPjORf5uAI+UUu4jorsAfHRRdgXXXnvtUbqUgr29vaNf\n7pPS1quKTfuPFVpXXvOt665i3ubrfu2Wpy0R43FxKXwSDZ1E4t+9k59a6EUTdynmH5n07HXkI8It\nlU0R80idd5xIe60tv/ZIPurWI/XJ9hISbSK6AcB7AfwigJ9fFN8K4JZF+n4AjyEo2vWnWPf393Hq\n1KljosyFWxLyvb09VfCre2/3b89d4atVuBuP4DlYK4wSEfVe1+3ltZi5FTuXvglJoaXFOGnHjJiO\nuExP2D0hjwq299Ley8u1fGQr9clonVY2sk+EvDGsh6jT/mUA/xzAG5uy06WUSwBQSnmZiK7T3sxF\nWxJrSZi1fCvQddv+Jncta+sl4WgdtxQnb/eVPqCW45Ri49JLelJTK+sR9t7YuCTq/BuEFTu3Qiq9\nYZW233k+6sCnvqzjRc4ltY2X8escceMjTnwOsU3BXh+uaBPR3wNwqZTyNBG909hV/St+5jOfOUq/\n9a1vxVve8hYxDt2mNYftCTkXc/7Sji21hTtvzY1rjic62RkNsXihFp4eCaNYgq7Fy70Qi7SKRcv3\nOEVPuHk+IrTtOeYSfS8tbXvKtH2svosy6uaTfp5//nk8//zz7n4Rp/0OAD9GRO8F8HoA30FEnwTw\nMhGdLqVcIqLrAXxVO8C5c+eO0lXw2jhzK5ataPa8JEHXxFtz63wrue/aRiIKfXABe8LTmhysIRYv\n5i2J9ahT185lrZ7RRNwSdC3UYvVZO4astCaQvczl0q28l5auzeoDKa+VWeW9x1kGV9vN4cYbb8SN\nN954lP/c5z4n7ueKdinlHgD3AAAR3QLgn5VSfoqI7gNwB4B7AdwO4GHtGO3qEQDm0j3PgUdF3Qqn\naC5cEnfrfNZSREs4JEHXluJJSw+1ZYdaPJuLuSTumnhb/3jCio1b8XIptMIFXOoTDS9mXm+wfF9t\nn0pdrTQq3G2bpDbyc/aKtPfeaForm9OxT+FqE2+PKeu0fwnAg0R0J4CXAJzTduRL/lqhq+l2K7lw\nSdQlUeV5LRQiOW8pb7lx7YbSLjvUJjr5B1wTEy5q0lYLr3Dx1R7wsdIjoRUtNq658FonOXDNjWvz\nCNoHvO7PBdNyu3O/rGNLdbVM2mrUG4R2Y5COERXrUfFctuhebaLeJdqllN8B8DuL9J8CeHfoJKeO\nn0YbuK2YeW5cEkrNjVsxbC9OfvnyZdGFX758OeT4rXbzyc42Xbd1H0mg6tYTRO/VCvrrXvc61Y17\nK1W8MI73jcGbuLW+jVgOWspborhpAt7mtXQVZ7614IZBy2t9Ke3TW5/0s5bfHgH0+GPdcgHn26iY\naw68J7TiuW+r3vtG4F2n1GdS/1nxYikW3W6jq1emOPI2Lm45csmVt3n+LUMKrwD2b5a3aUlUNCFv\n03OI8+j7tTZ5aWkbTWt57WbJbwBa/Rz0TGKfBNYm2haRwWs9wl6KHGLhZVHh9gTdm/TUYuzcjUvt\njHyYAd2R1/63Jgd73XhEwK20J+rWZGdkwtN6VSxBr32r5UdFWRrDU0Tfa0vbbk/M+b7adXtpjiek\nU+tbRr4JbCMrEW2pM7WvrRVvUFmDP+LILQculUVCLNbWip9bx9W+ZbRlkf5q32NNdmpOd9SJR1er\nSE5cc+XWtwbNkbdhAy7ekiu3xrLnXqMiPvUlnStS1ualtDeerPQcZZE6i5Mo1C0r/z1tjud4JCID\nT5vkbIXPcrlRUfee4LTCJ57IS6EUz41L7k3q77bu4ODgmMDVp0pbwetx45Zo9y45jK5ekRx5xI3z\nSU9vXGqxYl7W44ijdXOKuybUlpDzfbU6Ky3ltTKJyH4p2jMgCUcpV8bDrNgjT1vHr2XWKzLRaYml\n58IjYt4TG4+EV6xwUdsnPM37i4iOHHntd08ILWHuceiWYHtOPerCtVcpr80NSONRC7VoeKLXK+jL\neFnjwhPxaJ32ObXEdQ5xPqnivRLR5hNq1qDvEXQp3+INJm85Hne0WohlND4eDa14+3rhnfa6pInP\nyOBvV7pUQecu1QplSLHxkXi4tYIlOtnZirhUJoVWuAvnY1RLR4Qp4niXKejScSPnjnzGeB9E0lIf\neWW9+24zaxVtzbVowt3W8zQvi/xhvQ9NFSgtnuyFUqLC7rnwSJm3gsX7RuF9oNt+qX+z1o23YRUu\nfvUVfUpTe9qzR8B7XtoqFu0bRr3+ti/aNB+XEUGPjks+Pnl+GQIe2X+OrdQXXj4i1MsW7lXfGDZG\ntDUhb8W7FWNP0Hti5PU8PO+95l566K0ZH0l7k51c1KXraj+gvK9qugp4fb802Xn58mXXjfcsKezN\n9wi75cSlZYY8Vs7HniTi3oc9IupWem5h144jnZ/va7Uxer1W2iqL1E3Zd5TRc2y0aEvibYkwd+i9\nsXEJbzBWkdLSUoilzUeFNRojjwi4JvTatwbLjbcx87Zva33N18nO/f191Y33xrRHQyie0EdCLNGJ\nTi7qfLz2jEnNQUqOdS4R9/bl9V4+IvaRa/P6xCuzyqP1vcx1vLVMRPK6qHBL5a0wtOkWL5wi7Re5\nhoiL0EIrNe3Fo3vE3JqsHBH1NtQSDa+0fcX7ozpwyYnzUITmdiUXPhI6GV1Pbn1TsGLj7Tiv8XHP\nXHhjUhqXmvBF0yNCPXJTkPbRxo60HUlH8tE6C00resyixdqcdt16otwr5O2HIyLobXusMikfGUhc\nvPnLCqv0Tnj2CHdkNYu0n7YEkt+UpK1E7Yc2tFIdOQ87RGPUUyc05wqtWCtXomObj9tRUfduqFaZ\nJ8a95REht9ombaN1Vtoqs8qj+7bjHJDn3qQyzkpEe39//1heE15etywRlz4MvH1t/WiYpcfdWCEW\nLpi8LPL05hRXHg2xSEsj27QnBDzN+78KtyWK1iRmVIxHRF8qk0I+tay9IUkxcj6upf6wxiD/8Hvp\nqLh743hUsCM3Cik9so2mpbxXHq2fwsqdNhfPWla3noBr5VEhb/PSudt0e2eMOhzpA2UNBG8Qe+GV\nKU58iphb+2kv6xraJYW8n/gH1npARgtZ1LQnwlMmOyOP6kfceBX13jGujW2e9saj1/9aOpIfEfTI\nua129lxj5D2cKQI+Iu4rj2lLQii53zmcON+/J6+lNSGXRFz68Gh90+NiJEGXnO3IksPRuLj2Xm+i\nk7dXu8Z2DLVLDXn/aqGI6GTnqvNSm2p52+6II7fEnI9LS8S18alte8VcG9feeO85rpSWttFrlPpE\n66tIeS2ztEFjLRORXPi4WNftqIBrdbzcylvnkuparD+Eta82OCKD2oona3HoKSGVHrcdeWKUu3Ep\nts/deN225dokZ+vKLSc+EhaZ25VHfzxLEnBt6aEl4u2DQxKW4xwR6Whdr8B7ZVK6Z6tdv5evfd7W\njYh1Za2PsfN9IkJet1Ehl8rmdOERsW+30nW1/aLlI07HEnRvwnMkvOKJ9sjj/dZNxbsxWYLS9qn1\nAJDmfCOv6GP3U0Is7Utq+7InO3ldRPBG0qOCrn0mIueR9pG2PWnvG41U5rFy0eaCLd2JLLHWtpaw\nzy3qVl2kTdKNSbquFu3OHP0Q8NCDNNnpTXqOhFKik6RRsZfCKu3Wu5kBOObGpb+/5nJbUR9dXTJ1\n0tMKr0TcuCXm2hi0XKEmSnOKea+4R8XbO75VJm152su3OtDDyn/lT3I+WqO5uPWI+qg79sR8VLQj\nLjz6gdHq+c2vTXsixl1rdLKzret9ajPqvHsnOqX2a2687ZuKtORQi41raU/M515HLpXPsfTQGpN8\nbPJ9JEPmCbsn8r0Cbo393htEtM3Wdbb9YX1eLdb6cE3PXUa64LnFXEtHBvSUG0W07bzPPPFu+87b\nti9pQpCHUiLryEdXqoyGWqSbibT1rrfGyInoSLzbvvWeiOSrVHpDKqNPdEppa9mhFCZqr69nnLdY\nQu85z8g45WNWq5PGdsTE9B5bSltbXiYJucXan4jkWEKuCXf7Xp6eKuS9rjvy/sh5pbS0la5d64uI\nE7AGds9kZ1TMp6wd73HhmiOXHHp77cDxiU5rLLTC1+PER8Itcz00JD11ai09lNw5H5ft51Qbl215\nZFzyv8mokEfH+mhZpL3Slqct1vpPEKQ7Cy9rPyht3jsef19U0LWyKSLdK+Q9Aq5dn3fXjoq65Uy4\nQ9VCFNHJzuivGs4t5jyWb92ceB9J/VzFTZvoXMaE59yP6XNRl66jFXFpzLfjk49bS9i9canlpzhj\nT5x7yqPt0so81v6fa6JIjlq6SE/0+TEiW09IvTprUEfLrbZFyngfSB8U7cMTdS58cpOnvbDFyI9c\n9TzkM7L0UEp7H1xvLFSxk5xuT1w8IsLeQ0Pez+Fak7LViWsTnZqg889oz7j0jIU0Tq0ybSxPEfLI\nZ0Wri7B20Y4Ib8USa+0PW+u0TqnnksS91+n2CLvnTnqOHU3ztkt9J5VJacnt7OzsuM7cmiSMhFZ6\nhbxXrCOO3JvolD6EbXn9m9TfH28Frg2veBOePULulUWeDrWcuDbpKd2stDFubSNjVXKsllC3fxct\nHXXj2pjvfW+EtYu2RlR4I0Lu3RhawZYcfa+bHRX1iHB7ea2uZyv1oZVu+9ByOVzAeZkXb/YEeQ43\nHhV9a5mkdY3AazFyouMTne3fzJvonBpO6QmjRH7DpWfpofZgkDem+biUyrjhiwi4NG6nirpVbh0j\nwsaKNid6QdofLFrPz9kKueVYNZHsFfWIoEf2j5xLa6N0fTxt7S+5zLr1BrT0wIwUXtHSlgufsiql\n15lbsX0pNg4cn/Dkf1/pyUceqpDCF3M88DM1by07rCEW7TotA9I7RkedeGQMS+M5Oub5Ph4h0Sai\nNwL4VQB/GcABgDsBPAfgAQBnALwI4Fwp5dXQWWekV8x7Hbl0LknIPaH2tqsW8mh7LSHnrkZ6j9S/\nkQ+EFjv2BDEqtqt4stO6wUgx/3bL+6V1520/a5OD1lK/kXDKSHiFx8W1F78OKR8Z93y8WWltbEbH\naTteI+ke5+0RddofB/AfSin/kIhOAXgDgHsAPFJKuY+I7gLwUQB3B4+3NDRR5uXSHViql+oksWpd\nuSaCWjriJCIiHBX0keNL1yPdvKQ+lK6X912bjgxyTfSsiU5LpKfGxHtWsezt7alhoJrWPtStS/f+\nztYkYiukIzHwkX20yU9vslP6liGltfFqpaWyEUdujV0treU9XNEmou8E8DdKKXcsTrAH4FUiuhXA\nLYvd7gfwGDZAtCWinQHEHLlW3wo3398aSJFtRNQjgh4p5/t4541el9afmrhbTtx6aeEIa4LTi5nP\n8SSn58Z5eXstUrpFEhGi12Lm2gM0kgOXnPEcwj7H0kPtFTUukfFpifoyHLn0t7OIOO3vAfA/iejf\nAvhrAJ4A8E8AnC6lXFqc7GUiui50xg1hLkduHdcTsh5xjwq655i9QT2Sj7RRuj7r+nm67VtpoHtu\nXBNzyeUu46dqe/fxvjFo4SMuBFpYBdAnOyU37on1qt04T2s3pt4xL33meHkdbz1CLo1bazxbRET7\nFIAfBPAzpZQniOiXceio+dHVs124cOEoffbsWZw9ezZw2tUTvdN5Qi79MaX9I+5UEreogLdlI2Iu\n1Uvnsc6rtUtqP09b+0sflJrWQgs8lKI5WckB83xkonN06WFk0tOb6PQmOwEcuXBJ4KSlh5pgzina\nI3kp9MNvQNp1auPbGrfWmB0JrdTtE088gSeffBIe5AkVEZ0G8F9KKd+7yP8IDkX7+wC8s5RyiYiu\nB/BoKeUKNSaicv78ebch24J0951SF3GZEcHrFfSoMEfcSTQfaVfkmi2iX1U1pxoRQ09YR3+TfFTI\ntQlP6Vql8IolSFY4onfp4VxCHllHLoV8apulCU9LyL2xKI1Lz7BZjrzy9re/HaWUKw7uOu2FKH+R\niN5aSnkOwLsA/NHidQeAewHcDuBh71gnAe0mx91fLbPeZznGUdG29hsRcl43IubROuu62r6KuB6t\n7614ohZasZYeStvok51TQyjeDcNa7+5Ndta093eXYuNarDwq2FPCKlbampC13Pjck51eaMUjunrk\nwwB+nYheB+B5AB8EsAvgQSK6E8BLAM4Fj3UiiXY4YIdPtD8mF62eeHh0G017gh4t986llfH69kMg\n9Sd/n/ShsWLjmqD3LD30BLZ32eFIWEUK+WhhIu7G277h/atNcta0FVLZ2ZF/EXEVP2fL2zF1sjPy\nueNp3pcRQqJdSvlDAD8kVL07dJarlIi7rmXae7w7sjdQtHi4lO4R9Ei+V7ij54tsteu30m3fepNG\nWoglMtHpCfeIaI+GWbwJT+C1OLi0JXpthYr29+OiKOU9cV3WWvLIk52WG5cmPb0xLX3metiaJyJP\nCr2TnZYb1+olUR8V94jz9YRbKuvZXzqv1x7vmrX+lPqE96nmwGtdZLme54CXsZY84sYld259s9Am\n3viEpxV+0ESSp9c1yVnLvZU20rV5Y94bmxIp2huA5K6joRPvPZHQiiduvULuOeaRcl7mncsq49fF\n+9FLL2uyk7vyXhFethvXHLk1F1D7o7pxqT8joQkvvDElNq45by82Lt1stEnO3vCKRYr2BjJXfLzW\nS46Rl/e61F4hb8si+R6X0nPT8No8MtnJ+16Lj/MyKYbcirokkm26Z8JzDjfu/XCWdgOy+iQ62amF\nK9rtMp7i9PLa8sPRyc4IKdpbRG9oRXuP5srb/T3HebU5cm+yk+f5N5y2b6OTnVZooifEsYxH9C03\nzvPeDaolMtkppSXhjDryKUI+5z+V2NnZEccXJ0V7y5lDyCUR1yZEe13qMoRcKht15L1OXLq+yHVL\nfTsy2SmJuBa20Jz31JUqvY7da6d3A4v8LTWH64VWRn8FMerAJXG3lkZGSNE+gcw92SntM+rIo/v0\nCq8n5t7+Peno9Wj9KaWl/o24cW3C03LkbfnIb5DPERvn5dY3i97JTm2Vh+TCNVGfO6wSXRmjfZtr\nSdG+SljWZKc2wbkpk509It+Tj7SZpyPXz/tVm+QEYE50agKuCfrU31WZ83dWuBPnIRXpWwjw2mP5\n+/v74tiZOtHZK+K9bjxFOzGZa7JTqptjsnNZjpznIwIfOZeU1rZzTXbWGKg24WnFketWEknJDc8R\nZokKecQq9kEVAAAL4klEQVSN90x2VqwVHL1Pdo66ce99KdpJN3M5cqmOT871um+pbA4hHymPHF/a\nx7oubbKzPQbfv+1nfuPsmey0liF6rznDKtEfy5JuNtYNiveL1MfSBKEm6HM68hTtZFbmmOhs6yVH\nzvfpdalThbxNR4W8x5H3OnF+fVIfRfbVJjd5nq+tlkRdcuPaZOeyf85WK9PaKsXGRyY7vd8uSdFO\nNhrNRXMijtw6Jt9vxH3zsoiojzjyUSG32sW3WmhFKpOcPO/bqBuXJgclx7sKN94TVmnbpl2DNdnZ\nivrOzg729+1/KmFNckZ+ZyVFO1kpq57s1Mp4vVY2p5Bb9dYxInWRrdbH0vXzfpUmOa3JTs+NS463\nx2HPJeTR31lpwyvStdb+qIJtjSHvUXcvnJKinayVZU92SgLf41I9QZ8i5FKZtY91TKs90Wvl1ynt\n39a17jLiyr3lh1Y8etN+Z0W6AUmTnQCOlVt/a2vpIX9i0iNFO1k5y3Lkdd9WhCwX3qa1sohARkQ4\nIujW8XpuJtbWC2PxvDbZ2Rte6Z3s7HXSy3iyk7dNukHxbx5tn0j9qT0R2Qq7R4p2snaWNdmpfXDa\n9Lodea+YW8eU6rxrlfBuahUtrDIy2emFLyQXvqzfILfWk3MXzh25dv3SmNL+1h4p2snGMvdkp+bM\n+X5Rsbb2jYh6r/NelSOPTnZK6TknOy03rgl6ND1F7LVvB96kreTAeVnEwKRoJ1vFsic72/f2um/t\nvZZARhzzVOG2jq+1S7tW65qlfZc52dmmR1enLCNerk10SgIuuXGPFO1kq5l7srOtn3OyUxP2aHrE\nkfccz0pHr0fK876R+jky2SnFkfnSQ+6AJUe8rCc7PSHn4q29IqRoJyeSTZnsjAif58S1dI9gR8u8\ndGTLxVlz7Lzveyc7NTeuOXNPVOf8JcSoG+f5DI8kyYJVTHZqohR13db7tDJP1EfFXDumtp/Vduum\npvUB71ttgq/mpYnOtkwKW2hryb0liHMLeVuWop0kDnNOdmofOCnM0iPkczhyS9Cj5ZFzae3i5ZHJ\nzqgzr+nIZKc16emFNlbhxlO0k2SASGgl4sitVQFT3Tcv6xXyNt3jyq28VtezlfrISlckJ96mrdBK\nK+BtOIWne4U7RTtJ1oQV1+ZEHHlbb7nwNh0JO/D3bYMjj1wrT2v52ifS30AKqViTnVLai0l7Yp6i\nnSRrZlmTna0Yb8Jkp+fMI/tbxx9pdyvOEUfeO9mpLUO0HnOPhFS8/GyiTUT/FMA/BnAA4BkAHwTw\nBgAPADgD4EUA50opr0aOlyQnkagoS2Ii1XuTnT1ibe0/IuQ8PyLm0eN71xMRbWlf6WbZE1qJPtnJ\n09aE5yyiTUR/EcDPAfiBUsq3iOgBAD8B4CYAj5RS7iOiuwB8FMDd7hmT5CoiumoFiIVW+Fd+7jzX\n7cinlvN9rDbwMl6uhbCk92n97Am4N+HZO9kZIRoe2QXwBiI6APB6AF/GoUjfsqi/H8BjSNFOEpc5\nHfnIZOeouC/Lkfcewzt35Fqs6+fptn/nnOyUJj1ncdqllK8Q0b8A8CcAvg7gs6WUR4jodCnl0mKf\nl4noOvdsSZKIzO3Ipfo5Jzt7HHmvmGt1keNr+0SvWevL9liRyU5N1L3fJ4kQCY/8WQC34jB2/SqA\n3ySinwTAR4w66i5cuHCUPnv2LM6ePRtqXJJczUQdeS3X3mMduxWhTQ6tjNwIpPNZ7eLX1xJx5FMn\nO5988kk89dRTV5z7irZ4f2Qi+gcAfrSU8tOL/E8B+GEAfwvAO0spl4joegCPllKuUGMiKufPn3cb\nkiTJGBGXGK2X3CUvj4Yg2rJlOPIesfbqotdlpTmStnpOvOXmm29GKeWKE0Ri2n8C4IeJ6M8A+CaA\ndwF4HMDXANwB4F4AtwN4OHCsJElmZs7QihRS4S58GY7cE9GoI4+U832s8/I6Xu5NdlrlmiP3iMS0\nP09EDwF4CsDlxfbfAPgOAA8S0Z0AXgJwLnTGJEmWziZOdlputUfQvboRRx49X2SrXb+V7iG0eqSU\n8gsAfoEV/ymAdw+dNUmSlRN1chEh1+qjQhUR9x5Hrom3VNazv3cuq0y7Zq0vo+QTkUlyFSM5516h\n9o7b6749kde2cwj5FEcupb3raokKeIp2kiRHjMTHtffVeh6z5TeGHrG2yiJuOJKPiLl1jGh7tOvy\nSNFOksQl6shrnfQe7X2jk529jjyajgi6VM7LvHNpZR4p2kmSdGMJtvYQj/Q+SeB7Jjv5OTXn6rnd\nHtFdhiOXrkUjRTtJktnQXDTHc+PeMfl+nuvm7Rh15FGXPcWRe6RoJ0myVHonO3vCKvy9kuu2xFoq\nmyu00iPmvE0WKdpJkqyUZf7Oiva7IJogj9aNOHKe18TcI0U7SZKNYFmOvO4rufA2vSmO3CNFO0mS\njWTuh4G0yc6oMPesaBkV8ggp2kmSbBVzT3Zqzpzvt4rQSoQU7SRJtp5VT3ZqZVKdtV+KdpIkCbZ7\nstMjRTtJkquGdUx2jkyAWqRoJ0ly1bKKyU4uxJHJTosU7SRJEsack52aEHNnHiVFO0mSJMAyQiu9\ngg2kaCdJkgyxjMnOCCnaSZIkMzKHI7dI0U6SJFkiI5OdFinaSZIka2Akng2kaCdJkmwEURHfWXI7\nkiRJkhlJ0U6SJNkiUrSTJEm2iJWL9sWLF1d9yiGynfOS7ZyPbWgjkO1cFinaCtnOecl2zsc2tBHI\ndi6LDI8kSZJsESnaSZIkWwSNLvAOn4BouSdIkiQ5oZRSrnhMcuminSRJksxHhkeSJEm2iBTtJEmS\nLSJFO0mSZItYqWgT0XuI6Fkieo6I7lrluXsgoheJ6A+J6Cki+vy621Mhok8Q0SUi+q9N2ZuI6LNE\n9N+J6LeI6I0b2MaPEdGXiOgPFq/3rLONizbdQET/kYj+iIieIaIPL8o3rT95O39uUb5RfUpE1xLR\n7y8+M88Q0ccW5RvTn0YbN6ovPVY2EUlEOwCeA/AuAF8B8DiA20opz66kAR0Q0fMA3l5K+d/rbksL\nEf0IgK8BOF9K+auLsnsB/K9Syn2LG+GbSil3b1gbPwbg/5VS/uW62sUhousBXF9KeZqIvh3AkwBu\nBfBBbFZ/au38cWxen35bKeXrRLQL4HcBfBjA38dm9afUxr+LDetLi1U67ZsB/HEp5aVSymUAv4HD\nwbeJEDYwdFRK+c8A+I3kVgD3L9L3A3jfShvFUNoIHPbpxlBKebmU8vQi/TUAFwHcgM3rT6mdb15U\nb1qffn2RvBaHP/tcsHn9KbUR2LC+tFilML0ZwBeb/Jfw2uDbNAqA3yaix4nop9fdGIfrSimXgMMP\nOIDr1twejZ8loqeJ6FfXHXLgENF3A3gbgN8DcHpT+7Np5+8vijaqT4loh4ieAvAygN8upTyODetP\npY3AhvWlxca5yQ3hHaWUHwTwXgA/s/jKvy1s4sL7fw3ge0spb8Phh2VjvoYuQg4PAfjIwsny/tuI\n/hTauXF9Wko5KKX8dRx+Y7mZiP4SNqw/hTbehA3sS4tVivaXAXxXk79hUbZxlFL+x2L7CoB/j8PQ\nzqZyiYhOA0fxz6+uuT1XUEp5pbw2efIrAH5one2pENEpHArhJ0spDy+KN64/pXZuap8CQCnl/wJ4\nDMB7sIH9CRxv4yb3pcQqRftxAN9PRGeI6BoAtwH49ArPH4KIvm3hakBEbwDwdwD8t/W26hiE4/G3\nTwO4Y5G+HcDD/A1r4FgbFx/WyvuxOf35awC+UEr5eFO2if15RTs3rU+J6C/UsAIRvR7A38Zh/H1j\n+lNp47Ob1pceK32MfbGU5uM4vFl8opTySys7eRAi+h4cuuuCw4mKX9+UdhLRvwPwTgB/HsAlAB8D\n8CkAvwngRgAvAThXSvk/G9bGv4nDWOwBgBcBfKjGOdcFEb0DwH8C8AwO/9YFwD0APg/gQWxOf2rt\n/EfYoD4lor+Cw4nGncXrgVLKLxLRn8OG9KfRxvPYoL70yN8eSZIk2SJyIjJJkmSLSNFOkiTZIlK0\nkyRJtogU7SRJki0iRTtJkmSLSNFOkiTZIlK0kyRJtoj/D5VBLwjQVpEBAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "synth = np.apply_along_axis(lambda t: np.convolve(t, w, mode='same'),\n", " axis=0,\n", " arr=rc)\n", "\n", "plt.imshow(synth, cmap=\"Greys\", aspect=0.2)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "That's it! And it only needed 9 lines of Python! Not incldung boring old imports and plotting stuff.\n", "\n", "Here they are so you can count them:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true }, "outputs": [], "source": [ "length, depth = 40, 100\n", "model = 1 + np.tri(depth, length, -depth//3)\n", "model[:depth//3,:] = 0\n", "rocks = np.array([[2700, 2750], [2400, 2450], [2800, 3000]])\n", "earth = np.take(rocks, model.astype(int), axis=0)\n", "imp = np.apply_along_axis(np.product, -1, earth)\n", "rc = (imp[1:,:] - imp[:-1,:]) / (imp[1:,:] + imp[:-1,:])\n", "w = bruges.filters.ricker(duration=0.100, dt=0.001, f=40)\n", "synth = np.apply_along_axis(lambda t: np.convolve(t, w, mode='same'), axis=0, arr=rc)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "
\n", "

© Agile Geoscience 2016

\n", "
" ] } ], "metadata": { "anaconda-cloud": {}, "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.6.8" } }, "nbformat": 4, "nbformat_minor": 1 }