I think my question is answered in different posts but I can't really make sense of them.

Background: We are designing a survey of a sample population of farmers to identify the seeds they use. I have a vector layer indicating the land use (irrigated or rainfed). This layer is composed by 16000 polygons spread across the country.

What I need: I want to randomly generate points within this polygons where to survey. More specifically, I want to generate 66 points in every province (provinces are indicated by another vector layer). In other words I want to generate 66 points under two conditions: 1) within province boundaries, 2) within the land use polygons.

I already tried QGIS plugin "Random points" but it seems to crash the program...

How should I proceed?


2 Answers 2


How many provinces are in your dataset? I threw together a python script wich handles your task. However, it might be too slow. Below you can see a sample output with mockup data. The green polygons are farmland, the black lines indicate the border of a province, the points are what you want.

enter image description here

I wonder if your criteria should also be that the points are not too close together. Also, you should keep in mind that, depending on province size and percentage of farmland per province, the point density will vary (as you can see in the screenshot).

Maybe a better idea would be to create a point layer with equal distance between each point (think of grid square centroids), then sample from these via a simple random pick.

The python script (which is not optimized at all, just quickly thrown together and possibly too slow for your use case! Uses shapely and ogr):

    from shapely import wkb, geometry
import random
import ogr
import sys

# setup shapefile driver
driver = ogr.GetDriverByName('ESRI Shapefile')

# import shapefiles
sample_provinces_shapefile = driver.Open("sample_provinces.shp")
sample_farmplots_shapefile = driver.Open("sample_farmland.shp")

# get layers
sample_provinces = sample_provinces_shapefile.GetLayer()
sample_farmplots = sample_farmplots_shapefile.GetLayer()
number_of_provinces = sample_provinces.GetFeatureCount()

# array for found points
found_points = []

# create random points for each province, points have to be inside farm plots
for i in xrange(number_of_provinces):

    # convert to shapely
    province_tmp = sample_provinces.GetFeature(i)
    province = wkb.loads(province_tmp.GetGeometryRef().ExportToWkb())

    # reset random point count
    random_point_count = 0

    # visual output
    print("processed "+str(i)+" of "+str(number_of_provinces)+" "),

    # fill random points
    while random_point_count < 66:
        # create random point from polygon boundary
        maxx, maxy, minx, miny = province.bounds
        random_point = geometry.Point(minx + (random.random() * (maxx-minx)), (miny + (random.random() * (maxy-miny))))

        # check whether it is inside the polygon
        if not province.contains(random_point):
            # check whether it is inside a farm plot (brute force, will be slow with large datasets)
            # while loop to prevent segmentation fault
            z = 0
            while z != sample_farmplots.GetFeatureCount():
                plot_tmp = sample_farmplots.GetFeature(z)
                plot = wkb.loads(plot_tmp.GetGeometryRef().ExportToWkb())

                # first, compare to bounding box, then to actual geometry (speeds things up a little)
                if (plot.bounds[0] < random_point.x < plot.bounds[2]) and (plot.bounds[1] < random_point.y < plot.bounds[3]):
                    if plot.contains(random_point):
                        random_point_count += 1
                z += 1

# output shapefile with points
output_file = driver.CreateDataSource("random_points.shp")
output_layer = output_file.CreateLayer("point_out", None, ogr.wkbPoint )
for point in found_points:
    feat = ogr.Feature(output_layer.GetLayerDefn())
    pt = ogr.Geometry(ogr.wkbPoint)
    pt.SetPoint_2D(0, point.x, point.y)

output_layer = None
print "\n--> finished"

Many frameworks can give you a point on a surface (i.e. polygon). For example, you can do this with SpatiaLite (using SELECT PointOnSurface(polygon blob)), which is a wrapper around GEOSPointOnSurface(). Other packages such as postgis, and GDAL/OGR also use the GEOS code.

I would caution that I'm not at all sure how "random" they will be for multiple calls though, which might skew the results.

  • Thanks. So do you suggest to drop QGIS to do this?
    – Francesco
    Aug 11, 2012 at 0:14
  • I'm not sure how QGIS does it either. Perhaps a query to their mailing list (qgis.org/community/mailing-lists.html) would be a better place than this kind of forum.
    – BradHards
    Aug 11, 2012 at 0:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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