Hi John On Tue, Oct 10, 2006 at 07:48:04PM +0100, John Travers wrote:

As far as I can tell from the documentation (and the fitpack author's book) surfit shouldn't be used for interpolation. It is designed for smoothing scattered data. I say this as the reason I implemented the regrid function is because interp2d wouldn't work for my data. regrid does (and *is* designed for interpolation or smoothing of rectangular data). Can I suggest that interp2d uses regrid instead (or do people need to interpolate scattered data, in which case a different library alltogether would be required)?

I'm willing to work out the required patches to make these changes if people think I should.

As you say, we should probably use the function purpose written for a rectangular grid in interp2d. However, I wouldn't like to lose the ability to smooth scattered data. What do you mean by "or do people need to interpolate scattered data, in which case a different library altogether would be required"? Isn't that functionality currently exposed via the surface fitting algorithm? Regards Stéfan

On Wed, Oct 11, 2006 at 10:59:37PM +0100, John Travers wrote:

Of course. But interp2d should be for interpolation only, and I think it is OK just to have it only for a rectangular grid. regrid is designed specifically for this (though it can also smooth data). The interface to surfit (which is bisplprep and SmoothBivariateSpline) will still be there for smoothing of scattered data, though it shouldn't be used for interpolation.

[...]

...

What do you mean by "or do people need to interpolate scattered data, in which case a different library altogether would be required"? Isn't that functionality currently exposed via the surface fitting algorithm?

No, the surfit routine can only smooth scattered data, to interpolate through scattered data points would require a different library from dierckx fitpack. Something like the other netlib fitpack or natgrid would be suitable.

OK, I understand what you meant now. I've applied your patch (btw, please keep line lengths < 79 chars). I still need to modify interp2d, which I will do on Friday, unless you submit a patch before then. Cheers Stéfan

On Thu, Oct 12, 2006 at 01:42:44PM +0100, John Travers wrote:

Thanks! Sorry for the long lines, I'll be stricter on the format next time. I can send a patch for interp2d this evening (London time) if you want.

That's great, I'll keep an eye out for it. If anyone on the list currently uses interp2d on irregularly gridded data, now is the time to speak up. Cheers Stéfan

I've added a patch to the ticket #286 that makes interp2d use regrid. I think this must be the best compromise for now until we sort out a suitable routine for scattered data interpolation. Regards, John On 12/10/06, John Travers <jtravs@gmail.com> wrote:

I started to look at this and think some extra work is required.

The problem is that regrid requires the data to be evenly spaced on a rectangular grid, rather than just a rectangular grid (which could have, say, increasing spacing in the y direction). I think this may limit its use more than users of interp2d would expect.

The problem then is:

1. The currently used function (dierckx surfit) is not suitable for interpolation and needs to be replaced (it segfaults my simple test).

2. regrid works very well (fast, low mem etc.) for regular rectangular grids, but can't handle varying rectangular grids, let alone scattered data.

3. uses probably expect that at least varying rectangular grids are handled

Solutions:

a) just use regrid and say we can't work with no regular data.

b) use a regridding method to get regular data if required, then use regrid

I favour the second. As it happens, this is what MATLAB does which I guess is a reasonable example to follow. The are several regridding methods. MATLAB uses the freely available (I think scipy compatible) Qhull code.

So I'll have a look at packaging the qhull code into scipy, then we can move on from there.

In the mean time, do we use regrid over the restricted domain, or use the current surfit (which only works in very few cases).

Regards, John

On 12/10/06, Stefan van der Walt <stefan@sun.ac.za> wrote:

...

On Thu, Oct 12, 2006 at 01:42:44PM +0100, John Travers wrote:

...

Thanks! Sorry for the long lines, I'll be stricter on the format next time. I can send a patch for interp2d this evening (London time) if you want.

That's great, I'll keep an eye out for it. If anyone on the list currently uses interp2d on irregularly gridded data, now is the time to speak up.

Cheers Stéfan _______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-dev

-- My mobile (which I rarely answer): (+44) (0) 7739 105209 If you want to see what I do: www.femto.ph.ic.ac.uk

-- My mobile (which I rarely answer): (+44) (0) 7739 105209 If you want to see what I do: www.femto.ph.ic.ac.uk

On Thu, Oct 12, 2006 at 08:15:46PM +0100, John Travers wrote:

I've added a patch to the ticket #286 that makes interp2d use regrid. I think this must be the best compromise for now until we sort out a suitable routine for scattered data interpolation.

I have a couple of concerns: 1. I recently made some changes to the interp2d interface in order to make the coordinate specification more consistent. Your patch simply reverts to the older behaviour. I don't like this because it assumes that x values increase with rows and y values with columns, even though this isn't mentioned anywhere in the docstrings. I'm not suggesting that the alternative is better, but it does make fewer assumptions -- it simply flattens the x and y inputs and pass them to the fitpack routine (unless x*y != z.size, in which case it makes up a meshgrid). 2. The docstrings were updated to be fairly descriptive, whereas after patching we'd be back to short, non-descriptive documentation. I'm also not happy about the way the output shape is determined (but this is a shortcoming of the current version, not a problem with your patch). Again, do x values increase along rows or columns? One option which overcomes this problem is for the output to be the same as the argument to __call__, i.e. x([[0,1,2],[3,4,5],[6,7,8]],[[0,1,2,],[3,4,5],[6,7,8]]) returns a (3,3) array of values, whereas x([0,1,2],[0,1,2]) returns a (3,) array of values, unlike the current (3,3). Unfortunately, this doesn't allow the usage of ogrid. Ultimately, it doesn't really matter as long as we choose an axis system and document it clearly. 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. We can probably sort out these issues with minimal effort. Thanks a lot for working on interp2d -- it has been long overdue for an overhaul. Cheers Stéfan

