Find raster pixel coordinates on shapefile polygon corners using Python

I have a raster file and shapefile with multiple polygons.

How do I find raster pixel coordinates on shapefile polygon corners?

Not the geographical coordinates based on CRS, but the pixel row and column values on polygon corner location.

Later I will need to use these pixel coordinates to clip the raster image.

I tried something like this, but extracting all the px as points is not really a solution as it takes a lot of time.

And I don't believe that this is even the right way how to solve this.

``````import rasterio
import geopandas as gpd

geoms = shapefile.geometry.values
geometry = geoms[0]

from shapely.geometry import mapping
geoms = [mapping(geoms[0])]

with rasterio.open(r'W:\raster.tif') as src:
out_image, out_transform = mask(src, geoms, crop=True)

no_data=src.nodata
data = out_image[0,:,:]

for idx, row in shapefile[0:1].iterrows():
row, col = np.where(data != no_data)
elev = np.extract(data != no_data, data)
from rasterio import Affine # or from affine import Affine
T1 = out_transform * Affine.translation(0.5, 0.5) # reference the pixel centre
rc2xy = lambda r, c: (c, r) * T1

d = gpd.GeoDataFrame({'col':col,'row':row,'elev':elev})
# coordinate transformation
d['x'] = d.apply(lambda row: rc2xy(row.row,row.col)[0], axis=1)
d['y'] = d.apply(lambda row: rc2xy(row.row,row.col)[1], axis=1)
# geometry

from shapely.geometry import Point
d['geometry'] =d.apply(lambda row: Point(row['x'], row['y']), axis=1)
d = d.set_crs('epsg:32635')

point_with_polygons = gpd.sjoin(left_df=d, right_df=shapefile, how='inner')
``````

Well after some digging I have found some solution with a help from https://gis.stackexchange.com/a/409625/178904:

``````import os
import geopandas as gpd
import pandas as pd
import rasterio
import shapely

raster_path = 'your_shape.tif

# This will give us a dataframe of pixel coordinates according to your shape file

def shapefile_to_annotations(shapefile, rgb):
"""
Args:
shapefile: Path to a shapefile on disk. If a label column is present, it will be used, else all labels are assumed to be "Tree"
rgb: Path to the RGB image on disk
Returns:
results: a pandas dataframe
"""

# get coordinates
df = gdf.geometry.bounds

# raster bounds
with rasterio.open(rgb) as src:
left, bottom, right, top = src.bounds
resolution = src.res[0]

# Transform project coordinates to image coordinates
df["tile_xmin"] = (df.minx - left) / resolution
df["tile_xmin"] = df["tile_xmin"].astype(int)

df["tile_xmax"] = (df.maxx - left) / resolution
df["tile_xmax"] = df["tile_xmax"].astype(int)

# UTM is given from the top, but origin of an image is top left

df["tile_ymax"] = (top - df.miny) / resolution
df["tile_ymax"] = df["tile_ymax"].astype(int)

df["tile_ymin"] = (top - df.maxy) / resolution
df["tile_ymin"] = df["tile_ymin"].astype(int)

# Add labels is they exist
if "label" in gdf.columns:
df["label"] = gdf["label"]
else:
df["label"] = "Tree"

df["image_path"] = os.path.basename(rgb)

# select columns
result = df[[
"image_path", "tile_xmin", "tile_ymin", "tile_xmax", "tile_ymax", "label"
]]
result = result.rename(columns={
"tile_xmin": "xmin",
"tile_ymin": "ymin",
"tile_xmax": "xmax",
"tile_ymax": "ymax"
})

# ensure no zero area polygons due to rounding to pixel size
result = result[~(result.xmin == result.xmax)]
result = result[~(result.ymin == result.ymax)]

return result

results=shapefile_to_annotations(shapefile, raster_path)

# print(result)
#   xmin    ymin    xmax    ymax
# 0 4023    5033    4288    5398
# 1 3966    4880    4331    5027
# 2 3146    4567    3398    4879

from osgeo import gdal

"""
Now, we can use 'results' pixel coordinates to clip some part of whole raster image
srcWin = [xmin, ymin, (xmax-xmin), (ymax-ymin)]

"""

out = 'out_raster.tif'
ds = gdal.Open(raster_path)
ds = gdal.Translate(out, ds, srcWin = [5000,4523, 3500, 1452])
ds = None
``````

And the result with blue as an original raster image (`raster_path`) and green as clipped one (`ds`)

Are you looking for something like this (i can not comment yet, thus as answer...): Extract corner coordinates of polygon extent

I had a similar issue earlier on and that solved it for me.

• I'm not sure this solution gives me raster pixel index values for the polygon corners. I see that they simply extracting polygon bounding box coordinates and not raster pixel coordinates Aug 24, 2022 at 8:39

You can get the points coordinates of polygon and then extract raster values of these points. see this link Extract corner coordinates of polygon extent in PyQGIS / GDAL