Find raster pixel coordinates on shapefile polygon corners using Python

Question

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.

enter image description here

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
from rasterio.mask import mask
import geopandas as gpd
shapefile = gpd.read_file(r'W:\shape.shp')
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')

Answer 1

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
 
shapefile = gpd.read_file('your_shape.shp')
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
 """
 # Read shapefile
 gdf = gpd.read_file(shapefile)
 
 # 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"
 
 # add filename
 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)

enter image description here

Answer 2

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.

• g123456k: 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

Answer 3

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

• Community: Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.