Python: fastest way to map continuous coordinates to discrete grid?

Question

What's the fastest way to discretize continuous (x, y) or (x, y, z) coordinate data onto their nearest grid coordinates?

Said differently: what's the fastest way to turn the continuous, Euclidean position (6.778, 9.201) into a corresponding grid coordinate (7, 9) -- for example, the pixel coordinate in an image? Let's imagine I have 100 continuous (x, y) pairs and a known 150 x 150 discrete grid (e.g. 150^2 pairs of (0, 0), (1, 1), ..., etc. -- what's the fastest way to map each (x, y) pair to their nearest (grid1, grid2) coordinate?

Googling this didn't seem to give much results, but I'm sure this has been asked elsewhere, so apologies in advance if this is redundant.

Answer 1

Your best bet is to describe the grid in terms of its spacing and offset:

grid_offset = np.array([0.1, 0.3])
grid_spacing = np.array([1., 1.5])

In this notation, we take grid points to lie at

grid_offset + n * grid_spacing

For a given (x, y) pair, you can now trivially compute the nearest grid index with

point = np.array([x, y])
index = np.round((point - grid_offset) / grid_spacing)

To convert back into grid coordinates:

grid_offset + index * grid_spacing

If you have a large number of pairs, make the array shape (N, 2) (or 3, or however many dimensions you want). That way, you can broadcast directly with the (2,) grid arrays:

points = np.array([[x1, y1], [x2, y2], ..., [xn, yn]])
gpoints = grid_offset + np.round((points - grid_offset) / grid_spacing) * grid_spacing

However, your greatest speed boost would probably come from transposing the arrays (or defining them in Fortran order). Broadcasted operations on a (2, N) and (2, 1) array are going to do much fewer cache loads because of the contiguity of the dimensions being worked on.