For loop select by attribute and location...can't get it to work

418
3
01-26-2023 08:22 AM
by Anonymous User
Not applicable

I have a counties feature class and dozens of water well feature classes in a project geodatabase. I want to find all the water wells that are assigned to a county but have coordinates that place them outside of that county. My chunk of code just creates a new feature class that's a copy of the original water well feature class. I don't see what's not why this isn't working. I've played around with the selection and inversion clauses to try to get it right. If I invert both I can find wells that are not located within ANY Michigan county, but that's not particularly helpful. Have any ideas why this isn't working right?

```

print("Finding water well outliers")
arcpy.env.workspace = r"C:\Users\BarnbyJ\Documents\ArcGIS\Projects\Water Well Outliers Script\Water Well Outliers Script.gdb"
fcs = arcpy.ListFeatureClasses('*Wells*')
if fcs is not None:
for fc in arcpy.ListFeatureClasses("*_WaterWells"):
countyName = fc.split(".")[0] #assuming that the feature classes all have the same pattern
sel_counties = arcpy.management.SelectLayerByAttribute("Counties", "NEW_SELECTION", "NAME = '{}'", "NON_INVERT").getOutput(0)
sel_wells = arcpy.management.SelectLayerByLocation(fc, "INTERSECT", sel_counties, None, "ADD_TO_SELECTION", "INVERT").getOutput(0)
arcpy.management.CopyFeatures(sel_wells, "{}_WaterWells_Outliers".format(countyName), '', None, None, None)
arcpy.management.SelectLayerByAttribute("Counties", "CLEAR_SELECTION").getOutput(0)
print(sel_counties)
print("Outliers identified")

 

0 Kudos
3 Replies
DavidSolari
Occasional Contributor III

Did you forget the format call on the SQL string in line 7 in the script? If that's just a copy/paste typo, try wrapping the Counties polygons and each Well FC in an arcpy.management.MakeFeatureLayer call and work with the output layer, there's a good chance your selection is silently failing on a non-layer input. Another method is to merge all your wells into one big FC, spatially join it to the Counties, then select every record where the county name doesn't match the original FC name (the "add_source" parameter on the Merge tool makes this easy if your version of Pro supports it). This avoids fiddling with dynamic selections and should run much faster due to batching all the geometry calculations.

JoshuaBixby
MVP Esteemed Contributor

I don't understand these steps:

fcs = arcpy.ListFeatureClasses('*Wells*')
if fcs is not None:
    for fc in arcpy.ListFeatureClasses("*_WaterWells"):

Anything that Line 3 lists, Line 1 will also list, but the opposite isn't true.   Since the posted code does nothing to a feature class if it is listed in Line 1 but not Line 3, why bother checking?

0 Kudos
by Anonymous User
Not applicable

'ADD_TO_SELECTION' works on the layer that is getting selected, if it already has a selection. Since you are creating the selection on the counties layer first (by attribute), this wells layer is getting a NEW_SELECTION based on the intersect instead of the ADD_TO. Hope that makes sense.

I think you are missing a step though, as @DavidSolari and @JoshuaBixby brought up-making a Featurelayer using the attribute selections on both county and wells first, can set you up for the select by location and invert. *Wells* is the equivalent of sql's %Wells%, so it will include featureclasses with _WaterWells in the name.

fcCounties = r'CountiesLayer'

for fc in arcpy.ListFeatureClasses("*_WaterWells"):
    countyName = fc.split(".")[0] #assuming that the feature classes all have the same pattern
    # create selection on the county layer by sql
    sel_county = arcpy.management.MakeFeatureLayer(fcCounties, "county", f"NAME = '{countyName}'")
    # select the wells that are 'assigned' to the county by sql(?)
    sel_wells_by_attribute = arcpy.management.MakeFeatureLayer(fc, 'selwells', f"COUNTY = '{countyName}'")
    # Intersect the filtered county layer and filtered well layer, removing the intersecting (wells that are within the county) leaving the outliers
    outliers = arcpy.management.SelectLayerByLocation(sel_wells_by_attribute, "INTERSECT", sel_county, '', "NEW_SELECTION", "INVERT")
    # Copy outliers
    arcpy.management.CopyFeatures(outliers, f"{countyName}_WaterWells_Outliers")

 

0 Kudos