Latest Contributions by KimOllivier

If all you want is the bearing differences of vertices within a single polyline then a cursor can step through the vertices and calculate the angles like this: # vertex_angle.py # find change in angle along a polyline at each vertex # assume : # in a projected coordinate system, not dd # don't know what to do with output, maybe add value to a point # Kim Ollivier # 29 October 2010 # import arcgisscripting import math import sys gp = arcgisscripting.create(9.3) try:     inlay = sys.argv[1] except :     inlay = "d:/work/test.gdb/wiggle" desc = gp.Describe(inlay) shapefield = desc.ShapeFieldName cur = gp.SearchCursor(inlay) row = cur.next() n = 0 m = 0 p = 0 while row:     feat = row.getValue(shapefield)     print "Feature",n     for partNum in range(feat.partCount) :         part = feat.getPart(partNum)         ptLast = None         bearingLast = None         print "Part",m         for ptNum in range(part.count):             pt = part.next()             if ptLast:                 bearing = math.atan2((pt.Y - ptLast.Y),(pt.X - ptLast.X))                 if bearingLast:                     delta = bearing - bearingLast                     print p,delta/math.pi*180.0                 bearingLast = bearing             ptLast = pt             p+=1         m+=1     n+=1     row = cur.next() del cur,row Otherwise getting the vertices required is a lot harder if you do not have arc-node topology available. You would have this in a coverage if you built for nodes. This provides a unique list of nodes and tags on each "polyline" for the node at each end in the attribute tables AAT and NAT. I could then run a valence script to count the number of segments at each node to assist with the angle calculations. A simple extension of the valence concept can tag nodes as sources, sinks or passthrough. I know that a network dataset has all these capabilites, but the relationships and status are not accessible for python scripting, so I build my own network attributes. To build arc-node topology that is accessible for scripting in ArcGIS there is a very good tool in ET GEOTools called RenodePolylines. To get implied nodes at a T junction, clean the layer to a temporary layer using ETGeotools CleanPolyline in the absence of workstation clean. I would make sure there were no multi-parts to complicate the calculations. It is easy to then select each node and do some angle calculations between each segment. Rather than run a buffer I would just look up the last or first point of a polyline and then step back one vertex. A useful python function math.atan2() gives the absolute bearing of two points. The next issue is to get direct access to the polyline verticies. Opening a cursor on a featureclass is not enough, so I open it once and load the parts of the featureclass into Python list and dictionary structures. After processing in double-quick time, then open an update cursor and write out the results. It all sounds like reinventing Network in Python, and it is. Maybe this could be done with C#.NET and ArcObjects? I challenge anyone with idle time to a competition! ... View more

I suggest that you have a look at the logging standard library module. This is easy to setup to record each use of the tool into a separate file for later review. ... View more

The easiest way is to use the existing [ARCINFO Licence] ArcTool CalculateEndTime which moves the next rows datetime to the previous row. Then you can use simple subtraction to get the difference. If you look up the python datatype for datetime you can see that you can show the interval in any units. If you have not got an ArcInfo licence, then a short python script would get around this meanness. Read the dates in in sort order with a key and put in a dictionary. Then write out the dates using the dictionary as a lookup. ... View more

Use the tool DataManagement>Features>MinimumBoundingGeometry. This will create polygon rectangles with attributes for the max length and orientation without coding. ... View more

Use the numpy module , already included in ArcGIS 10 import numpy a = numpy.array((xa,ya,za)) b = numpy.array((xb,yb,zb)) diagDist = numpy.linalg.norm(a-b) ... View more

Try turning of the ArcGIS Desktop Index Service. This has been a strange fix, but proven for this error message. (Use Google desktop indexing for searching if you need to search) ... View more

There is a soft interrupt. Click on the Pythonwin icon in the tray at the bottom of the screen and you will get an option "break into running code" that will do a soft interrupt, very useful if you have got into a loop by forgetting the row = cur.next(). ... View more

Try leaving out the underscore in the field name. ... View more

You must have a strange installation. The os module works for me without any hacking. My pythonpath is set to C:\arcgis\bin (which I am not all that comfortable with) but I don't seem to need to set it to c:\python25\lib as well to work as a script. The os module behaves at 10 for me as well. ... View more

