Generating random coordinates in multipolygon in Python?

Question

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:
    data = json.load(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:
    data = json.load(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

Answer 1

With Shapely: https://shapely.readthedocs.io/en/latest/manual.html#polygons

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

Answer 2

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:
    data = json.load(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)
        point.AddPoint(random.uniform(xmin, xmax),
                       random.uniform(ymin, ymax))

        if point.Within(polygon):

            multipoint.AddGeometry(point)

            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:

Answer 3

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.

Answer 4

import geopandas as gd
from shapely.geometry import Point
import random

plgeo = gd.read_file('POL_2.gpkg')
plgeo = plgeo.sort_values('Population', ascending=False)
top_regions = plgeo.head(20)

def rand_range(a, b, cnt):
    return [random.uniform(a, b) for _ in range(cnt)] 

def gen_points(row):
    max_cnt = int(row.Population / 1000)
    [xmin, ymin, xmax, ymax] = row.geometry.bounds 
    bbox-points = zip(rand_range(xmin, xmax, max_cnt), rand_range(ymin, ymax, max_cnt))  
    points = [p for p in bbox-points if row.geometry.contains(Point(p))]
    return {'id': row.HASC_2, 'name': row.NAME_2, 'points': points}

with_points = top_regions.apply(gen_points, axis=1)