1

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')

3 Answers 3

2

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

0

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.

1
  • 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
    – g123456k
    Aug 24, 2022 at 8:39
-1

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

1
  • 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.
    – Community Bot
    Aug 24, 2022 at 7:58

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.