Unable to identify point features intersecting polygons using ArcPy in ArcGIS Pro

Description

In ArcGIS Pro, when retrieving point features intersecting polygon features using ArcPy, both intersecting and non-intersecting point features are returned. In this example, the script returns all point features displayed on the map instead of only the intersecting point features.

Cause

The CopyFeatures function references the initial point features instead of the intersecting point features identified through the spatial selection created by the SelectLayerByLocation function. As a result, the CopyFeatures function returns the initial point features instead of only the intersecting point features.

Solution or Workaround

Modify CopyFeatures to use the spatial selection from SelectLayerByLocation

1. Open the project in ArcGIS Pro.
2. On the View tab, in the Windows group, click Python Window.
3. In the Python window, copy and paste the code below to identify point features intersecting polygon features.
   1. Import the ArcPy library.

import arcpy

2. Define the main function and set the geoprocessing environment.

def main():
    arcpy.env.workspace = arcpy.env.scratchGDB
    arcpy.env.overwriteOutput = True

3. Define the paths of the input point and polygon feature classes. Replace <point_path> with the path of the point feature class and <polygon_path> with the path of the polygon feature class.

poly = r"<polygon_path>"
point = r"<point_path>"

4. Create feature layers from the input feature classes. Replace <point_layer> with the point feature layer name and <polygon_layer> with the polygon feature layer name.

arcpy.MakeFeatureLayer_management(point, "<point_layer>")
arcpy.MakeFeatureLayer_management(poly, "<polygon_layer>")

5. Retrieve the total number of point features.

total_points = int(arcpy.management.GetCount(point)[0])

6. Select the non-intersecting point features and calculate the intersecting point features.

total_points = int(arcpy.management.GetCount(point)[0])
  
    arcpy.SelectLayerByLocation_management(
        "point_lyr",
        "WITHIN",
        "poly_lyr",
        "",
        "NEW_SELECTION",
        "INVERT"
    )
    non_intersecting_points = int(
        arcpy.management.GetCount("point_lyr")[0]
    )
    intersecting_points = total_points - non_intersecting_points

    print(f"Total point features: {total_points}")
    print(f"Total intersecting point features: {intersecting_points}")

7. Export the non-intersecting point features.

point_select = "point_select"

    arcpy.CopyFeatures_management(
        "point_lyr",
        point_select
    )
if __name__ == "__main__":
    main()

The code below shows the sample script.

import arcpy

def main():

    arcpy.env.workspace = arcpy.env.scratchGDB
    arcpy.env.overwriteOutput = True

    poly = r"test.gdb\test_poly"
    point = r"test.gdb\test_point"

    arcpy.MakeFeatureLayer_management(point, "point_lyr")
    arcpy.MakeFeatureLayer_management(poly, "poly_lyr")

    total_points = int(arcpy.management.GetCount(point)[0])

    arcpy.SelectLayerByLocation_management(
        "point_lyr",
        "WITHIN",
        "poly_lyr",
        "",
        "NEW_SELECTION",
        "INVERT"
    )

    non_intersecting_points = int(
        arcpy.management.GetCount("point_lyr")[0]
    )

    intersecting_points = total_points - non_intersecting_points

    print(f"Total point features: {total_points}")
    print(f"Total intersecting point features: {intersecting_points}")

    point_select = "point_select"

    arcpy.CopyFeatures_management(
        "point_lyr",
        point_select
    )
if __name__ == "__main__":
    main()

4. Press Enter twice to run the code.

The output message displays the number of point features intersecting polygons.