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.
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
sys.stdout.write("\r")
sys.stdout.flush()
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):
pass
else:
# 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):
found_points.append(random_point)
random_point_count += 1
break
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)
feat.SetGeometry(pt)
output_layer.CreateFeature(feat)
output_layer = None
print "\n--> finished"
