This is an archival dump of old wiki content --- see scipy.org for current material

Look at Cookbook/Interpolation before reading this. SciPy probably has the interpolation functions you need.

Here are some homemade interpolating functions I use that may be usefull to others. They may still contain a lot of bugs, but...

```   1 from numpy import *
2
3
4 def interpn_nearest( z, targetcoords, bincoords=None ):
5         '''Interpolate some data, z, located at bincoords in the target locations targetcoords
6         using nearest neighbor interpolation.
7                 z - is the data array
8                 targetcoords - are the coordinates where we want to evaluate the data.
9                         It must have a format so that z[targetcoords] would make sense
10                         if targetcoords where integers.
11                         mgrid and ogrid can be used to get such a format.
12                 bincoords - is a list of arrays. Each array should containt.
13         '''
14         coords = interpn_check_data(z,targetcoords,bincoords)
15         indices = [ c.round().astype('i') for c in coords ]
16         return z[indices]
17
18
19 def interpn_linear( z, targetcoords, bincoords=None ):
20         '''Interpolate some data, z, located at bincoords in the target locations targetcoords
21         using linear interpolation.
22                 z - is the data array
23                 targetcoords - are the coordinates where we want to evaluate the data.
24                         It must have a format so that z[targetcoords] would make sense
25                         if targetcoords where integers.
26                         mgrid and ogrid can be used to get such a format.
27                 bincoords - is a list of arrays. Each array should containt.
28         '''
29         coords = interpn_check_data(z,targetcoords,bincoords)
30
31         indices = [ x.astype('i') for x in coords ]
32         weights = [ x-i for x,i in zip(coords,indices) ]
33
34         res = z[indices]*0
35
36         for selector in ndindex( *z.ndim*(2,) ):
37                 weight = 1
38                 for w,s in zip(weights,selector):
39                         if s: weight = weight*w
40                         else: weight = weight*(1-w)
41                 value = z[ [ i+s for i,s in zip(indices,selector) ] ]
42                 res += weight*value
43
44         return res
45
46
47
48 def rebin( a, newshape ):
49         '''Rebin an array to a new shape, without interpolation.
50         This can be usefull if the array type is not a vector space.
51         '''
52         assert a.ndim == len(newshape)
53
54         slices = [ slice(0,old, float(old)/new) for old,new in zip(a.shape,newshape) ]
55         coordinates = ogrid[slices]
56         indices = [ c.astype('i') for c in coordinates ]
57         return a[indices]
58
59
60
61 #################################
62 ## more internal functions
63
64 def array_coordinates( targetx, binx ):
65         '''Computes the coordinates in an array reference.
66                 targetx are the coordinates of the points that we want to get in the array reference.
67                 binx    are the coordinates of the array bins
68         Example:
69                 >>> array_coordinates( [2, 10], [1, 3, 13] )
70                 [0.5, 1.7]
71         '''
72         tol = 1e-12
73         binx = asarray(binx)
74         nx = clip( targetx, binx.min()*(1+tol), binx.max()*(1-tol) )
75
76         i = searchsorted( binx, nx )      # index of the biggest smaller bin
77
78         bx = empty((len(binx)+2,), 'd')   # avoid probles at the border
79         bx[0] = binx[0] - 1
80         bx[1:-1] = binx
81         bx[-1] = binx[-1] + 1
82
83         return i-1 + ( nx - bx[i] )/( bx[i+1] - bx[i] )
84
85
86 def interpn_check_data( z, targetcoords, bincoords ):
87         '''Check the validity of the input data for intepn_x functions
88         and returns the target coordinates in the array reference
89         '''
90         dim = z.ndim
91         if bincoords:
92                 if len(bincoords) != dim:
93                         raise ValueError, 'bincoords shape mismatch.'
94                 for i in range(dim):
95                         if prod(bincoords[i].shape) != z.shape[i]:
96                                 raise ValueError, 'bincoords shape mismatch.'
97
98                 coords = [ array_coordinates(targetcoords[i],bincoords[i].ravel()) for i in range(dim) ]
99         else:
100                 coords = targetcoords
101         return coords
102
103
104 def interp2_linear( z, tx, ty, binx=None, biny=None  ):
105         '''Toy function just like interpn_linear in 2 dimensions.
106         This function exists just to help the understanding and mantaining of interpn_linear.
107         '''
108         if not binx is None: tx = array_coordinates( tx,binx )
109         if not biny is None: ty = array_coordinates( ty,biny )
110
111         ix = tx.astype('i')
112         iy = ty.astype('i')
113         wx = tx-ix
114         wy = ty-iy
115
116         return    (1-wy)*(1-wx) * z[iy  ,ix  ] \
117                 + (1-wy)*   wx  * z[iy  ,ix+1] \
118                 +    wy *(1-wx) * z[iy+1,ix  ] \
119                 +    wy *   wx  * z[iy+1,ix+1]
120
121
122
123
124
125
126 ###########################################
127 ## some visual testing
128 ##
129 import pylab
130
131 def test():
132         x,y = ogrid[ -1:1:10j, -1:1:10j ]
133         z = sin( x**2 + y**2 )
134
135         binx = (x,y)
136         tx = ogrid[ -2:2:100j, -2:2:100j ]
137
138         pylab.subplot(221)
139         pylab.title('original')
140         pylab.imshow(z)
141         pylab.subplot(223)
142         pylab.title('interpn_nearest')
143         pylab.imshow( interpn_nearest( z, tx, binx ) )
144         pylab.subplot(222)
145         pylab.title('interpn_linear')
146         pylab.imshow( interpn_linear( z, tx, binx ) )
147         pylab.subplot(224)
148         pylab.title('interp2_linear')
149         pylab.imshow( interp2_linear( z, tx[0],tx[1], x.ravel(),y.ravel() ) )
150         pylab.show()
```

SciPy: PauGargallo/Interpolation (last edited 2015-10-24 17:48:23 by anonymous)