How to efficiently retrieve xy-coordinates from image

I have an image img with 1000 rows and columns each.
Now I would like to consider each pixel as x- and y-coordinates and extract the respective value.

An illustrated example of what I want to achieve:

In theory, the code snippet below should work. But it is super slow (I stopped the execution after some time).

img = np.random.rand(1000,1000)
xy = np.array([(x, y) for x in range(img.shape[1]) for y in range(img.shape[0])])
xy = np.c_[xy, np.zeros(xy.shape[0])]
for i in range(img.shape[0]):
    for j in range(img.shape[1]):
        xy[np.logical_and(xy[:,1] == i, xy[:,0] == j),2] = img[i,j]

Is there a faster way (e.g. some numpy magic) to go from one table to the other?

Thanks in advance!

A possible solution:

pd.DataFrame(img).stack().reset_index().to_numpy()

Method pd.DataFrame creates a dataframe where rows represent y-coordinates and columns represent x-coordinates. stack then compresses the dataframe's columns into a single column, turning the x-values into part of a hierarchical index alongside the y-values. reset_index flattens this multi-level index into columns, resulting in a dataframe with columns [y, x, val]. Lastly, to_numpy converts the dataframe into a numpy array.

Alternatively, we can use only numpy, through np.meshgrid and np.hstack (to horizontally concatenate the vertical vectors):

h, w = img.shape

y, x = np.meshgrid(np.arange(w), np.arange(h))

np.hstack([
    x.reshape(-1, 1),
    y.reshape(-1, 1),
    img.reshape(-1, 1)
])