Get coordinates of all pixels in a raster with rasterio

Question

I have a n by m raster image in rasterio and I would like to get the ground coordinates of all the pixels in this raster (ideally in an array of tuples). What would be an efficient way to do this?

This is what I have so far which is not efficient at all

file_name = 'sample.tif'
image = rasterio.open(file_name)
height, width = image.height, image.width
X_coord = np.zeros((image.height, image.width))
Y_coord = np.zeros((image.height, image.width))
for row in range(height):
    for col in range(width):
        X_coord[row,col] = image.xy(row, col)[0]
        Y_coord[row,col] = image.xy(row, col)[1]

This piece of code loop through rows and columns in the image and uses the .xy(row, col) method to store x and y coordinates separately in two arrays. obviously, time-consuming and inefficient.

Answer 1

The function rasterio.transform.xy also accepts a sequence of rows and columns. So to avoid loops, create a meshgrid of every row and column.

       file_name = 'sample.tif'
       with rasterio.open(file_name) as src:
            band1 = src.read(1)
            print('Band1 has shape', band1.shape)
            height = band1.shape[0]
            width = band1.shape[1]
            cols, rows = np.meshgrid(np.arange(width), np.arange(height))
            xs, ys = rasterio.transform.xy(src.transform, rows, cols)
            lons= np.array(xs)
            lats = np.array(ys)
            print('lons shape', lons.shape)

lats and lons should each have the same shape as band1.

Thanks to MikeT's post for showing how to answer this.