Script somehow deleting all points in a shapefile (pour point automation)

[RyanUtz]

Hello,

I am very new to Python scripting- and used to program in Arc/INFO (how I miss it so these days). All I am trying to do is write a script that will 1) import a list of integers from a CSV file, 2) use this list to select points in a shapefile with matching values, and 3) run the hydro tool "Snap pour point" & save the pour point. I thought that I had made great progress yesterday and had the code completed (thanks to Curtis). However, the code below produces nonsensical, non-existent pour points and, for some reason, the code deleted all of the entries in the point shapefile I pointed it to!!! The shapefile still exists, but is completely empty. And the SnapPourPoint command seems reasonably easy enough; thus I have no idea why it results in completely empty entities...

Please help a newcomer!! This is driving me crazy and I am really not becoming a Python fan... -CSV file: "E:/Weathering_GIS/Sheds.csv", point shapefile [that gets emptied for some reason]: "E:/Weathering_GIS/Working.gdb/USGS_Points", flow accumulation grid: "FAC_NORTHEAST"

Thanks,
Ryan

Code:

import csv, arcpy
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")


env.workspace = 'E:/Weathering_GIS/Working.gdb'

try:
   f = open('E:/Weathering_GIS/Sheds.csv')
   sheds = csv.reader(f)
   for row in sheds:
      arcpy.MakeFeatureLayer_management('E:/Weathering_GIS/Working.gdb/USGS_Points', 'Shedpoints')
      where = "GageCode = %s" % int(row[0])
      arcpy.SelectLayerByAttribute_management("Shedpoints","NEW_SELECTION",where)  
      snapped = SnapPourPoint("Shedpoints","FAC_NORTHEAST",0.0001)
      out = "X%s" % str(row[0])
      snapped.save("E:/Weathering_GIS/Working.gdb/"+out)

except:
    del sheds
    f.close()
    raise

[ChrisSnyder]

I see you are using a distance of "0.0001" as your pour point snapping tolerance... What units are these in? Is the projection of "USGS_Points" defined?

I suspect you may be using data in a geographic coordinate system... If so, I would recomend projecting your input data... Something like Albers if your data covers the entire U.S. - otherwise if your data is more local, use State Plane or UTM.

Some of the Grid functions end up a bit screwy if you use geographic coordinate systesm.

[StacyRendall1]

Hi Ryan,

I don't know anything about Spatial Analyst, but I did notice that you remake your feature layer every row of the CSV - which isn't necessary as you have defined your Select by Attributes to make a new selection. I have taken your code and tidied it out a bit; you didn't use the workspace so I deleted those lines, got rid of the try...except stuff and added an explicit delete if your output already exists, rather than using Overwrite Output (which, for some reason, may be allowing it to mess with your input file). Here is the code, let me know if you have any more success with it, or if you have any questions about the changes I made:

import csv, arcpy
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")

# Make the feature layer (once)
arcpy.MakeFeatureLayer_management('E:/Weathering_GIS/Working.gdb/USGS_Points', 'Shedpoints')

# open CSV and iterate over rows
f = open('E:/Weathering_GIS/Sheds.csv')
sheds = csv.reader(f)
for row in sheds:

  # generate where clause, then make selection
  where = "GageCode = %s" % int(row[0])
  arcpy.SelectLayerByAttribute_management("Shedpoints","NEW_SELECTION",where)

  # snap pour points and generate output path
  snapped = SnapPourPoint("Shedpoints","FAC_NORTHEAST",0.0001)
  out = "X%s" % str(row[0])

  # if output already exists, delete it, then save
  if arcpy.Exists("E:/Weathering_GIS/Working.gdb/"+out):
    arcpy.Delete_management("E:/Weathering_GIS/Working.gdb/"+out)
  snapped.save("E:/Weathering_GIS/Working.gdb/"+out)

# close open files and delete variables
f.close()
del sheds, snapped, out, where