Integrating 2D data over an irregular grid in python

Question

So I have 2D function which is sampled irregularly over a domain, and I want to calculate the volume underneath the surface. The data is organised in terms of [x,y,z], taking a simple example:

def f(x,y):
    return np.cos(10*x*y) * np.exp(-x**2 - y**2)

datrange1 = np.linspace(-5,5,1000)
datrange2 = np.linspace(-0.5,0.5,1000)

ar = []
for x in datrange1:
    for y in datrange2:
        ar += [[x,y, f(x,y)]]


for x in xrange2:
    for y in yrange2:
        ar += [[x,y, f(x,y)]] 

val_arr1 = np.array(ar)

data = np.unique(val_arr1)


xlist, ylist, zlist = data.T

where np.unique sorts the data in the first column then the second. The data is arranged in this way as I need to sample more heavily around the origin as there is a sharp feature that must be resolved.

Now I wondered about constructing a 2D interpolating function using scipy.interpolate.interp2d, then integrating over this using dblquad. As it turns out, this is not only inelegant and slow, but also kicks out the error:

RuntimeWarning: No more knots can be added because the number of B-spline
coefficients already exceeds the number of data points m.

Is there a better way to integrate data arranged in this fashion or overcoming this error?

It seems like you're trying to fit an exact solution to inexact (sampled) data, and scipy is recursing down the rabbit hole trying to find an exact solution. Honestly my first approach would just be to use `interp2d(kind='linear')` and re-sample a linear grid with a fixed dx, dy, add it up, and see if that gets you the precision you need. — Aaron, Mar 06 '19 at 15:15
I think this might be problematic as there is a very sharp feature (almost singular) at the origin of the function I want to sample. I was hoping that, so long as the interpolating function was built with data with high enough resolution, then `dblquad` could sample this efficiently. Fixing the grid step would require painfully small dx and dy. — Jiles, Mar 06 '19 at 15:21
will you always know the location of the sharp feature in order to sample it more thoroughly? — Aaron, Mar 06 '19 at 15:24
My next thought is to generate a voronoi diagram using the x,y points to get an area under each sample for `sum(z/area)`, though I'm not quite sure how the edges of the graph are handled by `scipy.spatial.Voronoi`... working on it... — Aaron, Mar 06 '19 at 15:29
I wouldn't use interpolation since you have a very sharp feature at the origin. Why not apply `dblquad` directly? — Tarifazo, Mar 06 '19 at 16:54
The actual function I need to evaluate is extremely slow and convoluted to find each point. Applying `dblquad` directly would be expensive — Jiles, Mar 06 '19 at 17:30

Aaron · Accepted Answer · 2019-03-15T18:46:39.227

If you can sample the data with high enough resolution around the feature of interest, then more sparsely everywhere else, the problem definition then becomes how to define the area under each sample. This is easy with regular rectangular samples, and could likely be done stepwise in increments of resolution around the origin. The approach I went after is to generate the 2D Voronoi cells for each sample in order to determine their area. I pulled most of the code from this answer, as it had almost all the components needed already.

import numpy as np
from scipy.spatial import Voronoi

#taken from: # https://stackoverflow.com/questions/28665491/getting-a-bounded-polygon-coordinates-from-voronoi-cells
#computes voronoi regions bounded by a bounding box
def square_voronoi(xy, bbox): #bbox: (min_x, max_x, min_y, max_y)
    # Select points inside the bounding box
    points_center = xy[np.where((bbox[0] <= xy[:,0]) * (xy[:,0] <= bbox[1]) * (bbox[2] <= xy[:,1]) * (bbox[2] <= bbox[3]))]
    # Mirror points
    points_left = np.copy(points_center)
    points_left[:, 0] = bbox[0] - (points_left[:, 0] - bbox[0])
    points_right = np.copy(points_center)
    points_right[:, 0] = bbox[1] + (bbox[1] - points_right[:, 0])
    points_down = np.copy(points_center)
    points_down[:, 1] = bbox[2] - (points_down[:, 1] - bbox[2])
    points_up = np.copy(points_center)
    points_up[:, 1] = bbox[3] + (bbox[3] - points_up[:, 1])
    points = np.concatenate((points_center, points_left, points_right, points_down, points_up,), axis=0)
    # Compute Voronoi
    vor = Voronoi(points)
    # Filter regions (center points should* be guaranteed to have a valid region)
    # center points should come first and not change in size
    regions = [vor.regions[vor.point_region[i]] for i in range(len(points_center))]
    vor.filtered_points = points_center
    vor.filtered_regions = regions
    return vor

#also stolen from: https://stackoverflow.com/questions/28665491/getting-a-bounded-polygon-coordinates-from-voronoi-cells
def area_region(vertices):
    # Polygon's signed area
    A = 0
    for i in range(0, len(vertices) - 1):
        s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
        A = A + s
    return np.abs(0.5 * A)

def f(x,y):
    return np.cos(10*x*y) * np.exp(-x**2 - y**2)

#sampling could easily be shaped to sample origin more heavily
sample_x = np.random.rand(1000) * 10 - 5 #same range as example linspace
sample_y = np.random.rand(1000) - .5
sample_xy = np.array([sample_x, sample_y]).T

vor = square_voronoi(sample_xy, (-5,5,-.5,.5)) #using bbox from samples
points = vor.filtered_points
sample_areas = np.array([area_region(vor.vertices[verts+[verts[0]],:]) for verts in vor.filtered_regions])
sample_z = np.array([f(p[0], p[1]) for p in points])

volume = np.sum(sample_z * sample_areas)

I haven't exactly tested this, but the principle should work, and the math checks out.

+1 for an elegant solution. I am, however getting the error `ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()` when run your example code. Do you get a similar issue? — Jiles, Mar 07 '19 at 10:10
I'm not sure I completely understand what's going on here. Why is it necessary to generate the lists `points_left` etc? Can I not just calculate the Voronoi diagram from the the sample points direclty? — Jiles, Mar 07 '19 at 10:59
@Jiles About the error: I didn't actually really test this at all, I assume it's on the line: `points_center = xy[np.where((bbox[0] <= xy[:,0] ...`? you may need to edit the syntax a little bit. Basically just a check to eliminate points outside the desired bounding box. — Aaron, Mar 07 '19 at 18:02
The general idea of the code is that voronoi cells define the area around a point for which that is the closest point. By definition, the regions on the perimiter extend out to infinity, so in order to bound the region, we use a symmetry trick to bound the region in a square. — Aaron, Mar 07 '19 at 18:10
The idea to use symmetry to bound the region is very clever, and I had no part in coming up with it ;P — Aaron, Mar 07 '19 at 18:13

Integrating 2D data over an irregular grid in python

1 Answers1

Linked