Friday, 13 July 2018

python - Best way to find points outside state boundary that are supposed to be within specific State


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



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...)


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...