arcgis-10.2arcpymergepolygonshapefile
I have a shapefile of polygons. I do not want any polygons with an area less than 25. For polygons with an area less than 25 I want to merge them with another polygon with an area greater than 25 that also has the same values for att1, att2 and att3.
I was able to write some code which found the polygons with areas less than 25 but I am stuck after that. can figure this out in ArcMap or ArcPy.
I want to automate this, I understand it can be done one by one but.
First of all, never use the regular cursors when you can use the data access (da) cursors introduced at 10.1. They are 10 times faster at least. Second, the Spatial Join can do the linking of polygon attributes to points faster than any series of multiple queries. A single Spatial Join and a da search cursor read into a dictionary and then transfer of the dictionary to the original points with a da update cursor will finish a quite large set of points and polygons in about 3 to 5 minutes or less. This assumes you do not want to replace your original points with a new point feature class directly created by the Spatial Join tool. It also assumes that the fsapoly field and fsapoint field are not the same field name. It they are the same name then "1" or "_1" will have to be appended to the fsapoly field name after the spatial join occurs.
This code could transfer many fields if you include the fields in the field lists that create a match between the source and the target. The more fields transferred using this code the more time it will save over other processes. For example: sourceFieldsList = ["TARGET_FID", fsapoly, "FSA_NAME"] updateFieldsList = ["OID@", fsapoint, "FSA_NAME"] will work to fill in the two fields with no other changes to the code.
The use of the da cursors and dictionary transfer data 10 times faster than an attribute join and the field calculator, even without any attribute index on either the source or the target. I have rewritten all of my scripts that used an attribute join and the field calculator to use similar code to what is written below and saved at least 2 hours of processing time so that these steps now process in about 20 to 30 minutes.
# Import arcpy module
import arcpy
# Allow overwrite
arcpy.env.overwriteOutput = True
# Script user input parameters
polygonLayer = arcpy.GetParameter(0)
pointLayer = arcpy.GetParameter(1)
workspace = arcpy.GetParameter(2)
fsaPoly = arcpy.GetParameterAsText (3) #FSA in FSA_Boundaries
fsaPnt = arcpy.GetParameterAsText (4) #FSA in Addresses
sjpoints = arcpy.GetParameter(5) # output of Spatial Join
#Run the Spatial Join tool, using the defaults for the join operation and join type
arcpy.SpatialJoin_analysis(pointlayer, polygonlayer, sjpoints)
# define the field list from the spatial join
sourceFieldsList = ["TARGET_FID", fsapoly]
# define the field list to the original points
updateFieldsList = ["OID@", fsapoint]
# populate the dictionary from the polygon
valueDict = {r[0]:(r[1:]) for r in arcpy.da.SearchCursor(sjpoints, sourceFieldsList)}
with arcpy.da.UpdateCursor(updateFC, updateFieldsList) as updateRows:
for updateRow in updateRows:
keyValue = updateRow[0]
if keyValue in valueDict:
for n in range (1,len(sourceFieldsList)):
updateRow[n] = valueDict[keyValue][n-1]
updateRows.updateRow(updateRow)
del valueDict
Use my script from Checking if polygon fits inside another polygon using ArcGIS or QGIS?
to create points inside your polygons. Each point is a centre of maximum inscribed circle of the parcel. Please note, script tested on shapefiles only and it does not handle polygons with holes, i.e. donut like polygons.
Table of centre point contains field theDist, that is the radius of above circle.
Select all points with theDist=>20, and create 20m buffers around them:
Use Minimum Bounding Geometry to create rectangles from buffers
Attach this script to tool with no parameters. It is a tool to run from mxd
import arcpy, traceback, os, sys,math
from math import radians,sin,cos
mxd = arcpy.mapping.MapDocument("CURRENT")
parent=arcpy.mapping.ListLayers(mxd,"parent")[0]
child=arcpy.mapping.ListLayers(mxd,"child")[0]
# NUMERIC ! common field
linkField="PAR_ID"
d=arcpy.Describe(parent)
SR=d.spatialReference
g=arcpy.Geometry()
gr=(math.sqrt(5)-1)/2
try:
def showPyMessage():
arcpy.AddMessage(str(time.ctime()) + " - " + message)
#golden section to find minimum
def gss(a,b,tol):
c=b-gr*(b-a)
d=a+gr*(b-a)
while abs(c-d)>tol:
fc=f(c);fd=f(d)
if fc<fd:
b=d
d=c
c=b-gr*(b-a)
else:
a=c
c=d
d=a+gr*(b-a)
return (b+a)/2
# rotate polygon
def ShapeMake(pGon,angle):
ar=arcpy.Array()
a=radians(angle)
part=pGon.getPart(0)
for p in part:
x,y=p.X-pGon.centroid.X,p.Y-pGon.centroid.Y
xN=cos(a)*x+sin(a)*y
yN=-sin(a)*x+cos(a)*y
pN=arcpy.Point(xN+pGon.centroid.X,yN+pGon.centroid.Y)
ar.add(pN)
pgonRotated=arcpy.Polygon(ar,SR)
return pgonRotated
#function to minimise
def f(a):
pgonRot=ShapeMake(square,a)
intR=pgonRot.difference(bigPolygon)
return intR.area
result=arcpy.GetCount_management(child)
nF=int(result.getOutput(0))
initFidList=[]
arcpy.SetProgressor("step", "", 0, nF,1)
with arcpy.da.UpdateCursor(child, ("SHAPE@",linkField)) as rows:
for row in rows:
square=row[0]
PID=row[1]
quer='%s%s%s=%i'%('"',linkField,'"',PID)
parent.definitionQuery=quer
bigPolygon=arcpy.CopyFeatures_management(parent,g)[0]
intR=square.difference(bigPolygon)
# square fits inside parent without rotation
if intR.area==0:initFidList.append(PID)
# find minimum area outside parent at differenr rotations
else:
angle=gss(0.0,90,1e-3)
# square fits if rotated
if f(angle)==0:
row[0]=ShapeMake(square,angle)
rows.updateRow(row)
initFidList.append(PID)
arcpy.SetProgressorPosition()
quer='"FID" IN '+str(initFidList)
quer='%s%s%s in %s'%('"',linkField,'"',str(tuple(initFidList)))
parent.definitionQuery=""
arcpy.SelectLayerByAttribute_management(child, "NEW_SELECTION", quer)
except:
message = "\n*** PYTHON ERRORS *** "; showPyMessage()
message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
message = "Python Error Info: " + str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()
Script works on shapefiles and assumes:
It child fits completely into original, script adds this square to selection. If not it attempts to fit it by rotating around centre points, and adds to selection if successful
Best Answer
If you have an Advanced license, you might try the Eliminate Tool, in Data Management, Generalization. Here is the tool description, "Eliminates polygons by merging them with neighboring polygons that have the largest area or the longest shared border."