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
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
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)?
Errr... I used it for making those example images. I just made it up. This is all there is to it!
Thanks again then
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!
Did you ever figure out your second issue? I'm having the same problem. Thanks!