Performing many means in numpy
Question: Good Morning, I am implimenting a Cressman filter for doing distance weighted
averages in Numpy.. I use a Ball Tree implimentation (thanks to Jake
VanderPlas) to return a list of locatations for each point in a request
array.. the query array (q) is shape [n,3] and at each point has the x,y,z at
point I want to do a weighted average of points stored in the tree.. the code
wrapped around the tree returns points within a certain distance so I get an
arrays of variable length arrays.. I use a where to find non-empty entries (ie
positions where there were at least some points within the radius of
influence) creating the isgood array...

I then loop over all query points to return the weighted average of the values
self.z (note that this can either be dims=1 or dims=2 to allow multiple co-
gridding)

so the thing that complilcates using map or other quicker methods is the
nonuniformity of the lengths of the arrays within self.distances and
self.locations... I am still fairly green to numpy/python but I can not think
of a way to do this array wise (ie not reverting to loops)

    
    
    self.locations, self.distances = self.tree.query_radius( q, r, return_distance=True)
    t2=time()
    if debug: print "Removing voids"
    isgood=np.where( np.array([len(x) for x in self.locations])!=0)[0]
    interpol = np.zeros( (len(self.locations),) + np.shape(self.z[0]) )
    interpol.fill(np.nan)
    for dist, ix, posn, roi in zip(self.distances[isgood], self.locations[isgood], isgood, r[isgood]):
        interpol[isgood[jinterpol]] = np.average(self.z[ix], weights=(roi**2-dist**2) / (roi**2 + dist**2), axis=0)
        jinterpol += 1
    

so... Any hints of how to speed up the loop?..

For a typical mapping as appied to mapping weather radar data from a
range,azimuth,elevation grid to a cartesian grid where I have 240x240x34
points and 4 variables takes 99s to query the tree (written by Jake in C and
cython.. this is the hard step as you need to search the data!) and 100
seconds to do the calculation... which in my opinon is slow?? where is my
overhead? is np.mean efficient or as it is called millions of times is there a
speedup to be gained here? would I gain by using float32 rather than the
default64... or even scaling to ints (which would be very hard to avoid wrap
around in the weighting... any hints gratefully recieved!


Answer: You can find a discussion about the relative merits of the Cressman scheme vs
using a Gaussian weight function at:

<http://www.flame.org/~cdoswell/publications/radar_oa_00.pdf>

The key is to match the smoothing parameter to the data (I recommend using a
value close to the average spacing between data points). Once you know the
smoothing parameter, you can set an "influence radius" equal to the radius
where the weight function falls to 0.01 (or whatever).

How important is speed? If you wish, rather than calling an exponential
function to determine the weight, you can make up a discrete table of weights
for some fixed number of radius increments, which speeds up the calculation
considerably. Ideally, you should have data outside the grid boundaries that
can be used in the mapping of the values surrounding the gridpoints (even on
the boundary points of the grid). Note this is NOT a true interpolation scheme
- it won't return the observed values at the data points exactly. Like the
Cressman scheme, it's a low-pass filer.

