topic Re: Randomly moving, and rotating, polygons in Geoprocessing Questions

Randomly moving, and rotating, polygons

MikeKnowles — Wed, 29 Feb 2012 20:33:59 GMT

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

Re: Randomly moving, and rotating, polygons

ChrisSnyder — Wed, 29 Feb 2012 22:35:17 GMT

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

Re: Randomly moving, and rotating, polygons

DavidStrip — Mon, 05 Mar 2012 03:23:45 GMT

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.

Re: Randomly moving, and rotating, polygons

MelanieHamel — Wed, 08 Feb 2017 06:01:10 GMT

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

Re: Randomly moving, and rotating, polygons

DarrenWiens2 — Sun, 12 Dec 2021 00:42:20 GMT

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

Re: Randomly moving, and rotating, polygons

MelanieHamel — Wed, 08 Feb 2017 22:07:47 GMT

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

Re: Randomly moving, and rotating, polygons

DarrenWiens2 — Wed, 08 Feb 2017 22:21:46 GMT

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

Re: Randomly moving, and rotating, polygons

MelanieHamel — Wed, 08 Feb 2017 22:31:02 GMT

Thanks again then

Re: Randomly moving, and rotating, polygons

MelanieHamel — Thu, 09 Feb 2017 23:33:51 GMT

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!

Re: Randomly moving, and rotating, polygons

PaulJensen — Fri, 19 May 2017 20:58:14 GMT

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!

Re: Randomly moving, and rotating, polygons

PaulJensen — Fri, 19 May 2017 21:00:15 GMT

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

Re: Randomly moving, and rotating, polygons

Yin-hsuenChen — Sun, 12 Dec 2021 00:42:23 GMT

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"

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Randomly moving, and rotating, polygons

SimonAllard2 — Sat, 29 Sep 2018 22:22:30 GMT

This is so much easier to do using FME...

Random rotation of polygons from a set of anchor points. - FME Knowledge Center

Re: Randomly moving, and rotating, polygons

BernardOlivier — Tue, 25 Feb 2020 11:20:27 GMT

Thanks both Yin-hsuen Chen‌ and @Darren for the code.

I am using this same code for a very large amount of polygons, and it is still working great. However, as expected, it is taking very long. Do you have any suggestions that might speed up the execution?

Thank you!