5

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
  • 1
    What does this have to do with openstreetmap? – BradHards Aug 21 '16 at 11:20
8

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
3

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:

enter image description here

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:

enter image description here

  • 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

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.