Select to view content in your preferred language

Random polygons

2761
5
01-17-2018 10:37 AM
FernandoCoelho
New Contributor

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?

0 Kudos
5 Replies
DanPatterson_Retired
MVP Emeritus

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)

FernandoCoelho
New Contributor

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"?

0 Kudos
DanPatterson_Retired
MVP Emeritus

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.

0 Kudos
XanderBakker
Esri Esteemed Contributor

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.

DanPatterson_Retired
MVP Emeritus

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