3

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.

1 Answer 1

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.

1
  • perfect answer!
    – Shawn
    Apr 6 at 3:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.