3

I have been trying for some time now to create a script to extract the polygons within an extent, but am unable to get the same result as when I use the user interface manual method.

Maybe someone can see the fault of my script, I have stripped it down to the minimum for simplicity.

import arcpy

layer = arcpy.GetParameterAsText(0)
top = arcpy.GetParameterAsText(1)
left = arcpy.GetParameterAsText(2)
bottom = arcpy.GetParameterAsText(3)
right = arcpy.GetParameterAsText(4)
outputFeature = arcpy.GetParameterAsText(5)

# Determine if Layer is a Layer
desc = arcpy.Describe(layer)
if hasattr(desc, "layer"):
     lyr = layer
     arcpy.AddMessage('Parameter was layer: {0}'.format(lyr))
else:
     lyr = arcpy.MakeFeatureLayer_management(layer,"in_memory\\source_layer")
     arcpy.AddMessage('Parameter made to layer: {0}'.format(lyr))
     
sr = arcpy.CreateSpatialReference_management("",lyr,"","","","")

lowerLeft = arcpy.Point(left,bottom)
upperLeft = arcpy.Point(left,top)
upperRight = arcpy.Point(right,top)
lowerRight = arcpy.Point(right,bottom)

a = arcpy.Array([lowerLeft,lowerRight,upperRight,upperLeft,lowerLeft])
thepoly = arcpy.Polygon(a,sr)

arcpy.AddMessage('Selecting form lyr by polygon: {0}'.format(thepoly))
selection = arcpy.SelectLayerByLocation_management(lyr, "WITHIN", thepoly, "", "NEW_SELECTION")

matchcount = int(arcpy.GetCount_management(selection).getOutput(0)) 
if matchcount == 0:
    arcpy.AddMessage('no features matched spatial and attribute criteria')
else:
    if matchcount < 1000:
        arcpy.CopyFeatures_management(selection, outputFeature)
        arcpy.AddMessage('{0} objects that matched criteria written to {1}'.format(matchcount, outputFeature))

Follow on after some responses

Whatever I try (and I've tried lots of different combinations) I don't get the desired results. Unless I use the GUI interface to SelectByLocation

Screen shot of GUI, just to clarify

The selected features produced by the GUI interface is exactly what I would expect. I've change the code a bit to output thepoly as feature class (BoundingAOI) so I could check that input parameters were correctly used.

import arcpy

try:

    layer = arcpy.GetParameterAsText(0)
    top = arcpy.GetParameter(1)
    left = arcpy.GetParameter(2)
    bottom = arcpy.GetParameter(3)
    right = arcpy.GetParameter(4)
    outputFeature = arcpy.GetParameterAsText(5)
    
    # Determine if Layer is a Layer
    desc = arcpy.Describe(layer)
    if hasattr(desc, "layer"):
         lyr = layer
         arcpy.AddMessage('Parameter was layer: {0}'.format(lyr))
    else:
         lyr = arcpy.MakeFeatureLayer_management(layer,"in_memory\\source_layer")
         arcpy.AddMessage('Parameter made to layer: {0}'.format(lyr))
         
    sr = arcpy.CreateSpatialReference_management("",lyr,"","","","")
    
    lowerLeft = arcpy.Point(left,bottom)
    upperLeft = arcpy.Point(left,top)
    upperRight = arcpy.Point(right,top)
    lowerRight = arcpy.Point(right,bottom)
    
    a = arcpy.Array([lowerLeft,lowerRight,upperRight,upperLeft,lowerLeft])
    thepoly = arcpy.Polygon(a)
    
    outFeatureClass = arcpy.CreateFeatureclass_management(arcpy.env.workspace, "BoundingAOI", "Polygon", "#","#","#", sr)
    arcpy.AddMessage('Created {0}'.format(outFeatureClass))
    cur = arcpy.da.InsertCursor(outFeatureClass, ["SHAPE@"])
    cur.insertRow([thepoly])
    
    arcpy.AddMessage('Selecting form lyr by polygon: {0}'.format(thepoly))
    selection = arcpy.SelectLayerByLocation_management(lyr, "COMPLETELY_WITHIN", outFeatureClass, "", "NEW_SELECTION")
    
    matchcount = int(arcpy.GetCount_management(selection).getOutput(0)) 
    if matchcount == 0:
        arcpy.AddMessage('no features matched spatial and attribute criteria')
    else:
        if matchcount < 1000:
            arcpy.CopyFeatures_management(selection, outputFeature)
            arcpy.AddMessage('{0} objects that matched criteria written to {1}'.format(matchcount, outputFeature))
            
except Exception, e:
    import traceback
    map(arcpy.AddError, traceback.format_exc().split("\n"))
    arcpy.AddError(str(e))

I'm new to ArcGIS and Python (background in Fortran, C, perl and C#).

Is there any way of step-by-step debugging ArcPy scripts?

  • Have applied service pack 1 to ArcGIS 10.1, so am now running 10.1sp1
  • Got PyScripter running and able to debug from ArcMap (did need to copy the contents of the rpyc.zip file from the PyScripter install to a rpyc directory under ArcGIS installation of Python C:\Python27\ArcGIS10.1\lib\site-packages\rpyc on my system), this maybe of use to someone
  • Cast the double input parameters as double as suggested
  • No errors from the script
  • Below is a screen shot of the script GUI as requested

result pane and script GUI

My input parameters are (on the comand line)

"D:\Current\PaddockGRASP VS2012\PaddockGRASP Rev 1\PaddockGRASP.gdb\QLD_CADASTRE_DCDB_20" -18.933367 145.278049 -19.735903 146.119615 "D:\Current\PaddockGRASP VS2012\PaddockGRASP Rev 1\PaddockGRASP.gdb\OutputFeatures"

It appears after all this that the input dataset had some sort of corruption (ArcSDE can't find anything wrong, not sure why or what is actually wrong), found another older version of the dataset and the problem went away.

0

2 Answers 2

2

I think @artwork21's suggestion was to try casting your left/top/bottom/right parameters to numbers, e.g.:

top = float(arcpy.GetParameterAsText(1))
left = float(arcpy.GetParameterAsText(2))
bottom = float(arcpy.GetParameterAsText(3))
right = float(arcpy.GetParameterAsText(4))

However, this is not necessary according to my testing (see below).

I would suggest stepping through your code in a Python debugger, such as PyScripter's. You can supply your arguments as command line parameters or hardcode them for testing purposes:

Command Line Parameters menu item Command Line Parameters dialog

There is more on debugging tools at Debugging ArcPy scripts.

I tested your original script as a script tool, and it works fine for me. Perhaps check the set up of your script tool parameters. This works for me (all but the last are Inputs, the last is an Output parameter). I am on 10.1 SP1.

Script tool parameters dialog

Also tested hardcoded parameter values as a standalone script in PyScripter:

layer = r"C:\GISData\12345.shp"
top = "394605.709522"
left = "13411007.111189"
bottom = "392946.646919"
right = "13413137.265642"
outputFeature = r"C:\GISData\12345_copyFeaturesByRectangle.shp"
0
1

The problem is that these variables (left, top, bottom, right):

lowerLeft = arcpy.Point(left,bottom)
upperLeft = arcpy.Point(left,top)
upperRight = arcpy.Point(right,top)
lowerRight = arcpy.Point(right,bottom)

need to be numeric not text. Surround your code with a try/except error handeling will help you figure out potential coding errors e.g.

try:

  # your code in here

except Exception, e:
  import traceback
  map(arcpy.AddError, traceback.format_exc().split("\n"))
  arcpy.AddError(str(e))
0

Not the answer you're looking for? Browse other questions tagged or ask your own question.