'Select Layer By Location' Not Working Properly

1615
3
10-11-2012 12:00 PM
AbbyFlory
New Contributor
(Note 1: I'm using Arc 10.1 and PyScripter)
(Note 2: Thank you, thank you, thank you in advanced for reading the monster explanation below!  I am extremely frustrated and would greatly appreciate any help!)

I'm not sure how to briefly sum-up my issue, but (very generally) I'm trying to assign values to a flag field in a point shapefile, where the flag field designates groups of points, and each group represents a collection of points that can be buffered to a given buffer size of interest without having their buffers overlap.  So basically, I'm sub-setting the point dataset into groups that, when buffered, have buffers that won't overlap. 

My logic follows a loop whereby:
1) All non-flagged features are selected,
2) Any features falling within "CombinedBuffers.shp" (see definition below) are removed from the selection,
3) First feature in remaining selection is chosen (i.e., basically random),
4) This feature is flagged (for ex., "1" - for group 1),
5) Then this feature is buffered, and the buffer is added (appended) to a "combined buffers" shapefile (each buffer added to this shapefile represents an area under which additional "group 1" points cannot fall; the shapefile is emptied when a new group starts), 
6) The code loops back to the top again - selecting all points not flagged, removing those under buffered areas, selecting the first feature, etc.
7) Once all possible points have been designated as group 1, CombinedBuffers.shp is emptied and the loop starts again on the subset of points not flagged (this time, assigning points to group 2, etc.)

My problem is introduced in #2 above: removing points from the current selection based on location.  (Note: Selection in #1 works fine.)  No error is thrown for selection #2, it simply does not remove any points when it should (proven by record count). 

What I've already tried/checked:
1) "Input" and "CombinedBuffers" datasets are shapefiles, but I've also tried them as feature classes in a personal geodatabase.
2) All data are in the same projected coordinate system; all data are stored in the same location.
3) Selection #1 followed by selection #2 works manually in ArcMap with tool interfaces AND the python window (using the same code!).
4) The appended buffer in CombinedBuffers.shp IS being recognized (feature counts are correct after each added buffer).
5) I have changed the type of selection, and it still will not select properly.
6) I have used similar selection sets in previous codes (and they worked), but I have not tested them using 10.1 yet.  Could this be a 10.1 issue?

My code is below (Note: sorry for the lengthiness - anyone willing to read through this is highly commended and your help incredibly appreciated!): 


#STEP 1: GENERAL SETTINGS-------------------------------------------------------

# Imports...
import time
print "Start time:",time.ctime()
StartSecs = time.clock()
import Functions
import arcpy
import sys
import traceback
import easygui
from arcpy import env
from arcpy.sa import *

# Variables...
BuffDist = 10
BuffUnit = "Kilometers"

Input = "C:/.../TestWells_TabAreaInput.shp"
Output = "C:/.../TestWells_TabAreaOutput.shp"
CmbdBuffs = "C:/.../CmbdBuffs.shp"
TempBuff = "C:/.../TempBuff.shp"

SpatialRef = arcpy.Describe(Input).spatialReference
Scratch_Path = "C:/.../Scratch"
CursorExists = "No"
RowExists = "No"

# Environments...
env.workspace = Scratch_Path
env.scratchWorkspace = Scratch_Path
env.overwriteOutput = "True"
env.extent = "C:/.../FC_Poly_IAState"
env.mask = "C:/.../FC_Poly_IAState"
env.maintainSpatialIndex = False

# Enable error handling...
try:

# STEP 2: CREATE POINTS FEATURE LAYER AND BUFFER FEATURE CLASS------------------

    # Make 'wells' feature layer and set flag field = 0...
    arcpy.MakeFeatureLayer_management(Input,"InputLyr")
    arcpy.CalculateField_management("InputLyr","BuffSelFlg",0)
    n = int(arcpy.GetCount_management("InputLyr").getOutput(0))
    print "Total (n) Records =", n

    # Create empty buffer feature class...
    Buffer = str(BuffDist*2) + " " + BuffUnit  # For these processes, buffers are defined as 2x buffer size of interest.
    if arcpy.Exists(CmbdBuffs):
        arcpy.Delete_management(CmbdBuffs)
    arcpy.Buffer_analysis("InputLyr",CmbdBuffs,Buffer)
    arcpy.DeleteFeatures_management(CmbdBuffs)

    BufferCnt = int(arcpy.GetCount_management(CmbdBuffs).getOutput(0))
    print "Buffer Count =", BufferCnt

    x = 1  # (assigned group value/flag)

