4

I'm getting a ZeroDivisionError from the following code:

#stacking the array into a complex array allows np.unique to choose
#truely unique points.  We also keep a handle on the unique indices
#to allow us to index `self` in the same order.
unique_points,index = np.unique(xdata[mask]+1j*ydata[mask],
                         return_index=True)
#Now we break it into the data structure we need.
points = np.column_stack((unique_points.real,unique_points.imag))

xx1,xx2 = self.meta['rcm_xx1'],self.meta['rcm_xx2']
yy1 = self.meta['rcm_yy2']
gx = np.arange(xx1,xx2+dx,dx)
gy = np.arange(-yy1,yy1+dy,dy)
GX,GY = np.meshgrid(gx,gy)

xi = np.column_stack((GX.ravel(),GY.ravel()))
gdata = griddata(points,self[mask][index],xi,method='linear',
                           fill_value=np.nan)

Here, xdata,ydata and self are all 2D numpy.ndarrays (or subclasses thereof) with the same shape and dtype=np.float32. mask is a 2d ndarray with the same shape and dtype=bool. Here's a link for those wanting to peruse the scipy.interpolate.griddata documentation.

Originally, xdata and ydata are derived from a non-uniform cylindrical grid that has a 4 point stencil -- I thought that the error might be coming from the fact that the same point was defined multiple times, so I made the set of input points unique as suggested in this question. Unfortunately, that hasn't seemed to help. The full traceback is:

Traceback (most recent call last):
  File "/xxxxxxx/rcm.py", line 428, in <module>
    x[...,1].to_pz0()
  File "/xxxxxxx/rcm.py", line 285, in to_pz0
    fill_value=fill_value)
  File "/usr/local/lib/python2.7/site-packages/scipy/interpolate/ndgriddata.py", line 183, in griddata
    ip = LinearNDInterpolator(points, values, fill_value=fill_value)
  File "interpnd.pyx", line 192, in scipy.interpolate.interpnd.LinearNDInterpolator.__init__ (scipy/interpolate/interpnd.c:2935)
  File "qhull.pyx", line 996, in scipy.spatial.qhull.Delaunay.__init__ (scipy/spatial/qhull.c:6607)
  File "qhull.pyx", line 183, in scipy.spatial.qhull._construct_delaunay (scipy/spatial/qhull.c:1919)
ZeroDivisionError: float division

For what it's worth, the code "works" (No exception) if I use the "nearest" method.

Community
  • 1
  • 1
mgilson
  • 300,191
  • 65
  • 633
  • 696
  • What are the relative ranges (e.g. max and min) of the data you're interpolating? I'm completely guessing, but I'm betting your problems stem from "thin triangles" in the delaunay triangulation used by `griddata` (Similar to: http://stackoverflow.com/questions/3864899/3867302 In extreme cases, you'll get division by zero instead of just wonky-looking interpolation.) Another work-around would be to use a radial basis function instead of delaunay triangulation: http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.Rbf.html#scipy.interpolate.Rbf – Joe Kington Jun 11 '13 at 17:56
  • @JoeKington -- if you're asking about `xx1`,`xx2` and `yy1`, the values are `-10,20,20` respectively, so I don't have a very extreme total aspect ratio. The actual data that is producing this error has a peak to peak difference of 17.68 and 21.21 in the `xdata` and `ydata` respectively. – mgilson Jun 11 '13 at 18:05
  • @JoeKington -- the `Rbf` might be what I end up using. My grid should be dense enough in the region where I'm interested to make it work. I would *rather* use delaunay tessellation as that will set values to numpy.nan outside my grid's domain (which is desirable) whereas I don't think the Rbf will do that for me... – mgilson Jun 11 '13 at 18:08
  • Well, in that case, I definitely guessed wrong. Yes, `scipy.interpolate.Rbf` interpolates everywhere instead of just inside the convex hull. You'd need mask interpolated points outside your domain after interpolation. The main advantage of an RBF is that it's "smooth" and in this case, that it shouldn't be particularly sensitive to closely-spaced points. I'd stick to `function='linear'` to make it even less sensitive to close-point-spacing issues. It's worth a shot. I'm still not sure why you're getting your original error, though... – Joe Kington Jun 11 '13 at 18:19
  • @JoeKington -- If you're like me and puzzles like these bug you until you understand them, I can post a link where you can get a very simple script/data (in `.npz` format) that reproduces the problem. However, if you want to pass on that one, that's fine with me too. – mgilson Jun 11 '13 at 18:41
  • @JoeKington -- Upon looking at the data more closely, it appears (to me at least) that the distribution of points may lead to "thin triangles" as you described it. I'll need to think about this a little more carefully I guess ... (I also tried `Rbf` and promptly ate all of the memory in my system ... so that approach probably won't work either ...) – mgilson Jun 11 '13 at 19:11
  • Scipy's RBF implementation is certainly memory-hungry. There are ways around it, (only use points inside of a search radius) but that adds a lot of complexity. At any rate, if you have a chance, post a link to the data or a subset of it, and I'd be glad to take a look at it. I won't be able to get to it until tonight, but my curiosity is certainly piqued! – Joe Kington Jun 11 '13 at 19:28
  • @JoeKington -- Here's a [link](http://iris.sr.unh.edu/~mgilson/reproduce_error/) to the data and script which can reproduce the error. If you're curious, the points are the intersection of magnetic fieldlines (starting on a simple grid at earth) with the equitorial plane. As the field gets bent and twisted by the solar wind, their mapped location in the equitorial plane gets scattered. – mgilson Jun 11 '13 at 19:37
  • 1
    Weirdly enough, It works perfectly for me... Here's the result: http://i.imgur.com/A4aJq2M.png – Joe Kington Jun 11 '13 at 19:57
  • 1
    Well ... That's obnoxious ... What version of numpy/scipy are you using? I suppose there could also be a difference in the Qhull library on the backend ... (I'm at scipy version 0.11.0, numpy 1.6.2) – mgilson Jun 11 '13 at 20:32
  • 3
    Turns out that upgrading to scipy 0.12.0 did the trick ... (weird)... – mgilson Jun 11 '13 at 21:15
  • Huh... Even wierder, because I'm using 0.11.0 and numpy 1.6.2. I'm completely stumped. Glad it finally worked! – Joe Kington Jun 11 '13 at 21:18
  • Perhaps pulling a fresh `scipy` also got me a fresh `qhull` -- Not sure if `qhull` can be installed independent of scipy, but if it can, there's a distinct possibility that I did it... Anyway, that might have been where the problem krept in at one point or another... – mgilson Jun 12 '13 at 01:15

0 Answers0