Voronoi map

voronoi.py

import numpy as np

def voronoi(nx, ny, dh, xp, yp, val, xrng, yrng):
    """
    Simple Voronoi 2D interpolation.

    :param n1: number of points in the first dimension of the output model
    :param n2: number of points in the second dimension of the output model
    :param dh: spatial sampling
    :param xp: x-coordinates of points to interpolate
    :param yp: y-coordinates of points to interpolate
    :param val: value of the points to interpolate
    :param xrng: tuple (xmin, xmax)
    :param yrng: tuple (ymin, ymax )
    """

    # Get the number of points
    npts = np.size(xp)

    # Declare output array
    model = np.zeros((nx, ny), dtype=np.float32)

    # Loop over dimensions
    for ix in range(0, nx):
        x = xrng[0]+float(ix)*dh
        for iy in range(0, ny):
            y = yrng[0]+float(iy)*dh
            # Initialize minimum distance
            dmin = 0.
            # Calculate the distances
            d = np.sqrt((x-xp[:])**2+(y-yp[:])**2)
            # Get the minimum distance
            imin = np.argmin(d)
            dmin = d[imin]
            # Populate model
            model[ix, iy] = val[imin]

    return model

if __name__ == "__main__":

    import matplotlib.pyplot as plt

    # Parameters
    xmin = -1
    xmax = 1
    ymin = -1
    ymax = 1
    vmin = -10
    vmax = 10

    
    nx = 201
    ny = 201
    dh = 0.01

    # Randomize
    npts = 2000
    for ipts in range(npts):
        x = np.random.random_sample(npts)*(xmax-xmin) + xmin
        y = np.random.random_sample(npts)*(ymax-ymin) + ymin
        v = np.random.random_sample(npts)*(vmax-vmin) + vmin
    
    print(x, y, v)
    # Voronoi
    model = voronoi(nx, ny, dh, x, y, v, (xmin, xmax), (ymin, ymax))

    plt.imshow(model, aspect='auto', cmap='jet', interpolation="gaussian")
    plt.show()