6

I have 2048x2048 mesh of irregular data zi = f(xi, yi) which are essentially three independent sets of 2048 real values. I need to smoothly interpolate (perhaps bicubic spline) that into a regular mesh of wi = f(ui, vi) where ui and vi are integer values from 0 to 2047.

I have tried griddata which seems to work well at images less than 1000x1000, but blows up as you get to 1500x1500 (memory qhull errors for the Delaunay Mesh evidently). I have looked at some of ndimage functions, namely geometric_transform, RectBivariateSpline and map_coordinates, but they all seem to take regularized data as input. I could be missing something and just implementing it wrong though too!

I am trying to use Python/SciPy to do what this Matlab script I have been doing using tformarray and makeresampler. Any suggestions as to what function I can use to process this large data set? Thanks!

Paul Floyd
  • 5,530
  • 5
  • 29
  • 43
Robert Kraml
  • 61
  • 1
  • 2
  • 4
    I would look at this question: http://stackoverflow.com/questions/1972172/interpolating-a-scalar-field-in-a-3d-space I've used Shepard interpolation before with success and it might work for you. – Yann Aug 15 '11 at 18:22

1 Answers1

2

I tried to reproduce your errors without success. Are you on a 32 bit system? I had problems with scipy/numpy and large arrays so switched to 64 bit, and have had no problems since then.

Here's the code I used to try to reproduce your error (it will generate nothing useful, but should at least experience the same errors):

y,x=indices([2048,2048],dtype='float64')
z = randn(2048,2048)
yr = y + randn(2048,2048)
xr = x + randn(2048,2048)
zn = griddata(xr.ravel(),yr.ravel(),z.ravel(),x,y)
zl = griddata(xr.ravel(),yr.ravel(),z.ravel(),x,y,interp='linear')

This works on my machine.

If you can't run a 64 bit version of python (which can be difficult depending on what OS you're using), could you break your 2048x2048 grid into 4 1024x1024 grids?

keflavich
  • 18,278
  • 20
  • 86
  • 118
  • I had to turn the input points and the output points each into tuples to get griddata() to run. On my 2011 Macbook with 12GB your code did 1024x1024 in 33 seconds, and 2048x2048 in 184 sec – Dave X Jan 05 '17 at 15:21
  • The code does not function. Aside from the imports , you need to wrap the griddata() inputs together into one tuple as `(xr.ravel(),yr.ravel())` and same with the outputs: `(x,y)` – Dave X Jan 05 '17 at 18:53