使用三角桶的Python Binning

Question

我正在尝试找到一个简单的python模块/包，它已经实现了2D三角桶，这样它就可以以类似于scipy binned_statistic_dd的方式使用。有人知道这样的工具吗？我已经搜索过，但没有发现任何东西:我发现的最接近的是matplotlib的hexbin。

如果我必须创建一个自制的解决方案，为三角网格生成顶点点是很容易的，但是如果可能的话，如何有效地(如果可能的话需要避免慢循环，因为数据集大约是100 K点)搜索一个三角形所在的三角形？

Answer 1

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)

​

​