# Generating random coordinates in 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, env, env, env
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
``````

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

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

def generate_random(number, polygon):
points = []
minx, miny, maxx, maxy = polygon.bounds
while len(points) < number:
pnt = Point(random.uniform(minx, maxx), random.uniform(miny, maxy))
if polygon.contains(pnt):
points.append(pnt)
return 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,env,env,env

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

counter += 1

print multipoint.ExportToWkt()
``````

I tried out above code with only one feature ('name':'A') of next "GeoJSON" vector layer for generating 1000 random points: After running the code, I used 'QuickWKT' plugin of QGIS for displaying the WKT multi point format printed at the Python Console of QGIS. Result was:

• Thank you so much . You can see the result HERE – Sadeq Sepehrnoush Aug 25 '16 at 7:13
• Careful if you use this in WGS84 or something similar. For Shapely there is no globe or sphere, all coordinates and coordinate ranges are equal to it. For geodesic uniformity you need to use something more sophisticated. For Shapely, you must use a projected coordinate system where distances are equal anywhere in it. – bugmenot123 Nov 20 '18 at 12:05

I had a similar use case and since I did not find any library available, I developed this python library: PyCristoforo. Github link: https://github.com/AleNegrini/PyCristoforo

Version 1.0.0 only supports European countries, but I plan to release other countries soon.