3

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.

  • I used the answer from Richard Fairhurst and it does a great job of creating a new shapefile with the information from the polygon. However, is there any way that the initial point file(pointLayer) can inherit the new field data instead of creating a whole new shapefile? – sbell423 Mar 9 '16 at 15:09
  • Yes - I wrote a script to do approximately this: pastebin.com/9LyBHX85 It requires having a blank field in the initial point file, but that's straightforward to make. – Andy Bradford Mar 9 '16 at 15:17
  • #Layer to be calculated InLayer = arcpy.GetParameterAsText(0) #InField: Layer which will receive final data InField = arcpy.GetParameterAsText(1) #SourceLayer: Layer which contributes data. SourceLayer = arcpy.GetParameterAsText(2) I'm a little confused, I only have one polygon and one point feature class, why are there 3 inputs? Should inField be a field within my points feature class? – sbell423 Mar 9 '16 at 15:29
  • ...thanks for catching that, sorry for the confusion! Yes, InField is the field in InLayer to be populated with the data from SourceLayer. – Andy Bradford Mar 9 '16 at 15:42
  • If you have a new question, please ask it by clicking the Ask Question button. Include a link to this question if it helps provide context. - From Review – jbchurchill Mar 9 '16 at 16:19
2

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
  • +1 For some great info on the cursor/dictionary techniques but personally I think the asker only needs to do a point in polygon overlay operation using something like Intersect or Spatial Join. I've added the results of a performance test that may interest you at gis.stackexchange.com/questions/115096/… – PolyGeo Sep 25 '14 at 3:49
  • Within SDE I am finding that I prefer not to replace data with entirely new versions. This is especially true of complex feature classes, like topologies, feature linked annotation, composite networks, network analyst datasets, etc. These are increasingly attractive and do not play well with complete replacements. So this code addresses that. In addition I just adapted it to replace about 10 tools that take about 10 minutes to process to create several levels of summary values and now processes in about 15 seconds. – Richard Fairhurst Sep 25 '14 at 4:45
  • The above script keeps returning the error: ExecuteError: ERROR 999999: Error executing function. The coordinates or measures are out of bounds. Failed to execute (SpatialJoin). I have checked that all the layers are in the same CRS and I do not get this error when I run SpatialJoin in model builder or as a tool. Could this error be caused by not defining my field_mapping? – wschwenger Oct 1 '14 at 17:00
  • The field list does not come into play with that error. Try setting up the Extent in the Geoprocessing menu Environment or adding an environment setting in the script. I would also export your modelbuilder model to a python script and see if it runs. It will add every parameter to the output and you could try pasting them into this script to see if including a particular parameter makes a difference. – Richard Fairhurst Oct 1 '14 at 17:42
  • Thanks Richard, my polygon layer was imported from CAD so it had a Z value. Once the z was removed it works. – wschwenger Oct 1 '14 at 17:55
1

I would do this using the much faster (than ArcPy) geoprocessing tool Intersect (Analysis) which:

Computes a geometric intersection of the input features. Features or portions of features which overlap in all layers and/or feature classes will be written to the output feature class.

Use your point and polygon features as the in_features and you should get an output point feature class which has the polygon attributes attached.

If you need to iterate through different cities then just do the first one from the tool dialog and then use Copy As Python Snippet in the Geoprocessing | Results window to copy/paste into your Python script that will then just need a for loop wrapped around it to iterate through a list of your cities.

  • Or Spatial Join. In this case I would use Near to populate the points with the polygon they're nearest to then do an attribute join followed by field calculation - much faster! – Michael Stimson Sep 25 '14 at 0:50
  • @MichaelMiles-Stimson I've added the results of a performance test that may interest you at gis.stackexchange.com/questions/115096/… – PolyGeo Sep 25 '14 at 3:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.