Find indices using for loop
If you are trying to make it fast, the last thing to do is drop into Python for-loops.  Keep it all in numpy as much as possible.  

The code below maps the particles to the mesh and calculates the mean v for all the particles in each mesh square.  With 400M particles having random x,y, and v on a 1600x900 size mesh, it runs in about 100 seconds on my laptop.

xi = np.digitize(x, np.linspace(xmin, xmax, nx, endpoint=False)) - 1
yi = np.digitize(y, np.linspace(ymin, ymax, ny, endpoint=False)) - 1

flat_index = nx*yi + xi

counts = np.bincount(flat_index)

vsums = np.bincount(flat_index, weights=v)

means = (vsums / counts).reshape(nx,ny)


np.digitize(x, bins) maps values the values in the array x, to ints based on the values in bins.  x &lt; bin[0] maps to 0; bin[0] &lt;= x &lt; bin[1] maps to 1; and so on.  The endpoints=False makes sure the indices go from 1 to 1440000 instead of 1 to 143999; and the -1 shifts it to 0 to 1439999.

np.linspace(start, end, nsteps) divides the range [start, end) in to nstep even intervals.  This is used to define the bins for np.digitize().

nx*yi + xi basically converts the x and y indices to a 1D index from 0 to 1439999, because bincount only works on 1D arrays

np.bincount(z, weights) if weights is omitted, this counts the number of occurrences of each integer in z (like a histogram).  If provided, weights is an array of the same dim as z, the values of the corresponding weights are summed instead of just counting the integers.  So setting weights to v, sums the velocity of each particle in a mesh square.

(sums / counts).reshape(nx,ny) calculates the means and reshapes the result into the 1600x900 mesh.