Random polygons
2849
5
01-17-2018 10:37 AM
by
by
by
XanderBakker
I thought I would have a go at it, but there are a number of things that would need to be resolved.
- rectangles should not be reused
- rectangles have common origin XY
- minimum area of 25% should be covered
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.
