[SciPy-dev] new function regrid for fitpack2

[John Travers]

Hi Stefan,

On 13/10/06, Stefan van der Walt <stefan at sun.ac.za> wrote:
> > > I'd also like to see some tests to verify the behaviour of the new
> > > interp2d that is *different* from the old version.  It's easier to add
> > > it now, while everything is still fresh in our memories.
> >
> > Will try and do this tomorrow.
>

Attached is a very simple script (will make it a test later) which
uses the test data from netlib->dierckx to check the use fo regrid and
surfit. The plot attached shows the results. There could be an error
in my dealing with meshgrid (I hate they way it swaps axes), but I
think it is right. The surfit part issues a warning on interpolation.
If you increase s to say 20 (0 is needed for interpolation) then the
warning goes, but then you are smoothing (as can be seen from the
corresponding output plot - which was down sampled for file size).

I'll make this into a test at some point and also improve the interp2d
interface to be more flexible to input array layout - but this will
have to be nest week due to time constraints (I need to work...)

Hope this demonstrates my point (and that yopu don't find an obvious
error in my code...)
John
-------------- next part --------------
A non-text attachment was scrubbed...
Name: intplot2.png
Type: image/png
Size: 22948 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20061013/e3e1d600/attachment.png>
-------------- next part --------------
# very simple script to test interpolation with regrid and surfit
# based on example data from netlib->dierckx->regrid
 
from numpy import *
from scipy.interpolate.fitpack2 import SmoothBivariateSpline, \
                                       RectBivariateSpline
import matplotlib
matplotlib.use('Agg')
import pylab

# x,y coordinates
x = linspace(-1.5,1.5,11)
y = x                        
# data taken from daregr
z = array([
[-0.0325, 0.0784, 0.0432, 0.0092, 0.1523, 0.0802, 0.0925, -0.0098, \
                                                    0.0810, -0.0146, -0.0019],  
[0.1276, 0.0223, 0.0357, 0.1858, 0.2818, 0.1675, 0.2239, 0.1671, \
                                                    0.0843, 0.0151, 0.0427],  
[0.0860, 0.1267, 0.1839, 0.3010, 0.5002, 0.4683, 0.4562, 0.2688, \
                                                    0.1276, 0.1244, 0.0377],  
[0.0802, 0.1803, 0.3055, 0.4403, 0.6116, 0.7178, 0.6797, 0.5218, \
                                                    0.2624, 0.1341, -0.0233],  
[0.1321, 0.2023, 0.4446, 0.7123, 0.7944, 0.9871, 0.8430, 0.6440, \
                                                    0.4682, 0.1319, 0.1075],  
[0.2561, 0.1900, 0.4614, 0.7322, 0.9777, 1.0463, 0.9481, 0.6649, \
                                                    0.4491, 0.2442, 0.1341],  
[0.0981, 0.2009, 0.4616, 0.5514, 0.7692, 0.9831, 0.7972, 0.5937, \
                                                    0.4190, 0.1436, 0.0995], 
[0.0991, 0.1545, 0.3399, 0.4940, 0.6328, 0.7168, 0.6886, 0.3925, \
                                                    0.3015, 0.1758, 0.0928],  
[-0.0197, 0.1479, 0.1225, 0.3254, 0.3847, 0.4767, 0.4324, 0.2827, \
                                                    0.2287, 0.0999, 0.0785],  
[0.0032, 0.0917, 0.0246, 0.1780, 0.2394, 0.1765, 0.1642, 0.2081, \
                                                    0.1049, 0.0493, -0.0502],  
[0.0101, 0.0297, 0.0468, 0.0221, 0.1074, 0.0433, 0.0626, 0.1436, \
                                                    0.1092, -0.0232, 0.0132]])

# plot original data
pylab.subplot(1,3,1)
pylab.imshow(z)
pylab.title('orig')

# check regrid
mrs = RectBivariateSpline(x,y,z)
zr = mrs(x,y)
print sum(abs(zr-z))
pylab.subplot(1,3,2)
pylab.imshow(zr)
pylab.title('regrid')

# check surfit
# x increases in columns, y increases along rows
ym,xm = meshgrid(y,x) # deal with meshgrid badness
mrs = SmoothBivariateSpline(ravel(xm),ravel(ym),ravel(z),kx=3,ky=3,s=0)
zr = mrs(x,y)
print sum(abs(zr-z))
pylab.subplot(1,3,3)
pylab.imshow(zr)
pylab.title('surfit')

pylab.savefig('intplot.png')
pylab.close()