from functools import reduce
def cartesian_product(arrays):
    broadcastable = numpy.ix_(*arrays)
broadcasted = numpy.broadcast_arrays(*broadcastable)
rows, cols = reduce(numpy.multiply, broadcasted[0].shape), len(broadcasted)
out = numpy.empty(rows * cols, dtype=broadcasted[0].dtype)
start, end = 0, rows
for a in broadcasted:
    out[start:end] = a.reshape(-1)
start, end = end, end + rows
return out.reshape(cols, rows).T
Numpy: cartesian product of x and y array points into single array of 2D points