Finding points outside state boundary that are supposed to be within specific State using ArcPy with ArcMap

Question

I have a nationwide featureclass called bakeries. In the bakeries FC, there is a field called "State" and the field is populated with state abbreviations(AL,NJ,CT,etc). I have another feature class called United States, and that FC also has a field called "State" but it is populated with full state name.

The bakeries dataset has some poorly geocodeded points that say that they are in a state(Alabama for example) but the location of the points falls outside the state boundary, sometimes very far away.

Here is a snippet of what I have so far, you'll see that the code works but returns me everything that falls outside of the state versus the state specific data.

My goal is to identify points that say they are within a state, but in actuality, fall outside the boundary.

for row in state_rows:

    arcpy.AddMessage(str(new_name + "=" + "'" + row.getValue("STATE") + "'"))        
    sqlExp = new_name + "=" + "'" + row.getValue("STATE") + "'" #making SQL Statement      

    States2FL = arcpy.MakeFeatureLayer_management(States, outworkspace + "\\" + row.getValue("STATE")+ "_FL_FL",sqlExp) #Making Feature Layer of States.shp        
    FC2FL = arcpy.MakeFeatureLayer_management(fc.CatalogPath,outworkspace + "\\" + row.getValue("STATE") + "_FC_FL") #Making Feature Layer of FC

    arcpy.SelectLayerByLocation_management(FC2FL,"",States2FL,"","NEW_SELECTION")#selecting all data inside respective selected state                                               
    arcpy.SelectLayerByLocation_management(FC2FL,"","","","SWITCH_SELECTION") #switch selection        

    record_countAF = int(arcpy.GetCount_management(FC2FL).getOutput(0)) #Get After Record count        
    percent = (float(record_countAF)/float(record_countb4)) * 100.00 #calc percent here

Answer 1

Using a spatial join is the easiest, but it creates an additional layer...

import arcpy

# Set overwrite option
arcpy.env.overwriteOutput = True

# USING SPATIAL JOIN
inputFC = "C:/YourFolder/Bakeries.shp"
joinFC = "C:/YourFolder/States.shp"
outputFC = "C:/YourFolder/BakeriesWithStates.shp"
arcpy.SpatialJoin_analysis(inputFC, joinFC, outputFC)

However, sometimes it's not desirable to create an additional layer. Maybe you just want to create another field in your original layer (bakeries.shp in this case) and store the "located state". Here is how that is done...

import arcpy

# Set overwrite option
arcpy.env.overwriteOutput = True

# KEEPING ORIGINAL LAYERS (NOT CREATING ADDITIONAL "JOINED" LAYER)
# Create FeatureLayers
arcpy.MakeFeatureLayer_management("C:/YourFolder/Bakeries.shp", "lyr_Bakeries")
arcpy.MakeFeatureLayer_management("C:/YourFolder/States.shp", "lyr_States")

# Add an "LOCATED_ST" field
if "LOCATED_ST" not in [f.name for f in arcpy.ListFields("lyr_Bakeries")]:
    arcpy.AddField_management("lyr_Bakeries", "LOCATED_ST", "TEXT", "", "", "20")

# Create a search cursor for the states
rows = arcpy.SearchCursor("lyr_States")
for row in rows:
    # What you'll do is select each state one at a time, and then select all the bakeries in that state and calculate the LOCATED_ST field
    # NOTE: If you are using not using shapefiles, then you'll have to change the FID in the line below to OBJECTID (or similar)
    arcpy.SelectLayerByAttribute_management("lyr_States", "NEW_SELECTION", "\"FID\" = " + str(row.getValue("FID")))
    arcpy.SelectLayerByLocation_management("lyr_Bakeries", "INTERSECT", "lyr_States", "", "NEW_SELECTION")
    arcpy.CalculateField_management("lyr_Bakeries", "LOCATED_ST", "'{0}'".format(str(row.getValue("STATEFIELD"))), "PYTHON_9.3", "")
    print "Finished processing " + str(row.getValue("STATEFIELD"))

Obviously you'll have to substitute your layer paths as needed... and if you are using a file geodatabase instead of a shapefile, then you'll need to change the FID to OBJECTID (as commented...) When you're done, any empty values in the LOCATED_ST field means that it was outside of all of the states.

You'll note that I used forward slashes in my path which looks "wierd" for windows users, but I like it better that way because you don't need to escape the backslashes... (the forward slash works that same way that you would expect a backslash to work...)

Answer 2

If you want to use python or just tools straight out of arc toolboxes, I would recommend using a spatial join. This will create a field containing the state which the points fall within. This can then be compared toursor the state field in the point fc. You can create a comparison in python as well, using a cursor.

Answer 3

You can use spatial queries, or a combination of spatial and attribute queries. For example, attribute query the points by their state attribute, this will be a New Selection. Next, using the select by location, select points by the state boundary (completely within or intersect). For this second query use Remove From Selection instead of New Selection. The difference are all the points outside of the boundary.