Finding geographic coordinates of pixels of known value in set of images

Question

I am trying to know the geographic coordinates of pixels whose pixel value is known from its corresponding GeoTIFF image/array of interest.

clip = [[[1.0 2.3 3.3
          2.4 2.6 2.7
          3.4 3.2 3.9]]]
       [[[3.0 4.3 7.3
          2.7 2.6 5.7
          3.4 4.2 8.9]]]

#interested pixels whose geographic coordinates are needed from their corresponding image. First value from first image and second value from second image...

values = (3.2, 2.6)

Here, the clip stores all the arrays of interest. The shape of each array of interest/subset image is (1,3,3).

I got stuck here. I am unable to find the geographical coordinates of pixels of interest within these arrays of interest, the pixels of interest are given in values (1st value for the first array and so on).

The output should be the coordinates of 3.2 from the first array/image and 2.6 from the second array.

Expected outcome:

# longitude and latitude of points 3.2 from the first image and 2.6 from the second image stored in the "clip".

Coordinates = (64.3, 18.6)
              (97.2, 11.2)

Answer 1

Search for the row and column index of the interested values using np.where and use the acquired indices for getting the respective geographic coordinates using rasterio.transform.xy.

code:

meta = clipdataset.meta
meta['width'], meta['height'] = 3, 3
meta['transform'] = rio.windows.transform(window, dataset.transform) 

for a, d in zip(clip, values):
    rc = np.where(a == d)
    rowcol = [num for elem in rc for num in elem]
    interest = rio.transform.xy(meta['transform'], rowcol[0], rowcol[1], offset='center')
    print(interest)

Answer 2

There are a couple of ways you could do this. Here are the first two I thought of:

1. rioxarray

Using isel_window you could do:

import rioxarray

rds = rioxarray.open_rasterio(image)
...
subset = rds.rio.isel_window(window)

Then, you can get the x & y coordinate values from:

subset.x.values
subset.y.values

2. rasterio

You can use the xy method:

from itertools import product
...
row_range, col_range = window.toranges():
for row, col in itertools.product(range(*row_range), range(*col_range))
    x, y = dataset.xy(row, col)
...