Python Binning using triangular bins

Question

I'm trying to find a simple python module/package that has implemented 2D triangular bins so that it can be use in a similar fashion to scipy binned_statistic_dd. Is anyone aware of such a tool? I've searched but not found anything: the closest I've found is matplotlib's hexbin.

If I have to create a home-made solution, generating the vertex points for the triangular grid is easy, but how would you efficiently (need to avoid slow loops if possible as datasets are about 100K points) search which triangle a point lies in?

To get any kind of speed you'll need compiled code. You probably need to search outside of a python/numpy/scipy context for ideas. The delauny tringulation/convex hull code that `scipy` provides uses c code that was available as a standalone source years ago. — hpaulj, Jul 18 '22 at 20:26

score 0 · Accepted Answer · answered Jul 21 '22 at 16:31

import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
def plot_triangular_bin_freq(x,y,Vx,Vy):
    X, Y = np.meshgrid(x, y)
    Ny, Nx = X.shape
    
    iy,ix = np.indices((Ny-1, Nx-1))
    # max vertice is supposed to be 
    # max(iy)*Nx + max(ix) + (Nx+1)
    # = (Ny-2)*Nx + (Nx-2) + (Nx+1)
    # = Ny * Nx - 1
    assert iy.max() == Ny-2
    assert ix.max() == Nx-2
    
    # build square grid and split it in a lower-left, upper-right triangles
    # and construct the triangulation
    vertices = (((iy * Nx) + ix)[:,:,None] + np.array([0,1,Nx,Nx,Nx+1,1])[None,None,:]).reshape(-1, 3)
    triangles = tri.Triangulation(X.flatten(), Y.flatten(), vertices)
    
    # Normalized point coordinates
    Vx = (np.asarray(Vx).flatten() - x[0]) * ((Nx-1) / (x[-1] - x[0]))
    Vy = (np.asarray(Vy).flatten() - y[0]) * ((Ny-1) / (y[-1] - y[0]))
    
    m = (0 <= Vx) & (Vx < Nx-1) & (0 <= Vy) & (Vy < Ny-1)
    
    # get indices on the x,y boxes
    Ix, Rx = divmod(Vx[m], 1)
    Iy, Ry = divmod(Vy[m], 1)
    
    # (Rx+Ry)=1 is the boundary between the two triangles
    # w indicates the index of the triangle where the point lies on
    w = ((Rx+Ry)>=1) +  2*(Ix + (Nx-1)*Iy)
    assert max(Ix) < Nx-1
    assert max(Iy) < Ny-1
    assert max(Ix + Iy*(Nx-1)) < (Nx-1)*(Ny-1)
    
    # z[i] is the number of points that lies inside z[i]
    z = np.bincount(w.astype(np.int64), minlength=2*(Nx-1)*(Ny-1))
    plt.tripcolor(triangles, z, shading='flat')

x = np.arange(15)/2.
y = np.arange(10)/2.
Vx = np.random.randn(1000) + 3
Vy = np.random.randn(1000) + 1

plot_triangular_bin_freq(x,y,Vx,Vy)

Thanks, this is a nice solution that is much better than manually constructing a series of triangles and looping. Since it based on a rectangular grid, it may not suit my most common use of equilateral triangles, but I think there's a way probably modify this approach to make it work. — jpmorr, Jul 28 '22 at 15:00
Yes, if you do a non singular affine transformation to the data and to the triangles the count is preserved. So if you have data and a triangle mesh you apply a transform so that the triangles align to this lattice. Compute the counts and then apply the inverse transform to the triangles and you know the counts for each triangle. — Bob, Jul 28 '22 at 19:14

Python Binning using triangular bins

1 Answers1