# How to generate random coordinates in a multipolygon in python

How to generate random points in multi-polygon using geojson in python ? thus far I found a class in JavaScript named :

``````randomPointsOnPolygon(numberOfPoints, polygon)
``````

but I need to use something like that in python , my code is like

``````import json
import geojson

with open('my_city_boundaries.geojson') as f:

for feature in data['features']:
if feature['properties']["name:en"] == state_name:
# I need to generate some coordinates within that state
``````

I tried shapely like this and it didnt work . it doesnt recognize Points attribute

``````import shapely import Point
import json
import geojson
from osgeo import ogr

def generate_random(number, polygon):
list_of_points = []
minx, miny, maxx, maxy = env[0], env[2], env[1], env[3]
counter = 0
while counter < number:
pnt = Point(random.uniform(minx, maxx), random.uniform(miny, maxy))
if polygon.contains(pnt):
list_of_points.append(pnt)
counter += 1
return list_of_points

with open('ir_states_boundaries_coordinates.geojson') as f:
print("\n\nName your state from this list\n\n")
for feature in data['features']:  # print a list of valid state names
print(feature['properties']['name:en'])
state_name = raw_input("\n")
for feature in data['features']:
if feature['properties']["name:en"] == state_name:
geom = feature['geometry']
geom = json.dumps(geom)
polygon = ogr.CreateGeometryFromJson(geom)
env = polygon.GetEnvelope()
result = generate_random(10, polygon)
print result
``````
• What does this have to do with openstreetmap? – BradHards Aug 21 '16 at 11:20

With Shapely: http://toblerity.org/shapely/manual.html#polygons

The `polygon` in code below must be Polygon object.

``````import random
from shapely.geometry import Point

def generate_random(number, polygon):
list_of_points = []
minx, miny, maxx, maxy = polygon.bounds
counter = 0
while counter < number:
pnt = Point(random.uniform(minx, maxx), random.uniform(miny, maxy))
if polygon.contains(pnt):
list_of_points.append(pnt)
counter += 1
return list_of_points
``````
• i think it must be like this `from shapely import geometry` then use it like `geometry.Point` – Sadeq Sepehrnoush Aug 25 '16 at 7:11
• Note: These will not be randomly distributed! More points will be clustered towards the poles with the effect worse the closer you are to the poles. – Cramer Aug 6 '18 at 23:54
• @Cramer: They ARE randomly distributed in a uniform matter on the surface of the polygon. Shapely operates on a 2D euclidean coordinate system. If you need geodesic uniformity, you need to use something else. – bugmenot123 Nov 20 '18 at 12:03
• @bugmenot123 The data is coming from geojson, which uses coordinates. The space of coordinates is two dimensional but it's not Euclidean, even if shapely believes it is. The data I've used shapely for was in lon/lat, the edges were relatively short so treating it as Euclidean was a good approximation. BUT, if you're generating random points you need to consider the shape of the space. If the points are all within a few hundred km, the Euclidean approximation is fine, if you're generating points in the shape of Africa, you're going to have some weird effects. – Cramer Nov 21 '18 at 13:16
• We agree on the problem but I just want to stress the point. Shapely considers whatever coordinates to throw at it to be in 2D euclidean space. If you feed it coordinates that you know are geographic, it will be a plate carrée projection. Shapely will do all its calculation like that and it is correct because Shapely does not care about coordinate systems. If you need geodesic calculations (like distributing points uniformly on a sphere/ellipsoid, you need to use something else (geographiclib?) or implement the calculations yourself. – bugmenot123 Nov 21 '18 at 14:32

Next, I used the first part of your code with 'ogr' module; although this one has the possibility of "GeoJSON" driver. However, when it's used 'json' module, it is necessary to create a geometry for selected feature. Once created, you can get its bounding box to narrow generated random points. The 'Within' method corroborates if each point is into the feature.

``````import json
import geojson

import random

from osgeo import ogr

with open('/home/zeito/pyqgis_data/xwRcl.geojson') as f:

for feature in data['features']:
if feature['properties']['name'] == 'A':

geom = feature['geometry']
geom = json.dumps(geom)
polygon = ogr.CreateGeometryFromJson(geom)

env = polygon.GetEnvelope()
xmin, ymin, xmax, ymax = env[0],env[2],env[1],env[3]

num_points = 1000

counter = 0

multipoint = ogr.Geometry(ogr.wkbMultiPoint)

for i in range(num_points):
while counter < num_points:

point = ogr.Geometry(ogr.wkbPoint)
random.uniform(ymin, ymax))

if point.Within(polygon):