Randomly moving, and rotating, polygons

02-29-2012 12:33 PM
New Contributor
We have polygons for protected areas which are entirely contained within a study area polygon.  We would like to randomly move, and rotate, each protected areas polygon such that they do not overlap, and are fully contained, within the study area polygon.  We would like to repeat this randomization process x number of times.

We currently have a corridor network based on the the existing configuration of protected areas and would like to compare this corridor network with the random protected areas configurations.

Is this possible?  This is currently outside my expertise (but willing to learn) and am looking for direction/terms/concepts/etc. I should be aware of to get me started.

0 Kudos
13 Replies
Regular Contributor III
Anything is possible... You could conceivably accomplish this via Python scripting.

Feature Shifting: I don't think there is an out of the box way to do this, although you could use a "false easting" and "false northing" in a custom projection file to do this. Also, I have seen scripts posted on the forum that shift features by a systematic offset (but can't seem to locate them right now). The basic idea is to read through all the feature geometry and then write out some new features where all the vertices have a shift in their x/y coordinates (say like 20 feet to the north, 37 feet to the east).

Feature Rotating: Again, I don't think there is an out of the box tool to do this (other than the rotate tool available via an edit session). You could write some Python code to do this... although the math for doing that is above my immediate geometry skills for sure... There is a toolbox tool for rotating rasters: http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00170000007s000000 and I think there is some ArcObjects-type methods for handling vector rotation (like the edit sesion tool). This might be kind of hard in Python - and maybe better for ArcObjects if indeed there is already a method to do feature rotation there.

So I assume this is some sort of monte carlo optimization you are trying to devise? http://en.wikipedia.org/wiki/Monte_Carlo_method

I built a version of a monte carlo-based "optimizer" for estimating the number of animal territories (northern spotted owls) that could be supported given different scenarios about forest conditions over time. Basically (like what you are trying to do), it just tries a bunch of random permutations and builds up distribution of the results. It doesn't really seek to optimize anything exactly, but rather it generates a distribution of "likely outcomes". Some outcomes of course ended up with more territories than others - I guess you could think of more territories as being an optimization, although that's not what the project was really about. In retrospect, it was a pain and took me a long time to program, but in the end it was pretty cool and I learned a lot from doing it.

That said, I don't think you would be able to effectively do this kind of thing without knowing how to program. As far as I know, there are no "out of the box" tools for this sort of spatial optimization. Network Analyst is really a tabular optimization (although there is a spatial component sort of) and uses this sort of method I think: http://en.wikipedia.org/wiki/Linear_programming
0 Kudos
New Contributor II
The difficulty of your task depends a lot on how densely the protected area polygons fill the study area polygon. A second factor is whether you are looking to just perturb the protected area polygons or reshuffle them completely within the study area. As already noted, if you can write Python code, you can certainly do the random move and rotation steps - the math is pretty straightforward. However, the packing problem is extremely difficult. If your perturbations are small and the packing density is low, you can just generate a random trial, see if it satisfies the intersection/containment criteria and just throw away the ones that don't meet the rules. You can afford this strategy for sparse cases. Dense cases you'll throw away just about every trial you generate. Overall I'd say you're working on a pretty hard problem.
0 Kudos
New Contributor II

Hi Mike,

I know this is a long call since the subject is old but... I am trying to do something very similar and I'm wondering if you ever wrote that script and would like to share it? Was this work published?

Thanks in advance


MVP Honored Contributor

Fun Python task! I think it should go something like this:

import random
pa = 'protected_areas' # protected areas
sr = arcpy.Describe(pa).spatialReference # spatial ref
sa = 'study_area' # study area
extent = arcpy.Describe(sa).extent # study area extent
sa_geom = [row[0] for row in arcpy.da.SearchCursor(sa,'SHAPE@',spatial_reference=sr)][0] # study area geometry
new_polys = [] # placeholder
def rotate(poly,centroid): # magic function
    ang = random.random() * 2 * math.pi # random rotation angle in radians
    new_array = arcpy.Array() # placeholder
    rnd_x = random.uniform(extent.XMin,extent.XMax) # random x coordinate for new centroid
    rnd_y = random.uniform(extent.YMin,extent.YMax) # random y coordinate for new centroid
    rnd_centroid = arcpy.Point(rnd_x,rnd_y) # centroid point
    for part in poly: # for each polygon part
        for pnt in part: # for each vertex
            x_trans = pnt.X - centroid.X # normalize to zero
            y_trans = pnt.Y - centroid.Y # normalize to zero
            x_transprime = (math.cos(ang) * x_trans) - (math.sin(ang) * y_trans) # new x coord, from zero
            y_transprime = (math.sin(ang) * x_trans) + (math.cos(ang) * y_trans) # new y coord, from zero
            x_prime = x_transprime + rnd_centroid.X # move to new centroid x
            y_prime = y_transprime + rnd_centroid.Y # move to new centroid y
            new_array.add(arcpy.Point(x_prime, y_prime)) # add to array
    candidate_poly = arcpy.Polygon(new_array,sr) # make polygon from array
    for new_poly in new_polys: # compare to polygons already added to new_polys
        if not candidate_poly.disjoint(new_poly) or not sa_geom.contains(candidate_poly): # check if overlap with polygons and is fully inside study area
            candidate_poly = rotate(poly,centroid) # if the new polygon violates criteria, call function again
    return (candidate_poly) # return the final good polygon
with arcpy.da.SearchCursor(pa,'SHAPE@',spatial_reference=sr) as cursor: # for each polygon
    for row in cursor:
        centroid = row[0].centroid # calculate centroid
        new_polys.append(rotate(row[0],centroid)) # add returned good polygon to list
arcpy.CopyFeatures_management(new_polys,r'in_memory\new_polys') # write to disk

* rotation math from here: Rotating Features in ArcGIS for Desktop using ArcPy? - Geographic Information Systems Stack Exchange 

New Contributor II

Oh wow, thank you very much for this (and so quick), I'm looking forward to try it! May I ask what you used it for (and if the work can be read somewhere - just curious)?

0 Kudos
MVP Honored Contributor

Errr... I used it for making those example images.  I just made it up. This is all there is to it!

New Contributor II

Thanks again then 

0 Kudos
New Contributor II

So I think I understand your script and that kind of works thanks - but...

1. It looks like the small polygons (protected areas) need to be single part for it to work, even though it looks like the loop does look at multipart (?). I tried with sample data - did not work with multipart but worked fine with single.

2. My new random polygons always fall within the extent of the study area, but not always within the actual study area (which is a complex multipart polygon too). Although only using a large single part study area still did not work.

Note that I'm a Python newbie (but work well in R) so obvious things to you might not be obvious to me

I'll then have to loop that randomisation n times and calculate values based on intersection between new polygons and other layers for every randomisation... steep learning curve but will get there!

0 Kudos
New Contributor

Did you ever figure out your second issue?  I'm having the same problem.  Thanks!

0 Kudos