Finding pixel location in raster using coordinates?

Question

I'm not a GIS person. I haven't found or at least identified a solution to my problem, so decided to ask.

I have a big rasters and a list of thousands of coordinates. I need to find the pixel location corresponding to a given lat/lon pair in the raster, and then clip an NxN rectangle with the pixel located in the middle of it and save the clip on disk. I haven't for the life of me figured out yet how I could find the pixel location using coordinates, and I've been on it a few days now. I found that with snappy I should be able to do something like that, but unfortunately it doesn't run on my system (anaconda3 with Python 3.6).

How can I go about solving this problem?

• Loïc Dutrieux: Going from coordinates in a given CRS to array indices and vice-versa is what the affine library does. Any rasterio dataset has an Affine which can be accessed via the transform attribute
• Michael Stimson: You have a tag for GDAL, is that the API you're interested in? In GDAL you open your raster and retrieve the GeoTransform which is 6 double values see the tutorial gdal.org/gdal_tutorial.html, column is ( X - GeoTransform[0] ) / GeoTransform[1] and row is ( GeoTransform[3] - Y ) / GeoTransform[5] because the reference is top left and Y cell sizes are negative.
• Michael Stimson: That formula only works if the raster is georeferenced docs.qgis.org/2.18/en/docs/training_manual/forestry/… and is in the same spatial reference as your input points.
• GISjoe: Sorry, should've been more explicit about the APIs. I'm interested in both GDAL and rasterio, although I prefer the latter. In any case, Luke's answer is exactly what I've been trying to do, and it's with rasterio. Thanks anyways for your help!

Answer 1

A rasterio way of doing this is pretty simple. Note this requires your raster be in the same projection as your coordinates. You can of course project your coordinates on the fly, but that's another question...

import rasterio as rio
infile = r"C:\Temp\test.tif"
outfile = r'C:\Temp\test_{}.tif'
coordinates = (
 (130.5, -25.5) , # lon, lat of ~centre of Australia
 (146.0, -42.0) , # lon, lat of ~centre of Tasmania
)
# Your NxN window
N = 3
# Open the raster
with rio.open(infile) as dataset:
 # Loop through your list of coords
 for i, (lon, lat) in enumerate(coordinates):
 # Get pixel coordinates from map coordinates
 py, px = dataset.index(lon, lat)
 print('Pixel Y, X coords: {}, {}'.format(py, px))
 # Build an NxN window
 window = rio.windows.Window(px - N//2, py - N//2, N, N)
 print(window)
 # Read the data in the window
 # clip is a nbands * N * N numpy array
 clip = dataset.read(window=window)
 # You can then write out a new file
 meta = dataset.meta
 meta['width'], meta['height'] = N, N
 meta['transform'] = rio.windows.transform(window, dataset.transform)
 with rio.open(outfile.format(i), 'w', **meta) as dst:
 dst.write(clip)

• GISjoe: Thanks, Luke, your answer is exactly what I've been trying to do with no success. I just tried it with the raster and coordinates I have, and it seems to be working like a charm.
• Derek Eden: do you if this returns nearest neighbour? I cant figure it out..wondering if this requires exact coordinates (within a certain precision) or returns NN for a given coordinate
• user2856: What do you mean "nearest neighbour"? If a coordinate is outside the extent of a raster?

Answer 2

Another way to do this is to use the rasterio.transform.rowcol() method described in the rasterio transform docs.

Example:

import numpy as np
import rasterio
xs = np.array([130.5, 146.0])
ys = np.array([-25.5, -42.0])
with rasterio.open("my.tif") as src:
 rows, cols = rasterio.transform.rowcol(src.transform, xs, ys)

Answer 3

Another option (as per March 2020) could be using xarray.open_rasterio in combination with the method .sel(). Like this you have yet another very short solution:

import xarray as xr
# We define a location
lon1, lat1 = (-0.25, 39.95) 
# In this example I am reading a file with (time, x, y) as dimensions
xarr = xr.open_rasterio(path_to_tiff_file)
# Slice one of the bands
img = xarr[0, :, :] 
#Use the .sel() method to retrieve the value of the nearest cell close to your POI
val = img.sel(x=lon1, y=lat1, method="nearest") 
print("This is val: ", val)

Which returns the following xarray.DataArray description:

>>> This is val: <xarray.DataArray ()>
array(0.249235)
Coordinates:
 band int64 1
 y float64 39.98
 x float64 -0.2087
Attributes:
 transform: (0.13190025669672106, 0.0, -60.553065717372434, 0.0, -0.1...
 crs: +init=epsg:4326
 res: (0.13190025669672106, 0.13190025669672106)
 is_tiled: 0
 nodatavals: (nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, n...
 scales: (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1...
 offsets: (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0...
 descriptions: ('0[-] SFC="Ground or water surface"', '0[-] SFC="Ground ...
 AREA_OR_POINT: Area

And if we open such TIFF file in QGIS and explore the location (i.e. 'Band 001') we can see its value:

enter image description here

So it seems to be working as expected. The only thing here is that I am not sure how fast this compact solution would be in the event of making this operation, say, thousands of times.

Answer 4

Another way to do this is to use the rasterio.transform.rowcol() method described in the rasterio transform docs.

I have been using this, but while dealing with tens of thousand of latitudes and longitudes profiling showed it to be very time-consuming as it was going in a linear loop.

Vectorizing it gave a pretty decent benefit; About 60 percent improvement in time for my use case

 lons = np.asarray(lons)
 lats = np.asarray(lats)
 #lon_lat = zip(lons,lats)
 # Invert the transform
 inv_transform = ~transform
 
 # Apply the transformation
 px, py = inv_transform * (lons, lats)
 
 return zip(np.floor(py).astype(int), np.floor(px).astype(int))

compared to linear

 x, y = rasterio.transform.rowcol(transform, lon, lat) 

Note that I am getting the transform object from rasterio

#bbox_obj = (3612, 3612, BoundingBox(left=-82.001666666667, bottom=38., right=-80.998368, top=40.00163))
transform = rasterio.transform.from_bounds(*bbox_obj[2], width=bbox_obj[0], height=bbox_obj[1])