Randomly moving, and rotating, polygons

Discussion created by msknowles on Feb 29, 2012
Latest reply on Sep 29, 2018 by simonal_boprc

Like • Show 0 Likes0
Comment • 12

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.

Thanks,
Mike

Outcomes

12 Replies

csny490 Feb 29, 2012 2:35 PM
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
Like • Show 0 Likes0
Actions
drstrip Mar 4, 2012 7:23 PM
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.
Like • Show 0 Likes0
Actions
Lepidochelys75 Feb 7, 2017 10:01 PM
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
Mel
Like • Show 1 Like1
Actions
- dkwiens @ Lepidochelys75 on Feb 8, 2017 11:48 AM
  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
  3 people found this helpful
  Like • Show 4 Likes4
  Actions
  - Lepidochelys75 @ dkwiens on Feb 8, 2017 2:07 PM
    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)?
    Like • Show 0 Likes0
    Actions
    - dkwiens @ Lepidochelys75 on Feb 8, 2017 2:21 PM
      Errr... I used it for making those example images. I just made it up. This is all there is to it!
      Like • Show 1 Like1
      Actions
      - Lepidochelys75 @ dkwiens on Feb 8, 2017 2:31 PM
        Thanks again then
        Like • Show 0 Likes0
        Actions
  - Lepidochelys75 @ dkwiens on Feb 9, 2017 3:33 PM
    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!
    Like • Show 0 Likes0
    Actions
    - pgjensen99 @ Lepidochelys75 on May 19, 2017 2:00 PM
      Did you ever figure out your second issue? I'm having the same problem. Thanks!
      Like • Show 0 Likes0
      Actions
  - pgjensen99 @ dkwiens on May 19, 2017 1:58 PM
    Thank you for your script...huge help for someone that has never used Python! As Melanie mentioned below, the script mostly works but random polygons are only completely within my defined study area polygon about 30-40% of the time. For those that aren't, at least some part of the random polygon is within the study area. I could generate a large sample of the random polygons and simply discard those that aren't completely contained, but do you see anything obvious in your original script that would explain what's going on. Again, many thanks!
    Like • Show 0 Likes0
    Actions
  - simonal_boprc @ dkwiens on Sep 29, 2018 3:22 PM
    This is so much easier to do using FME...
    
    Random rotation of polygons from a set of anchor points. - FME Knowledge Center
    Like • Show 0 Likes0
    Actions

eisen520 Aug 19, 2018 11:12 AM

Hey all,

Thanks @Darren Wiens for this really cool code. I was helping my friend figuring out this code. The biggest issue to me is that the random polygon would not fall in the study area polygon. So I added a try/except and if/else statement to check and return until every polygon fall in the study area. Again, I have a simply polygons for test so it might not working if you have huge amount of polygon, and these polygons should be singlepart and no hollow inside. Also, I remove the "disjoint" function for checking if the result polygon is overlap or not. I hope this help

import arcpy
from arcpy import da
import random
import math

#set-up environment
workspace = r"...:/.../..." #location of your files
pa = workspace + '\\' + "polygons.shp"  # polygons to randomly rotate and move
sr = arcpy.Describe(pa).spatialReference # spatial reference
sa = workspace + '\\' + "study_area.shp" # 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

# magic function
def rotate(poly,centroid): 
    #ang = random.random() * 2 * math.pi # random rotation angle in radians if you want to have unified angle to result polygons 
    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) # random centroid point

    while True:
            try:
                for part in poly:#for each polygon part
                    ang = random.random() * 2 * math.pi # random rotation angle in radians # set as 0 if you don't need to rotate polygon
                
                    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

                temp_poly = arcpy.Polygon(new_array,sr) # make a temporarily polygon for test

                if sa_geom.contains(temp_poly) is True:
                    print 'inside'
                    break             
                else:
                    raise Exception
            except: 
                print 'outside sa, start over...'
                return rotate(poly,centroid) #here you start over and run the loop again... until all the polygon fall in study area.

            break
                
    good_poly = arcpy.Polygon(new_array,sr) # make all good polygons (indside study area) from array
    return (good_poly) # return the final good polygon

print "start processing..."

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, workspace + '\\' + "TheResult.shp") # write to disk


print "done"

1 person found this helpful

Actions

GeoNet

The Esri Community

Randomly moving, and rotating, polygons

Outcomes

Related Content

Recommended Content