Hello guys.
I need to create 20X30m polygons at random positions within a predefined boundary. And I also need the number of polygons created equal 25% of the area.
Any suggestion?
you could generate a fishnet within the area that more or less covers it. Assign a random number to each grid square then sample 25% and dump the rest.
Also, random polygons could also be accomplished by buffering 25 random points that are at a distance no-closer-than (I have code for this in python)
Do you want random shapes? (not going to happen easily)
I'd like the shape of the polygons to be rectangular (20 X 30 meters). And, where would the "fishnet" look for the area of the "base polygon" to know how many polygons (20X30m) to be 25% sampled area of the "base polygon"?
That might entail 'humanware' to provide some guidance. You could however, determine the extent of the base polygon then use that extent to set the bounds of the fishnet. You would just want to ensure that your fishnet has sufficient rows and columns that you could get a 25% sample of rectangles within the base polygon. That is not a problem since 'select by spatial location' with a 'completely within' control would ensure that you select only polygons within the base polygon, then continue to select polygons perhaps incrementally, until the 25% areal coverage is met.
I thought I would have a go at it, but there are a number of things that would need to be resolved.
This is the code I came up with, but it is far from fast and should not be used for a "large" set of data:
import arcpy
def main():
arcpy.env.overwriteOutput = True
# settings
fc_in = r'C:\GeoNet\RandRectangles\data.gdb\areas'
fc_tmp = r'C:\GeoNet\RandRectangles\data.gdb\rectangles_tmp'
fc_sel = r'C:\GeoNet\RandRectangles\data.gdb\rectangles_tmp_sel'
fc_out = r'C:\GeoNet\RandRectangles\data.gdb\rectangles_v01'
rect_width = 20
rect_height = 30
min_area_frac = 0.25
# Some asumptions:
# rectangles should not be reused
# rectangles have common origin XY
# minimum area of 25% should be covered
sr = arcpy.Describe(fc_in).spatialReference
ext = arcpy.Describe(fc_in).extent
# create a tmp featureclass with the rectangles
ext_fin = ReshapeExtent(ext, rect_width, rect_height)
feats = []
for c in range(int(ext_fin.width / rect_width)):
for r in range(int(ext_fin.height / rect_height)):
polygon = CreateRectanglePolygon(ext_fin.XMin + c * rect_width,
ext_fin.YMin + r * rect_height,
ext_fin.XMin + (c+1) * rect_width,
ext_fin.YMin + (r+1) * rect_height,
sr)
feats.append(polygon)
arcpy.CopyFeatures_management(feats, fc_tmp)
# now select only the ones that overlap any area polygon
arcpy.MakeFeatureLayer_management(fc_in, 'area_lyr')
arcpy.MakeFeatureLayer_management(fc_tmp, 'rect_lyr')
arcpy.SelectLayerByLocation_management("rect_lyr", "INTERSECT", "area_lyr", None, "NEW_SELECTION")
arcpy.CopyFeatures_management('rect_lyr', fc_sel)
# create a dictionary from the selected reactangles
dct_rect = {r[0]: r[1] for r in arcpy.da.SearchCursor(fc_sel, ('OID@', 'SHAPE@'))}
# loop through polygons and determine which to use
lst_oid_taken = []
lst_results = []
with arcpy.da.SearchCursor(fc_in, ('OID@', 'SHAPE@')) as curs:
for row in curs:
oid_area = row[0]
area_pol = row[1]
dct_overlap = GetDictOfRectanglesForArea(area_pol, dct_rect, lst_oid_taken)
lst_selection, lst_oid_taken = GetRandomRectanglesForArea(area_pol, dct_overlap, lst_oid_taken, min_area_frac)
lst_results.append([oid_area, lst_selection])
for result in lst_results:
oid_area = result[0]
lst_selection = result[1]
for oid_rect in lst_selection:
print "{0}\t{1}".format(oid_rect, oid_area)
def GetRandomRectanglesForArea(area, dct_overlap, lst_oid_taken, min_area_frac):
import random
lst_selection = []
lst_oid = dct_overlap.keys()
polygon_area = area.area
tot_rect_area = 0
fraction = tot_rect_area / polygon_area
cnt = 0
print "\nGetRandomRectanglesForArea"
while fraction < min_area_frac:
cnt += 1
print " - cnt:", cnt
i = random.randrange(0, len(lst_oid) - 1)
print " - i:", i
oid = lst_oid[i]
print" - oid:", oid
if not oid in lst_oid_taken:
rectangle = dct_overlap[oid]
overlap_area = GetOverlapArea(area, rectangle)
print " - overlap_area:", overlap_area
tot_rect_area += overlap_area
print " - tot_rect_area:", tot_rect_area
fraction = tot_rect_area / polygon_area
print " - fraction:", fraction
lst_oid_taken.append(oid)
lst_selection.append(oid)
if cnt > 100:
break
return lst_selection, lst_oid_taken
def GetOverlapArea(area, rectangle):
if rectangle.within(area):
return rectangle.area
else:
overlap = area.intersect(rectangle, 4)
return overlap.area
def GetDictOfRectanglesForArea(area, dct, lst_done):
dct_overlap = {}
for oid, rectangle in dct.items():
if not oid in lst_done:
if rectangle.within(area):
dct_overlap[oid] = rectangle
elif rectangle.overlaps(area):
dct_overlap[oid] = rectangle
return dct_overlap
def CreateRectanglePolygon(xmin, ymin, xmax, ymax, sr):
points = [arcpy.Point(xmin, ymin), arcpy.Point(xmin, ymax),
arcpy.Point(xmax, ymax), arcpy.Point(xmax, ymin),
arcpy.Point(xmin, ymin)]
return arcpy.Polygon(arcpy.Array(points), sr)
def ReshapeExtent(ext, width, height):
xmin_new = divmod(ext.XMin, width)[0] * width
ymin_new = divmod(ext.YMin, height)[0] * height
xmax_new = divmod(ext.XMax, width)[0] * width + width
ymax_new = divmod(ext.YMax, height)[0] * height + height
return arcpy.Extent(xmin_new, ymin_new, xmax_new, ymax_new)
if __name__ == '__main__':
main()
The result is list of rectangle OIDs with the corresponding area OIDs which I joined back to the rectangles (since I was too lazy to create a proper output featureclass.
The steps are more or less the following:
1) a tmp featureclass with rectangles is created using lower left XY which is a multiple of the width and height and covers the entire area:
2) The rectangle features that overlap the areas are selected and read into a dictionary:
3) random rectangles are selected until at least 25% overlap is obtained (taking only the overlapping area of the rectangles into account):
There are still some situations that are probably not desired. Such as rectangles that have a very small overlap or have a larger overlap with another polygon and that might not be what you want:
So it will probably need to be optimized further.
Nice! Now if you want Xander Bakker try option 4 (spaced points ) in my Point Tools for ArcGIS Pro if you are looking to give 2.1 a test ride. Point spacing should be guaranteed, then buffer by 1/2 the distance should also ensure no overlaps