I have two feature classes, a polygon of sales areas for a city, and point layer of addresses. I want to select the addresses based on the sales area they are in and update their attribute table to show the sales area. I have searched python+point+in+polygon for other posts and I know how to select the polygons but I can't figure out how to populate the address points.
I found the following script that works but takes almost two hours to run.
import arcpy
import os
# Allow overwrite
arcpy.env.overwriteOutput = True
# Script user input parameters
polygonLayer = arcpy.GetParameter(0)
pointLayer = arcpy.GetParameter(1)
workspace = arcpy.GetParameter(2)
copyField = = arcpy.GetParameterAsText (3)
# Add new field to target feature class
fieldList = arcpy.ListFields(pointLayer)
fieldNameList = []
for field in fieldList:
fieldNameList.append(field.name)
fieldType = arcpy.ListFields(polygonLayer,copyField)[0].type
if not copyField in fieldNameList:
arcpy.AddField_management(pointLayer, copyField, fieldType)
arcpy.MakeFeatureLayer_management(polygonLayer, "PolygonLayer")
pntCursor = arcpy.UpdateCursor(pointLayer)
for pnt in pntCursor:
rowID = pnt.OBJECTID
arcpy.MakeFeatureLayer_management(pointLayer, "PointLayer", '"OBJECTID" = ' + str(rowID))
# Polygon selection based on CONTAINS selection method
arcpy.SelectLayerByLocation_management("PolygonLayer", "CONTAINS", "PointLayer")
# Copy field values from polygon to point layer
polyCursor = arcpy.SearchCursor("PolygonLayer")
for poly in polyCursor:
pnt.setValue(copyField, poly.getValue(copyField))
pntCursor.updateRow(pnt)
arcpy.Delete_management("PointLayer")
del pntCursor, polyCursor, pnt, poly
I don't think I need to create a new field so I have simplified my script too:
# Import arcpy module
import arcpy
import os
# 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
# Select all Address in FSA
polyCursor = arcpy.SearchCursor("PolygonLayer")
for poly in polyCursor:
# Select based on CONTAINS
arcpy.SelectLayerByLocation_management(poly, "CONTAINS", "PointLayer")
# Add FSA to Address
pntCursor = arcpy.UpdateCursor("pointLayer")
for pnt in pointCursor:
pnt.setValue(fsaPnt, fsaPoly)
pntCursor.updateRow(pnt)
del pntCursor, polyCursor, pnt, poly
I know my update cursor is wrong but i'm stuck on how to update one feature class using information selected from the other.
I know I can do this manually but I need to do it for multiple cities so I would like to automate the process.
I am using Arc 10.2.2 to run my python script.
Answer
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
No comments:
Post a Comment