Hi John On Fri, Oct 13, 2006 at 12:48:56AM +0100, John Travers wrote:

...

1. I recently made some changes to the interp2d interface in order to make the coordinate specification more consistent. Your patch simply reverts to the older behaviour. I don't like this because it assumes that x values increase with rows and y values with columns, even though this isn't mentioned anywhere in the docstrings. I'm not suggesting that the alternative is better, but it does make fewer assumptions -- it simply flattens the x and y inputs and pass them to the fitpack routine (unless x*y != z.size, in which case it makes up a meshgrid).

The reason I changed this is because 'regrid' can only take 1d x,y arrays which have the same number of points as the respective axis of z. Flattening doesn't work here. It probably isn't hard for interp2d to take either coordinate arrays, mgrid matrices, meshgrid matrices or ogrid arrays and for us to detect which axis they are varying in and then get things right for regrid. But as per my note above, when we get scattered data working too we can deal with the input issues better. In the mean time I think it is better to have just regrid working.

Right, the relevant comments in regrid.f for anyone else following: c given the set of values z(i,j) on the rectangular grid (x(i),y(j)), c i=1,...,mx;j=1,...,my, subroutine regrid determines a smooth bivar- c iate spline approximation s(x,y) of degrees kx and ky on the rect- c angle xb <= x <= xe, yb <= y <= ye. c z : real array of dimension at least (mx*my). c before entry, z(my*(i-1)+j) must be set to the data value at c the grid point (x(i),y(j)) for i=1,...,mx and j=1,...,my. c unchanged on exit.

...

2. The docstrings were updated to be fairly descriptive, whereas after patching we'd be back to short, non-descriptive documentation.

I tried to only remove things that were no longer relevant (like the input array format). But if we change as per above, we should re-instate your docs. (I don't mean to tred on toes...). I can make clearer docstrings for the temporary regrid only interp2d first.

No! Please go ahead and forget about the toes. I only ask these questions to make sure that we are on the same page. I want the documentation to clearly state exactly what we are doing, so that there is no confusion on how to use the new interface (I've never seen a ticket complaining about the detailed level of docstrings ;).

...

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.

OK, unless someone beats me to it, I will merge the patch tomorrow morning (ZA time). Cheers Stéfan

Hi John On Fri, Oct 13, 2006 at 03:42:13PM +0100, John Travers wrote:

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 think matplotlib may be throwing a spanner in the wheels here by using its own interpolation. If you use imshow(x,interpolation='nearest') you get a plot like http://mentat.za.net/results/interpolate.png (I changed to a more densely sampled grid). As for the meshgrid behaviour, take a look at numpy's mgrid, which doesn't do the argument swapping.

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...)

I'm sorry, I still havn't had time to merge the patch -- will hopefully be able to do that over the weekend at least.

Hope this demonstrates my point (and that yopu don't find an obvious error in my code...)

Thanks for the demo! Pretty pictures always make for a convincing argument. Cheers Stéfan

On 03/01/07, Tim Leslie <tim.leslie@gmail.com> wrote:

Hi All, I've been working through a bunch of tickets and noticed this one has stagnated. Can someone comment on the status of this.

http://projects.scipy.org/scipy/scipy/ticket/286

The second patch to this ticket has not been applied? Is there a reason not to? I noticed the other thread on this topic suggested that a complete solution might require moving delaunay out of the sandbox. Is this still the case?

Cheers,

Tim

Hi Tim, everybody, I am (slowly) trying to work on improving the interp2d function as well as spline fitting in general. The first patch attached to the ticket you mention is applied and provides the RectBivariateSpline class for fitting a spline to rectangular gridded data. The second patch changed interp2d to use this class rather than the firtpack surfit function (which absolutely should not be used for interpolation). This was not applied as it was decided that interp2d really must support non-uniform grid data. The only solution to this problem (surft still needs replacing) proposed to date is to use Rokert Kern's delaunay package from the sandbox. In order to achieve this I proposed that the spline *fitting* code be taken out of scipy.interpolate and put into scipy.spline. And that scipy.interpolate should import methods from both the spline and delaunay package to provide interp1d, interp2d, interpnd only. My first step towards this was to create a spline package in the sandbox. In doing so I got side tracked into cleaning it up (moving to purely f2py etc.) And then I got seriously side tracked at work (PhD thesis work), so my timing wasn't great. TODO: 1. finish the sandbox.spline package (I could skip cleaning it up for now, I'm about half way through on my machine (uncommitted) - do people think moving to clean f2py interface is really worth it?) 2. create a sandbox.interpolate package that contains the new interp functionality based on spline and delaunay only 3. after testing, move all three packages out of the sandbox What do people think about all of this? It will lead to quite large API changes (but on a high level, i.e. scipy.interpolate.RectBivariateSpline will become scipy.spline.RectBivariateSpline) Cheers, John