# STEP 3: START LOOP------------------------------------------------------------
# Note: The use of "group" below is defined as a group of wells that can be buffered
# to the buffer size of interest without having their buffers overlap.

    while x <= n:
        # Select wells that have not been added to a group (i.e., where flag = 0)...
        arcpy.SelectLayerByAttribute_management("InputLyr","NEW_SELECTION",'"BuffSelFlg" = 0')
        RcdCnt = int(arcpy.GetCount_management("InputLyr").getOutput(0))
        print "Record Count where flag is zero =", RcdCnt

        if RcdCnt > 0:  # (i.e., if there are wells that have not been added to any group...)

            # Remove (from selection) wells within the current group's buffers...
            arcpy.SelectLayerByLocation_management("InputLyr","INTERSECT",CmbdBuffs,"","REMOVE_FROM_SELECTION")
            OutsideBuffRcdCnt = int(arcpy.GetCount_management("InputLyr").getOutput(0))
            print "Record Count where flag is zero and wells are outside buffers =", RcdCnt
            # NOTE: THIS IS THE SELECTION THAT DOES NOT WORK

            if OutsideBuffRcdCnt > 0:  # (i.e., if there are wells left that are qualified for assignment into this group...)
                UpdtCur = arcpy.UpdateCursor("InputLyr","",SpatialRef)
                row = UpdtCur.next()
                row.BuffSelFlg = x
                UpdtCur.updateRow(row)

                UnqID = row.UnqRecordI
                print "UnqID =",UnqID
                arcpy.SelectLayerByAttribute_management("InputLyr","NEW_SELECTION","\"UnqRecordI\" = '" + UnqID + "'")
                arcpy.Buffer_analysis("InputLyr",TempBuff,Buffer)
                arcpy.Append_management(TempBuff,CmbdBuffs,"NO_TEST")
                arcpy.AddSpatialIndex_management(CmbdBuffs)

                BufferCnt = int(arcpy.GetCount_management(CmbdBuffs).getOutput(0))
                print "Buffer Count =", BufferCnt

                arcpy.SelectLayerByAttribute_management("InputLyr","CLEAR_SELECTION")
                arcpy.Delete_management(TempBuff)
                del UpdtCur
                del row

            else: # (Start a new group)
                arcpy.SelectLayerByAttribute_management("InputLyr","CLEAR_SELECTION")
                arcpy.DeleteFeatures_management(CmbdBuffs)
                x = x + 1

        else: # (All records have been flagged into groups - job done!)
            break

# STEP 4: COPY FEATURES AND DELETE FEATURE LAYER----------------------------
    arcpy.CopyFeatures_management("InputLyr",Output)
    arcpy.Delete_management("InputLyr")


# ERROR HANDLING----------------------------------------------------------------

except arcpy.ExecuteError:
    # Get the tool error messages:
    msgs = arcpy.GetMessages(2)

    # Return tool error messages for use with a script tool:
    arcpy.AddError(msgs)

    # Print tool error messages for use in Python/Python Window:
    print "Geoprocessing tool error:"
    print msgs


except:
    # Get the traceback object:
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0]

    # Concatenate info together concerning the error into a message string:
    pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1])
    msgs = "ArcPy ERRORS:\n" + arcpy.GetMessages(2) + "\n"

    # Return python error messages for use in script tool or Python Window:
    arcpy.AddError(pymsg)
    arcpy.AddError(msgs)

    # Print Python error messages for use in Python/Python Window:
    print arcpy.GetMessages()
    print pymsg


finally:
    if CursorExists == "Yes":
        del InputCur
    if RowExists == "Yes":
        del row
    if arcpy.Exists("InputLyr"):
        arcpy.Delete_management("InputLyr")

    print "End time:",time.ctime()
    TimeElapsed, Step = (time.clock()-StartSecs), "END"
    Functions.PrintTime(TimeElapsed,Step)
Tags (2)
0 Kudos
3 Replies
MathewCoyle
Frequent Contributor
So your problem is nothing is ever removed from your selected features? Can you post your data you are running this on? Or at least the printout of the results from one of the runs? I'd also try using meters instead of kilometers.

Overall this seems like a very convoluted work flow, but I don't understand enough of your problem to comment on a way to improve it.
0 Kudos
AbbyFlory
New Contributor
Thanks Matthew,

