1

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)

Any help in this regard is much appreciated.

2
  • In which way is your solution not working? – Stefan May 8 '20 at 7:24
  • I think you need to minimize the absolute difference ( min |a-b| ) – Stefan May 8 '20 at 7:25
2
+50

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)




1

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
  1. 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)
...
3
  • 1
    Can you please elaborate on the 2nd method to match best with my problem. I am unable to figure out how to give the interested values values in your partial code. – Sai Kiran May 11 '20 at 2:51
  • That complicates things a bit. What if the value appears more than once ... do you want the logic to pick only one or all of them? Do you want an exact match on the value or an approximate match within a tolerance? – snowman2 May 13 '20 at 1:15
  • I want the exact match on the value. I assure you that the values won't be repeated in each separate array. – Sai Kiran May 13 '20 at 1:38

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.