# Mapping spatial corrdinates (X,Y) back to pixel space using Rasterio

For the problem I am solving I have lat/lon values in DataFrame and need to check if these values coincide with some of the satellite images I have. This part I managed to solve following this post (https://stackoverflow.com/questions/64282097/how-to-check-given-a-coordinate-lat-long-if-is-in-a-raster-image-file-tif-fil).

Following this I need extract a portion from the image, which needs to be 256x256 pixels with the point at the center. In other words the satellite images I have are quite large and my task is to extract 256x256 image patch from this larger image where the point falls at the center of this 256x256 image.

``````import rasterio
import pyproj
import pandas as pd

# Lets just consider one image for now
dataset = rasterio.open("path to image")
points = pd.read_csv("path to csv containing gps coordinates")

pp = pyproj.Proj(init = "epsg:32650")

# loop through each row of the dataframe
for i,row in points.iterrows():
#Project lat/lon values to x,y
px,py = pp(row.GPS_X,row.GPS_Y)

#Check if point falls/coincide with the image
if (dataset.bounds.left < px < dataset.bounds.right) and (dataset.bounds.bottom < py < dataset.bounds.top):

# Need some sort of inverse transformation function here to convert px,py or GPS_X,GPS_Y values to pixel values
# Lets calls these corresponding values ppx, ppy
from_row, to_row, from_col, to_col = int(ppx - 128), int(ppx + 128), int(ppy -128) , int(ppy + 128)

# get rgb image values
rgb_img = dataset.read([1,2,3]) #shape : (3, 13044, 11137) for reference
# make a slice of 256x256 surrounding the point
rgb_img_slice = rgb_img[:, from_row:to_row, from_col:to_col]
``````

Any clue how I can achieve this ?

What I need is some sort of function that can transform the projected X,Y values to pixel values.

## 1 Answer

You could use the `rasterio.transform.rowcol` function and the `rasterio.windows.Window` class to do a windowed read of your slice to avoid reading the entire dataset.

``````import pandas as pd
from pyproj import Transformer
import rasterio as rio
from rasterio.windows import Window, transform as windows_transform
from rasterio.transform import rowcol

proj = Transformer.from_crs("EPSG:4326", "EPSG:32650")
points = pd.read_csv("coords.csv")

# Lets just consider one image for now
with rio.open("raster.tif") as dataset:
dataset_transform = dataset.transform
profile = dataset.profile.copy()
profile.update(width=256, height=256, count=3)

# loop through each row of the dataframe
for i, row in points.iterrows():
# input is Y, X as epsg4326 axis order is lat, lon
# returns X, Y as that's the axis order of the output crs
px, py = proj.transform(row.GPS_Y, row.GPS_X)
row, col = rowcol(dataset_transform, [px], [py])
row, col = row[0], col[0]

window = Window(col_off=col-128, row_off=row-128, width=256, height=256)
window_transform = windows_transform(window, dataset_transform)
profile.update(transform=window_transform)

rgb_img_slice = dataset.read([1,2,3], window=window, boundless=True)

with rio.open(f"raster{i}.tif", "w", **profile) as output:
output.write(rgb_img_slice)```
``````
• Any clue how to save this as a raster. I tried the following but when I load on QGIS it doesnt align with orignal image. I get x,y coordinates as (~518877,76715) for bottom left point for the original image and (0,0) as x,y coordinates of the bottom left of the new image. This is the code I used to save the new image. Commented May 15 at 22:58
• `with rasterio.open(os.path.join(save_folder, save_fname), 'w', dtype = "uint8", count = 3, width = 256, height = 256, crs = "epsg:32650") as new_dataset: new_dataset.write_band([1,2,3], rgb_slice)` Commented May 15 at 22:59
• See edit - Use `rasterio.windows.transform` to get the georeferencing of the window. Note use of `boundless=True` in `dataset.read` to ensure a full 256x256 tile even near the edge of the raster. Commented May 15 at 23:50
• Yep works well, thanks heaps for this :) Commented May 16 at 4:49