2

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 1

3

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)```
4
  • 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.
    – user144483
    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)
    – user144483
    Commented May 15 at 22:59
  • 1
    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.
    – user2856
    Commented May 15 at 23:50
  • Yep works well, thanks heaps for this :)
    – user144483
    Commented May 16 at 4:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.