Selecting parcels situated in different county than is written in attribute table using ArcPy [closed]

Question

I have two polygon feature classes: one with parcels, the other one with counties (40 counties)

I have to identify the parcels that in the attribute table belonging to county 'X', but spatially they are situated in another county.

Manually, it works with Data Reviewer extension in ArcMap: I create a BatchJob with a Geometry on Geometry Check where I write 40 expressions (because I have 40 counties!!!) like this: Feature Class1 = [PARCELS] where COUNTYID = 'X'; Feature Class2 = [COUNTIES] where COUNTYID = 'X'; SPATIAL RELATION: INTERSECT (checking the option 'Find features not in this relationship') --> and the output are the parcels that have COUNTYID = 'X', but spatially ARE NOT in the county 'X'.

Is there a way to do this automatically?

I thought about creating a list with the values of the COUNTYID field for both feature classes and then iterate through each value and select by location.

The problem is that I cannot figure it out automatically the option 'Find features not in this relationship' that DataReviewer has and SelectByLocation does not.

Answer 1

If creating a new shapefile using Spatial Join is not an option, below is some arcpy code (uses arcpy.da.SearchCursor so you'll need 10.1) that will identify which county each parcel is in. You will simply need to create a field called "ACTUAL_CTY" in your Parcels feature class that will hold the value of the county that the parcel is located in.

# The script below requires the Parcels feature class to have a field called ACTUAL_CTY
# You will need to add this field before running the script
# Adjust the location of your actual files in the first 2 lines of code

import arcpy

countiesFC = 'C:/YourDataFolder/Counties.shp'
parcelsFC = 'C:/YourDataFolder/Parcels.shp'

# Get a unique list of the counties
values = [row[0] for row in arcpy.da.SearchCursor(countiesFC, ("COUNTYID"))]
list_of_counties = set(values)

# Create a feature layer for the Parcels
arcpy.MakeFeatureLayer_management(parcelsFC, "lyr_Parcels")

for county in list_of_counties:
    print county
    # Create a feature layer that contains just the current county
    arcpy.MakeFeatureLayer_management(countiesFC, "lyr_" + county, "\"COUNTYID\" = '{0}'".format(county))

    # Select the parcels that have their center in the current county
    arcpy.SelectLayerByLocation_management("lyr_Parcels", "HAVE_THEIR_CENTER_IN", "lyr_" + county)

    # Write the current county to the ACTUAL_CTY field
    arcpy.CalculateField_management("lyr_Parcels", "ACTUAL_CTY", "'{0}'".format(county), "PYTHON_9.3")

# Select parcels that have an incorrect COUNTYID
arcpy.SelectLayerByAttribute_management("lyr_Parcels", "NEW_SELECTION", "\COUNTYID\" <> \"ACTUAL_CTY\"")

# Create a shapefile containing just the incorrect parcels
arcpy.CopyFeatures_management("lyr_Parcels", 'C:/YourDataFolder/Incorrect Parcels.shp')

Answer 2

Well, I don't know if you have to do this with the data in-place or if you can create a new output feature class with the results and then just link it back to your original dataset. If you can create a new feature class in the process of accomplishing this job, then the thought that comes to my mind would be to use the Spatial Join geoprocessing tool. What I would do in this case, first off, is check and make sure you have a unique identifier for each feature (other than the ObjectID or GlobalID fields); if you don't have such a field, no worries, just add a field to your main dataset with the field type of long integer and run the field calculate tool where the new field you created is set to the value of the objectid field. Then you would use that spatial join tool: join the parcels feature class and the counties feature class spatially, which if set up correctly should add the attributes, including name, of the counties feature class, to each parcel (leaving existing parcel attributes) and write all of that to a new output feature class. Then you could do a select by attribute on that new feature class to select all of the features where the parcel's county field is not equal to the county's name field. Then, if you really wanted to, you could extract those selected records to a new feature class or just a new table and then join that to your original parcel feature class so you could do a sort, select statement, etc... on your original feature class to get the specific records you need or do what you need with them. Personally, if you're going to try doing this and you're not sure about coding it in ArcPy directly, try using model builder (can be accessed from arccatalog or arcmap) and set up the operation that way, where you can drag-and-drop geoprocessing tools, feature classes, etc... . Then you can just select the option within model builder to export to python script.

Conceptually, it would probably look something like (this is just pseudo-code):

(Assuming unique ID field already exists, don't use OID or GlobalID as these can be different than their original values in an exported dataset or newly created dataset such as the output of spatial join)

Spatial Join (Target: Parcels, Join: Counties, Join Operation: One to Many - this
    allows parcels in multiple counties to be identified, field map: set up as 
    applicable for your data, others: defaults or as needed for your situation)
Select (Input: output from spatial join, Expression: ParcelCountyName <> CountyName)

This isn't necessarily perfect, depending on what you need, but figured I'd share in case it helps. Also, FYI: I'm using v10.1 to write this, so minor modifications may be needed if you're using a different version.