Load the ElementTree module. This is very helpful to parse XML files of metadata. Here is an example where I read in an old metadata file, fiddle with the dates and rewrite it. # LoadMetadata.py # with altered dates for current month # using element tree # create original Metadata with same name as layer # run this to alter to name_et.xml # reload into filegeodatabase # Note there is no tool to unload metadata # 15 March 2010 import arcgisscripting,sys,os import elementtree.ElementTree as ET import sys,os,datetime print print def alter(xmlfile,edDate,publishDate,createDate) :     """     read xml file for featureclass or table     change dates to today,loading date and  extract date     empty processing logs     write out file with _et suffix     return file name     """     print xmlfile     tree = ET.parse(xmlfile)     ## print tree.getroot().tag, tree.getroot().text,tree.getroot().tail,tree.getroot().attrib     # Edition Date     elem = list(tree.iter("resEdDate"))[0]     # print elem.tag,elem.text     elem.text = edDate     ## print elem.text     # Reference Date 001 (Creation)     elem = list(tree.iter("refDate"))[0]     # print elem.tag,elem.text     elem.text = createDate     ## print elem.text     # note there may be two of these dates     # DateTypCd 001 and 002     # Reference Date 002 (Publication)     if len(list(tree.iter("refDate"))) > 1 :         elem = list(tree.iter("refDate"))[1]         # print elem.tag,elem.text         elem.text = publishDate     else :         print "skipping publication date",xmlfile     ## print elem.text     # clear out lineag if it exists     try :         lin = list(tree.iter("lineage"))[0]         # print lin.tag         lin.clear()         lin.text = "Cleared"     except :         gp.AddMessage("skipping clear lineage")     outfile = xmlfile.replace(".","_et.")     tree.write(outfile)     return outfile # ---------------------- main ---------------------- try :     publishDate = sys.argv[1]     createDate  = sys.argv[2]     if createDate == '#' :         createDate = publishDate     gp.AddMessage(publishDate+type(publishDate)) except :     today = datetime.datetime.now()     firstSat = today.replace(day=1) + datetime.timedelta(5 - datetime.datetime.now().replace(day=1).weekday())     publishDate = firstSat.strftime("%Y%m%d")     createDate  = publishDate # override publishDate = '20100804' createDate = '20100710' gp = arcgisscripting.create(9.3) os.chdir("e:/crs/metadata") edDate = str(datetime.datetime.now().date()).replace("-","") edDate = '20100914' gp.AddWarning(edDate+" edit date") gp.AddWarning(publishDate+" publish date") gp.AddWarning(createDate+" create date") ws = "e:/crs/corax.gdb" metasrc = "e:/crs/metadata" gp.Workspace = ws os.chdir(metasrc) print print ws print metasrc print lstFC = gp.ListFeatureClasses("*") for fc in lstFC :     # print fc     fcxml = metasrc+"/"+fc+".xml"     if os.path.exists(fcxml) :         etxml = alter(fcxml,edDate,publishDate,createDate)         gp.MetadataImporter_conversion (etxml,fc)         print fc,"updated"     else :         pass         print etxml,"not found"         lstTab = gp.ListTables("*") for tab in lstTab :     # print tab     tabxml = metasrc+"/"+tab+".xml"     if os.path.exists(tabxml) :         etxml = alter(tabxml,edDate,publishDate,createDate)         gp.MetadataImporter_conversion (etxml,tab)         print tab,"updated"         gp.AddMessage(tab+" updated")     else :         pass         print etxml,"not found"         gp.AddError(etxml+" not found") ... View more

I have a rule that any process that takes longer than a "cup of coffee" should be interrupted and a better way found. Multiple ring buffer is a particularly difficult process because it is a script that repeats the problem many times. You problem is likely to be inappropriate data. Can you get a single buffer to work? Buffers can be miss-applied if you have the wrong tolerances and over complex source that you are buffering. You can easily end up with millions of intersections that have to be resolved, and eventually after an hour the process runs out of memory. There is a good presentation by Dale Honeycutt and al "Overlay and Proximity" technical workshop http://proceedings.esri.com/dvd/uc/2009/uc/tws/workshops/tw_715.pdf that explains what happens in a buffer_dissolve. Start with a very simple shape and get that working. Experiment with simplifying the source first and make sure it finishes in seconds. Then try the multiple buffer. ... View more

I am told that there was a deliberate decision to not create one for 10.0. We are supposed to use the intellisense for syntax  and the help documentation. I miss the 9.3 diagram, I found it a helpful quick reference. Maybe we can set up a request on ideas.arcgis.com and vote for it for 10.1? ... View more

You are reinventing ArcGIS in a Python script. There is rarely a need to loop through a set of features to do a selection. You will be running out of memory in double quick time. It is tempting to run a cursor with a few hundred records, but calling a tool while in a cursor is also very slow. I am not sure what your process is, but there must be a command in the toolbox that would do what you appear to be wanting to do with one overlay command that will attach the respective objectids as an attribute in one fast step. I have been doing a similar overlay process this week with 700,000 records on a route system made from 500,000 records and it completes in around 7 minutes. I did avoid CalibrateRoutes and used a modified PointsToPolylines to build a routesystem directly. ... View more

             pnt = part.Next()         partnum += 1 This doesn't look right either, setting a point object to the next part. ... View more

Vertices in Parts are not properly pythonic iterable. You have to use a next() operator or wrap the function in an iterator. See my example on the resources for stepping through parts to find and delete donuts in polygons. The same applies to polylines. http://resources.arcgis.com/gallery/file/geoprocessing/details?entryID=C4E10FE5-1422-2418-A06D-33952BB8D1D7 ... View more