# Sample all existing points of a polygon in Python

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
``````
• 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. Aug 12, 2020 at 12:14
• mask then read masked values - rasterio.readthedocs.io/en/latest/topics/… Aug 13, 2020 at 7:25
• Are you ultimately looking to calculate raster statistics within each polygon?
– Aaron
Jan 23, 2021 at 15:24

## 1 Answer

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
``````