-
Notifications
You must be signed in to change notification settings - Fork 3
/
interpolation-scipy.html
341 lines (280 loc) · 18.4 KB
/
interpolation-scipy.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
<!DOCTYPE html>
<html lang="en" itemscope itemtype="http://schema.org/Article">
<head>
<title>Interpolation methods in Scipy</title>
<meta charset="utf-8">
<meta property="og:title" content="Interpolation methods in Scipy">
<meta property="og:site_name" content="Modesto Mas | Blog">
<meta property="og:image" content="https://mmas.github.io/images/profile.jpg">
<meta property="og:image:width" content="200">
<meta property="og:image:height" content="200">
<meta property="og:url" content="https://mmas.github.io/interpolation-scipy">
<meta property="og:locale" content="en_GB">
<meta name="twitter:image" content="https://mmas.github.io/images/profile.jpg">
<meta name="twitter:url" content="https://mmas.github.io/interpolation-scipy">
<meta name="twitter:card" content="summary">
<meta name="twitter:domain" content="mmas.github.io">
<meta name="twitter:title" content="Interpolation methods in Scipy">
<meta name="description" content="Among other numerical analysis modules, scipy covers some interpolation algorithms as well as a different approaches to use them to calculate an inter...">
<meta name="twitter:description" content="Among other numerical analysis modules, scipy covers some interpolation algorithms as well as a different approaches to use them to calculate an inter...">
<meta property="og:description" content="Among other numerical analysis modules, scipy covers some interpolation algorithms as well as a different approaches to use them to calculate an inter...">
<meta name="keywords" content="interpolation,numerical-analysis,numpy,python,scipy">
<meta property="og:type" content="blog">
<meta name="viewport" content="width=device-width, initial-scale=1">
<meta property="og:type" content="article">
<meta property="article:author" content="https://github.com/mmas">
<meta property="article:section" content="interpolation">
<meta property="article:tag" content="interpolation,numerical-analysis,numpy,python,scipy">
<meta property="article:published_time" content="2015-10-28">
<meta property="article:modified_time" content="2015-10-28">
<link rel="stylesheet" type="text/css" href="/css/main.css">
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
CommonHTML: {
scale: 93,
showMathMenu: false
},
tex2jax: {
"inlineMath": [["$","$"], ["\\(","\\)"]]
}
});
</script>
<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-MML-AM_CHTML"></script>
</head>
<body class="entry-detail">
<header>
<div>
<img src="https://mmas.github.io/images/profile.jpg">
<a class="brand" href="/">Modesto Mas</a>
<span>Data/Python/DevOps Engineer</span>
<nav>
<ul>
<li><a href="/tags">Tags</a></li>
<li><a href="https://github.com/mmas/mmas.github.io/issues" target="_blank">Issues</a></li>
</ul>
</nav>
</div>
</header>
<section id="content" role="main">
<article>
<header>
<h1><a href="/interpolation-scipy">Interpolation methods in Scipy</a></h1>
<time datetime="2015-10-28">Oct 28, 2015</time>
<a class="tag" href="/tags?tag=interpolation">interpolation</a>
<a class="tag" href="/tags?tag=numerical-analysis">numerical-analysis</a>
<a class="tag" href="/tags?tag=numpy">numpy</a>
<a class="tag" href="/tags?tag=python">python</a>
<a class="tag" href="/tags?tag=scipy">scipy</a>
</header>
<aside id="article-nav"></aside>
<section class="body">
<p>Among other numerical analysis modules, <a target="_blank" href="https://www.scipy.org/scipylib/index.html">scipy</a> covers some interpolation algorithms as well as a different approaches to use them to calculate an interpolation, evaluate a polynomial with the representation of the interpolation, calculate derivatives, integrals or roots with functional and class-based interfaces.</p>
<p>Let's see some interpolation examples for one and two-dimensional data.</p>
<p>First of all, the required modules:</p>
<pre><code class="python">import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
</code></pre>
<h2>Univariate interpolation</h2>
<p>In the next examples, <code>x</code> and <code>y</code> represents the known points. We will need to obtain the interpolated values <code>yn</code> for <code>xn</code>. As a representation, <code>y0</code> will be the true values, generated from the original function to show the interpolator behavior.</p>
<p>If nothing about the plots is said, they will be generated as:</p>
<pre><code class="python">plt.plot(xn, y0, '--k', label='True values')
plt.plot(x, y, 'ok', label='Known points')
plt.plot(xn, yn, label='Interpolated values')
plt.legend()
plt.show()
</code></pre>
<h3>Linear interpolation</h3>
<p>The <a target="_blank" href="https://en.wikipedia.org/wiki/Linear_interpolation">linear interpolation</a> is easy to compute but not precise, due to the discontinuites at the points.</p>
<p>We'll do some examples with this values:</p>
<pre><code class="python">x = np.linspace(0, 4, 12)
y = np.cos(x**2/3+4)
xn = np.linspace(0, 4, 100)
y0 = np.cos(xn**2/3+4)
</code></pre>
<p><img alt="interpolation methods in scipy 001" src="/images/scipy_interpolate_001.png"></p>
<p>We can compute a linear interpolation with numpy:</p>
<pre><code class="python">yn = np.interp(xn, x, y)
</code></pre>
<p>It's time to introduce the scipy's one-dimension interpolate class. The (<code>interp1d</code>)[http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html#scipy.interpolate.interp1d] object will be created from the known points and we can obtain <code>yn</code> evaluating itself with the corresponding <code>xn</code>. <code>interp1d</code> offers different interpolation methods by the <code>kind</code> argument and the default is <code>linear</code>:</p>
<pre><code class="python">f = interpolate.interp1d(x, y, kind='linear')
yn = f(xn)
</code></pre>
<p><img alt="interpolation methods in scipy 002" src="/images/scipy_interpolate_002.png"></p>
<h3>Nearest-neighbor interpolation</h3>
<p>The univariate <a target="_blank" href="https://en.wikipedia.org/wiki/Nearest-neighbor_interpolation">nearest-neighbor interpolation</a> takes the same value of the closest known point:</p>
<pre><code class="python">f = interpolate.interp1d(x, y, kind='nearest')
yn = f(xn)
</code></pre>
<p><img alt="interpolation methods in scipy 003" src="/images/scipy_interpolate_003.png"></p>
<h3>Polynominal interpolation</h3>
<p><a target="_blank" href="https://en.wikipedia.org/wiki/Polynomial_interpolation">Polynominal interpolation</a> algorithms are computationally expensive and can present oscillator artifacts in the extremes due to the <a target="_blank" href="https://en.wikipedia.org/wiki/Runge_phenomenon">Runge's phenomenon</a>. Due to this, it is much better idea the use of <a target="_blank" href="https://en.wikipedia.org/wiki/Chebyshev_polynomials">Chebyshev</a> polynomials or interpolate using <a target="_blank" href="https://en.wikipedia.org/wiki/Spline_%28mathematics%29">splines</a> (more later).</p>
<p><a target="_blank" href="https://en.wikipedia.org/wiki/Lagrange_polynomial">Lagrange</a> or <a target="_blank" href="https://en.wikipedia.org/wiki/Newton_polynomial">Newton</a> are examples of polynomial interpolation. Just to mention and to introduce different interpolation problems approaches in scipy, let's see the Lagrange interpolation:</p>
<pre><code class="python">f = interpolate.lagrange(x, y)
yn = f(xn)
</code></pre>
<p>The <a target="_blank" href="https://en.wikipedia.org/wiki/Lagrange_polynomial#Barycentric_interpolation">barycentric interpolation</a> uses Lagrange polynomials. We can calculate the interpolated values directly with the interpolation functions:</p>
<pre><code class="python">yn = interpolate.barycentric_interpolate(x, y, xn)
</code></pre>
<p>Alternativelly, we can use the class-based interpolators to generate a polynomial from the known points and then, call this polynomial with our <code>xn</code> data:</p>
<pre><code class="python">f = interpolate.BarycentricInterpolator(x, y)
yn = f(xn)
</code></pre>
<p><img alt="interpolation methods in scipy 004" src="/images/scipy_interpolate_004.png"></p>
<p>The use of the class-based approach is recommended if we need to evaluate the <code>xn</code> data more than once, since we already have our polynomial calculated.</p>
<h3>Splines</h3>
<p>A spline is composed of polynomial functions connected by knots and, unlike the polynomial interpolation, does not present Runge's phenomenon, making the spline interpolation a stable and extended method of interpolation.</p>
<p>Let's change our data:</p>
<pre><code class="python">x = np.linspace(0, 2, 8)
y = 10*np.sinc(x*2+4)
xn = np.linspace(0, 2, 100)
y0 = 10*np.sinc(xn*2+4)
</code></pre>
<p><img alt="interpolation methods in scipy 007" src="/images/scipy_interpolate_007.png"></p>
<p>The easiest way to use splines in scipy is, again, with <code>interp1d</code>. Setting <code>kind</code> as <code>quadratic</code> or <code>cubic</code> we'll calculate the second and third order spline:</p>
<pre><code class="python">fq = interpolate.interp1d(x, y, kind='quadratic')
ynq = fq(xn)
fc = interpolate.interp1d(x, y, kind='cubic')
ync = fc(xn)
plt.plot(xn, y0, '--k', label='True values')
plt.plot(x, y, 'ok', label='Known points')
plt.plot(xn, ynq, label='Quadratic spline')
plt.plot(xn, ync, label='Cubic spline')
plt.legend()
plt.show()
</code></pre>
<p><img alt="interpolation methods in scipy 008" src="/images/scipy_interpolate_008.png"></p>
<p>Specifying an integer as a kind we'll set the order of the polynomials, taking into account that the order has to be lower than the number of known points:</p>
<pre><code class="python">f4 = interpolate.interp1d(x, y, kind=4)
yn4 = f4(xn)
f5 = interpolate.interp1d(x, y, kind=5)
yn5 = f5(xn)
f6 = interpolate.interp1d(x, y, kind=6)
yn6 = f6(xn)
f7 = interpolate.interp1d(x, y, kind=7)
yn7 = f7(xn)
plt.plot(xn, y0, '--k', label='True values')
plt.plot(x, y, 'ok', label='Known points')
plt.plot(xn, yn4, label='Spline order 4')
plt.plot(xn, yn5, label='Spline order 5')
plt.plot(xn, yn6, label='Spline order 6')
plt.plot(xn, yn7, label='Spline order 7')
plt.legend()
plt.show()
</code></pre>
<p><img alt="interpolation methods in scipy 009" src="/images/scipy_interpolate_009.png"></p>
<p><a target="_blank" href="https://en.wikipedia.org/wiki/Hermite_polynomials">Hermite polynomial</a> is related to Newton polynomial, it is a divided derivatives calculation. Matches de value of the n points and the and its first m derivatives, so the resulting polynomial will have a degree of, at most, n(m+1)-1.</p>
<p>The <a target="_blank" href="https://en.wikipedia.org/wiki/Cubic_Hermite_spline">cubic Hermite interpolation</a> consists in a spline of third-degree Hermite polymonials and the Hermite curves can be specified as <a target="_blank" href="https://en.wikipedia.org/wiki/B%C3%A9zier_curve">Bézier curves</a>, widely used in vectorial graphics design.</p>
<p>In scipy, the cubic Hermite interpolation has the two different approaches presented in the previous section, the functional interpolation:</p>
<pre><code class="python">yn = interpolate.pchip_interpolate(x, y, xn)
</code></pre>
<p>and the class-based interpolator:</p>
<pre><code class="python">f = interpolate.PchipInterpolator(x, y)
yn = f(xn)
</code></pre>
<p><img alt="interpolation methods in scipy 010" src="/images/scipy_interpolate_010.png"></p>
<p>As we can see, the interpolated values are quite different than the true values. Like in vectorial graphic design, the points needed for the Bézier curve has to be more significative, so in this case, we need to locate points in local maximimums or minimums. For the sinc function Asinc(Bx+C) we will find the peaks close to 1/2B, due to the C parameter, and then, close to every 1/2B, in this case {0.25, 0.75, 1.25, 1.75}. So, using the next points we will get a better result using the cubic Hermite interpolation:</p>
<pre><code class="python">x2 = np.array([0, .25, .75, 1.25, 1.75, 2])
y2 = 10*np.sinc(x2*2+4)
yn = interpolate.pchip_interpolate(x2, y2, xn)
</code></pre>
<p><img alt="interpolation methods in scipy 011" src="/images/scipy_interpolate_011.png"></p>
<p>Splines can also be calculated using the class-based interpolator, In this case with the <a target="_blank" href="http://www.netlib.org/fitpack/">FITPACK</a> interface:</p>
<pre><code class="python">spline = interpolate.InterpolatedUnivariateSpline(x, y, k=3)
yn = spline(xn)
</code></pre>
<p>And through the FITPACK functional interface:</p>
<pre><code class="python"># Get spline representation (knots, coefficients and degree).
tck = interpolate.splrep(x, y, k=3)
# Evaluate spline and its derivatives.
yn = interpolate.splev(xn, tck)
</code></pre>
<p><img alt="interpolation methods in scipy 012" src="/images/scipy_interpolate_012.png"></p>
<p>In both cases we will have many evaluation methods, for example, getting the roots of the function:</p>
<pre><code class="python">print spline.roots()
</code></pre>
<pre><code class="stdout">[ 4.60185838e-17 4.91329683e-01 1.00253142e+00 1.51200452e+00]
</code></pre>
<p>With FITPACK functional interface:</p>
<pre><code class="python">print interpolate.sproot(tck)
</code></pre>
<pre><code class="stdout">[ 4.60185838e-17 4.91329683e-01 1.00253142e+00 1.51200452e+00]
</code></pre>
<p>By the definition of the sinc function Asinc(Bx+C), roots will be found at every 1/B and at 0, due to the C parameter, which for the interval (0,2] will be {0, 0.5, 1, 1.5}.</p>
<h2>Multivariate interpolation</h2>
<p>Multivariate interpolation refers to a spatial interpolation, to functions with more than one variable. It is mainly used in image processing (<a target="_blank" href="https://en.wikipedia.org/wiki/Bilinear_interpolation">bilinear interpolation</a>) and geology elevation models (<a target="_blank" href="https://en.wikipedia.org/wiki/Kriging">Kriging interpolation</a>, not covered here).</p>
<p>First, let's go to define some data: <code>xn</code> and <code>yn</code> are the coordinates where we are going to interpolate our data, this coordinates are defined as well as a meshgrid (<code>numpy.mgrid</code>); <code>z0</code> corresponds to the true values for the coordinates; and <code>f0</code>, the function to define it (both, unknown in the practice).</p>
<pre><code class="python">def f0(x, y):
return np.sin(2*np.pi*x)**2 + np.sin(2*np.pi*y)**2
grid_xn, grid_yn = np.mgrid[0:1:200j, 0:1:200j]
xn = grid_xn[:,0]
yn = grid_yn[0,:]
z0 = f0(grid_xn, grid_yn)
plt.pcolor(grid_xn, grid_yn, z0)
plt.colorbar()
plt.show()
</code></pre>
<p><img alt="interpolation methods in scipy 013" src="/images/scipy_interpolate_013.png"></p>
<p>In the next examples, if nothing is said, the interpolation results will be displayed as a pseudocolor map as:</p>
<pre><code class="python">plt.pcolor(grid_xn, grid_yn, zn)
plt.colorbar()
plt.show()
</code></pre>
<p>For the next example we neeed to create some unstructured data, an array of points and the corresponding values:</p>
<pre><code class="python">points = np.random.rand(500, 2)
values = f0(points[:, 0], points[:, 1])
plt.scatter(points[:, 0], points[:, 1], c=values)
plt.colorbar()
plt.axis([0,1,0,1])
plt.show()
</code></pre>
<p><img alt="interpolation methods in scipy 014" src="/images/scipy_interpolate_014.png"></p>
<h3>Two-dimensional nearest-neighbor interpolation</h3>
<p>Multivariate nearest-neighbor is not computationally expensive, so it is used in real-time three dimensional rendering. The concept is equivalent to <a target="_blank" href="https://en.wikipedia.org/wiki/Voronoi_diagram">Voronoi diagram</a>.</p>
<p>Class-based interpolator:</p>
<pre><code class="python">f = interpolate.NearestNDInterpolator(points, values)
zn = f(grid_xn, grid_yn)
</code></pre>
<p><img alt="interpolation methods in scipy 015" src="/images/scipy_interpolate_015.png"></p>
<h3>Bilinear interpolation</h3>
<p>A bilinear interpolation is based in two linear interpolations in a 2D grid. It is used in image resampling.</p>
<p>Class-based interpolator:</p>
<pre><code class="python">f = interpolate.LinearNDInterpolator(points, values, fill_value=0)
zn = f(grid_xn, grid_yn)
</code></pre>
<p><img alt="interpolation methods in scipy 016" src="/images/scipy_interpolate_016.png"></p>
<h3>Bicubic spline interpolation</h3>
<p>It is based in two splines of order three in a two-dimensional grid. Used in image processing and GIS altitude data. It is smoother than bilinear and nearest-neighbor interpolation but more computationally expensive.</p>
<p>In this case we are going to use <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html#scipy.interpolate.griddata" target="_blank"><code>scpiy.interpolate.griddata</code></a> function, which creates directly a grid with the interpolated values in a similar way as <code>interp1d</code>. This function supports also linear and nearest-neighbor interpolations in two dimensions with the <code>method</code> argument.</p>
<pre><code class="python">zn = interpolate.griddata(
points, values, (grid_xn, grid_yn), fill_value=0, method='cubic')
</code></pre>
<p><img alt="interpolation methods in scipy 017" src="/images/scipy_interpolate_017.png"></p>
<h3>More about multivariate splines</h3>
<p>For the next examples we are going to need some structured data, values in a grid corresponding to ranked x an y coordinates:</p>
<pre><code class="python">grid_x, grid_y = np.mgrid[0:1:15j, 0:1:15j]
x = grid_x[:, 0]
y = grid_y[0, :]
z = f0(grid_x, grid_y)
plt.scatter(grid_x, grid_y, c=z)
plt.axis([0,1,0,1])
plt.colorbar()
plt.show()
</code></pre>
<p><img alt="interpolation methods in scipy 018" src="/images/scipy_interpolate_018.png"></p>
<p>Similar to the functional FITPACK interface for one-dimesional splines, here's the representation and evaluation of a bicubic spline using points and values as a grid data:</p>
<pre><code class="python">tck = interpolate.bisplrep(grid_x, grid_y, z, s=0, kx=3, ky=3)
zn = interpolate.bisplev(xn, yn, tck)
</code></pre>
<p>And the same with the class-based FITPACK interface. In this case we will need to calculate the spline with points not as a grid but as a ranked arrays:</p>
<pre><code class="python">spline = interpolate.RectBivariateSpline(x, y, z, s=0, kx=3, ky=3)
zn = spline(xn, yn)
</code></pre>
<p><img alt="interpolation methods in scipy 019" src="/images/scipy_interpolate_019.png"></p>
</section>
</article>
</section>
<footer></footer>
<script type="text/javascript" src="/js/article.js"></script>
</body>
</html>