scipy: RectBivariateSpline usage inconsistent with other interpolation functions and documentation
Let’s say that I want to interpolate data on a rectangular grid. For the purpose of illustration, I’ll define a function ‘f’ to generate the data:
x = arange(-2,10)
y = arange(3,12)
xx, yy = meshgrid(x,y)
f = lambda x,y: exp(x)-y
z = f(xx,yy)
If I want to interpolate this with interp2d, I’d do this:
exInterp = interp2d(x,y,z,kind='cubic')
print 'interp2d tests:'
print 'interpolated:', exInterp(2,3)
print 'exact:', f(2,3)
print 'interpolated:', exInterp(5,8)
print 'exact:', f(5,8)
This produces the expected results; the interpolated and exact answers are in agreement:
interp2d tests:
interpolated: [ 4.3890561]
exact: 4.38905609893
interpolated: [ 140.4131591]
exact: 140.413159103
However, RectBivariateSpline cannot be used in the same way:
In [2]: RectBivariateSpline(x,y,z)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-2-39b2819e680b> in <module>()
----> 1 RectBivariateSpline(x,y,z)
/usr/lib/python2.7/dist-packages/scipy/interpolate/fitpack2.pyc in __init__(self, x, y, z, bbox, kx, ky, s)
759 raise TypeError('y must be strictly ascending')
760 if not x.size == z.shape[0]:
--> 761 raise TypeError('x dimension of z must have same number of '
762 'elements as x')
763 if not y.size == z.shape[1]:
TypeError: x dimension of z must have same number of elements as x
Instead x and y must be swapped, both in creating and using RectBivariateSpline:
exInterp2 = RectBivariateSpline(y,x,z)
print 'RectBivariateSpline tests:'
print 'interpolated y,x:', exInterp2(3,2)
print 'interpolated x,y:', exInterp2(2,3)
print 'exact:', f(2,3)
print 'interpolated y,x:', exInterp2(8,5)
print 'interpolated x,y:', exInterp2(5,8)
print 'exact:', f(5,8)
This produces the output
RectBivariateSpline tests:
interpolated y,x: [[ 4.3890561]]
interpolated x,y: [[ 17.08553692]]
exact: 4.38905609893
interpolated y,x: [[ 140.4131591]]
interpolated x,y: [[ 2975.95798704]]
exact: 140.413159103
So, the (y,x) argument ordering produces the correct answer, while (x,y) does not.
Perhaps all this is intentional and is not a bug per se. However, it is inconsistent with x and y as defined by meshgrid. It is inconsistent with x and y as used by interp2d. And it is inconsistent with the documentation, which says that init and call both have arguments in the (x,y) ordering.
Because of these inconsistencies, I submit that this is a bug.
(FWIW, I’m using SciPy 0.12.0 on Kubuntu.)
About this issue
- Original URL
- State: closed
- Created 11 years ago
- Comments: 17 (7 by maintainers)
Commits related to this issue
- Added tests for plate_concentrated Fixed a bug when the number of elements was different in the two directions (see scipy/scipy#3164) and added tests based on the data of Timoshenko "Theory of Plates ... — committed to astrojuanlu/pfc by astrojuanlu 9 years ago
- Note in docs to address issue #3164. — committed to kyleaoman/scipy by kyleaoman 6 years ago
@codethief I’m glad the report was useful. However, it doesn’t look like the documentation has been changed to make this issue clear. That would be the real solution…
+1 to improving documentation of the issue, just fell prey to this inconsistency as well.
While this function is used in many contexts, would be possible to print a warning message when the larger axis (y, for example) becomes the lesser?
Here’s a situation to explain: Suppose x_before > y_before when interpolater was initialized. However, when interpolator runs, y_after > x_after.
interpolator = interpolate.RectBivariateSpline(y_before, x_before, my_data) scaled = interpolator(x_after, y_after)The implementation could detect that the larger dimension has now become the smaller. Additionally, when the order of parameters are transposed, the x/y ratios of the before and after would be different.
I realize this recommendation is specific to situations where the scaling code in question is maintaining aspect ratios, and that’s not always the case, but if that’s the most common one, the log warning message could save a lot of people time and trouble.
I stumbled across this issue for the second time, last time was 2 years ago. Remember to make a yourself note.
@lnmaurer: Yes, indexing=‘ij’ works perfectly for me. Another alternative is to transpose the coordinate grids returned by np.meshgrid, i.e.
xx, yy = xx.T, yy.T. I agree that this should be documented better. (Also, thank you for filing this bug because it saved me a lot of time searching for my mistake!)That sounds reasonable, but the the inconsistency between interp2d and RectBivariateSpline should be clearly noted in the documentation.
For example, the interp2d documentation has an example using “xx, yy = np.meshgrid(x, y)”, and the interp2d page suggests using RectBivariateSpline as a speedup. However, I see nothing that notes RectBivariateSpline must be used differently. That is misleading.
I ran into trouble when I tried to move from interp2d to RectBivariateSpline because of the noted performance increase. I was working with square arrays, so I didn’t get the TypeError; instead I had to pull my hair out figuring out why things weren’t working. (It took some time to realize that the interpolation was at issue because I wasn’t viewing that directly, and I made other changes before realizing that something was wrong; I searched everywhere else first because it looked like RectBivariateSpline was a drop-in replacement for interp2d.)
I see that meshgrid has two ordering options. If I switch to indexing=‘ij’, will I avoid this problem?
In short, while I now agree this is not a bug in the code, I still think it is a bug in the documentation.