How to get array indices from geo coordinates?

[MartinKist]

Hello,

For my project I want to calculate the sky-view-factor of different Points in Germany. For this I'm using a digital elevation model of Germany in Geotiff format. I also have a list of points (given as tuples of EPSG:3857 coordinates) of which I need to calculate SVF of. I already got my algorithm working, now I only need the matching array indices to my given geo coordinates.

I saw that there is the GeoTiff.get_coords(i, j) method to get geo coords from array indices, but I couldn't find a way to do it the other way around.

I already tried to implement it myself using the appropriate formulas of the Mercator projection, but I just couldn't get the latitude right :(

So my question is:
Is there already a way to achieve this using this package? And if not could this feature be added?

Best wishes,
Martin

[KipCrossing]

Hi @MartinKist, there are 2 hidden/internal methods that might help.

geo_tiff._get_y_int(lat)
geo_tiff._get_x_int(lon)

These methods round to the nearest int, but if you want a float, changes will have to be made.

Note this will only work if your as_crs code is the same as crs_code.

[MartinKist]

Thanks! I'll try that on Monday

[MartinKist]

thanks again, it worked flawlessly. But why is this not part of GeoTiff's public interface? 🤔

[KipCrossing]

Those methods have a very specific purpose, and as as part of the UI, you may get a lot of unexpected behaviour as they are not designed with the user in mind and to catch edge cases. Perhaps there could be a method that could help, something like:

x, y = geo_tiff.xy_indexes(lon, lat)

This would require some thinking as geotiffs may not be parallel to the lat, lon lines and the user would have to be given some control over the input crs.

I'm going to leave this thread here to see if there is further appetite for a feature like this.

[apwebber]

Just giving my two cents - I think this functionality would be very useful... pulling out point data from a geotiff is a common use-case for all sorts of things. I am making cross-section lines for example. I will also be extracting intersection intervals from horizon layers and borehole positions. I have subclassed and implemented my own solution (using the zarray) but a well-thought out feature in the package would be great.