Yes, my problem is that nothing is ever removed from my selected features.  I have tried changing the selection type (for example, to "WITHIN" the buffer, such that I'm selecting from a selection instead of removing from a selection), and it still does not work.  Below is a print out of the first few loops running on a subset of my data.  If this doesn't help I can post the data.  Meters vs. kilometers should make no difference (10 Kilometers vs. 10000 Meters).  The buffers are being created successfully (with proper distance value - proven when I open the shapefile in ArcMap), so the units do not seem to be an issue.  I agree that this might be a convoluted work flow, but I'm not sure how else to approach it. 

My real issue is this: I need to retrieve the area (in square units) of each land-cover type/category within each feature of a polygon feature class.  In my case, polygon features = buffers.  The "Tabulate Area" tool (Spatial Analyst > Zonal) does this, however there is a problem inherent in the tool: if any polygons overlap, one will get TabArea results for the full polygon while the other(s) will only get results for part of the polygon (the part that does not overlap).  I need areas for each full polygon. 

My original approach (and the approach of my predecessor) was to loop through each feature of interest (i.e. each polygon/buffer) one at a time, and run the TabArea tool individually on each feature.  This is quite taxing on implementation time - especially when you're dealing with 35,000+ features.  Originally written in VBA, this process took 3-4 seconds per feature.  I've gotten Python to mimic the VBA code and processing time (since ESRI is phasing-out VBA), but at 35,000 records, it will take 29-39 hours to complete a single run - and I need to run this process on multiple buffer sizes and multiple land-cover datasets (approx. 35 times).  While this code process does work, it's too inefficient for my current needs.  My predecessor and I decided that if we could place the features into non-overlapping groups, the tool could run on each group instead of each feature (drastically reducing implementation time).  This is the code I included with my original post.  Logically, it works (albeit, convoluted), and it can be used with any buffer size of interest.  However, as you're aware, I'm having unusual issues with the second selection set.  Like I originally posted, this process works in ArcMap - both with the selection tool interfaces and the Python window.

Let me know your thoughts, and thanks for your help!



BEGINNING STEP 1: GENERAL SETTINGS...
Start time: Fri Oct 12 09:42:10 2012
Step 1 COMPLETE! Time elapsed: 14.45 seconds

BEGINNING STEP 2:
Total (n) Records = 80
Buffer Count = 0
Record Count where flag is zero = 80
Record Count where flag is zero and wells are outside buffers = 80
UnqID = GTC_NO384734
Buffer Count = 1

Record Count where flag is zero = 79
Record Count where flag is zero and wells are outside buffers = 79
UnqID = GTC_NO385488
Buffer Count = 2

Record Count where flag is zero = 78
Record Count where flag is zero and wells are outside buffers = 78
UnqID = GTC_NO323237
Buffer Count = 3

Record Count where flag is zero = 77
Record Count where flag is zero and wells are outside buffers = 77
UnqID = USGS7073
Buffer Count = 4

Record Count where flag is zero = 76
Record Count where flag is zero and wells are outside buffers = 76
UnqID = GTC_NO357868
Buffer Count = 5

Record Count where flag is zero = 75
Record Count where flag is zero and wells are outside buffers = 75
UnqID = PWTS3399
Buffer Count = 6

Record Count where flag is zero = 74
Record Count where flag is zero and wells are outside buffers = 74
UnqID = PWTS753
Buffer Count = 7
0 Kudos
MathewCoyle
Frequent Contributor
You seem to be printing the same variable twice. Is this intentional or did you mean to print OutsideBuffRcdCnt the second time?

    while x <= n:
        # Select wells that have not been added to a group (i.e., where flag = 0)...
        arcpy.SelectLayerByAttribute_management("InputLyr","NEW_SELECTION",'"BuffSelFlg" = 0')
        RcdCnt = int(arcpy.GetCount_management("InputLyr").getOutput(0))
        print "Record Count where flag is zero =", RcdCnt

        if RcdCnt > 0:  # (i.e., if there are wells that have not been added to any group...)

            # Remove (from selection) wells within the current group's buffers...
            arcpy.SelectLayerByLocation_management("InputLyr","INTERSECT",CmbdBuffs,"","REMOVE_FROM_SELECTION")
            OutsideBuffRcdCnt = int(arcpy.GetCount_management("InputLyr").getOutput(0))
            print "Record Count where flag is zero and wells are outside buffers =", RcdCnt
0 Kudos