I'm trying to sample raster values to points based on the geometries of a shapefile using the below code but the raster values are only sampled at the coordinates of the polygon geometry. Is there a way to sample for all the points along the boundary and inside the polygon.

multi_values_points = pd.Series()
 with fiona.open(inputShape) as shp:
     for feature in shp:
        Name = feature['properties']['Name']
        coords = feature['geometry']['coordinates']
        with rio.open(inputRaster) as stack_src:
           value = [v for v in stack_src.sample([coords])]
        multi_values_points.loc[Name] = value
  • 1
    By definition, there are an infinite number of points inside a polygon. In the past, when I needed to make a heat map from a polygon, I 1) chose a grid size, 2) generated a multipoint from the centers of the hexagon grid covering the extent of the polygon, 3) applied a negative buffer of sqrt(2)/2.0 times the grid spacing to the source polygon, 4) intersected the points with the buffered geom 5) densified the boundary of the original polygon to three times the grid size, and 6) appended the densified boundary points to the interior point sample. The result looked pretty good.
    – Vince
    Commented Aug 12, 2020 at 12:14
  • 1
    mask then read masked values - rasterio.readthedocs.io/en/latest/topics/…
    – user2856
    Commented Aug 13, 2020 at 7:25
  • Are you ultimately looking to calculate raster statistics within each polygon?
    – Aaron
    Commented Jan 23, 2021 at 15:24

1 Answer 1


I tried to generate points inside my various polygon boundaries of the shapefile by using below code:

import ogr
import random
from osgeo import ogr
from shapely.geometry import Point
from shapely.prepared import prep
from shapely.wkb import loads

g = ogr.Open(inputShapeFile)# Open the shapefile
for layer in g.GetLayer(0):
    field = loads(layer.GetGeometryRef().ExportToWkb())
    fieldP = prep(field)
    lon_min, lat_min, lon_max, lat_max = field.bounds
    n_points = 0 # A counter
    field_lonlat = []
    while n_points < 10:
        r_lon = random.uniform(lon_min, lon_max)
        r_lat = random.uniform(lat_min, lat_max)
        pt = Point(r_lon, r_lat)
        if fieldP.contains(pt):
            field_lonlat.append([r_lon, r_lat])
            n_points += 1

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.