Latest Contributions by KimOllivier

It just means that you have to use Python to load the XY.txt table to a geodatabase table first. There are plenty of tools for that such as the CSV module. There are plenty of other issues with Windows XY events when trying to write it out to a featureclass. ... View more

It must be the different default environment settings. We need more detail of your setup to guess the issues. Where are you creating the featureclass? SDE, versioning, local geodatabase? Have you set your scratch workspace too? These will all have defaults when you add a tool What version and service pack? Why can't you use a template? ... View more

I finally reproduced your errors, and some more using your data. There seems to be something about data in proper metric units that works ok but not US Feet 😉 So I thought that it must be a data corruption. I have not found a fault that would explain that, and now after fiddling around have ended up with some similar results. I think the fundamental problem is similar to one I hit on another project where the shape object is only referenced in memory when it is extracted with a cursor. Later when we try to access it, there are no methods because it is not a proper Python object. Try doing a dir(shape). If it was a compliant Python object it would return all the methods listed on the Geometry object. Even with interactive mode you cannot do any useful request on a shape retrieved from the dictionary, you get a crash, even when editing a script and invoking a method. It only works while in the cursor. The only solution at 9.3 is to extract out all the vertices into an object that Python can use, like a string based GeoJSON dictionary! (At 10 this is practical because there is a function to convert a shape to a GeoJSON string.) See Bruce Harold's example trying to comprehensively detect the changes between two feature datasets on the Resource Centre. Although the script is in 10 format, the same applies at 9.3. My shortcut solution for a difference machine was to test the perimeter, length or centroid instead of every vertex because if even one changes it will affect these values. Can you save gp.Geometry objects in a dictionary? It appears to be intended at 10, but I know there were issues at 9.3. So I don't think that there is much point in proceeding with more testing at 9.3. ... View more

The code example worked for me. I could not get it to fail on a complex polygon featureclass such as the pine forest landcover for a small country. I changed the test to a random shape and printed all vertices and handled possible donuts. I did not compare the vertices with a normal cursor, but they look plausible. There are a few more convenient touches at 10 for iterating over parts and of course the spatial operators but otherwise its the same for 9.3 or 10. The featurecount is 47,713 and the script took 14 seconds. THIS CODE LOADS A FC INTO A PYTHON DICTIONARY IN RAM # http://forums.arcgis.com/threads/17956-Using-centroid-geometry-to-get-values-of-intersecting-feature#post60607 # edits by Kimo import arcgisscripting, sys import random,datetime,pprint start = datetime.datetime.now() gp = arcgisscripting.create(9.3) #THE INPUT FC fc = "d:/workspace/datetest.gdb/pineforest" #Process: Build a field list and a mini dictionary to look up the field order (aka index) given a field name fieldList = gp.listfields(fc) fieldIndex = 0 fieldNameDict = {} for field in fieldList:     fieldNameDict[field.name] = fieldIndex     fieldIndex = fieldIndex + 1 shapeFieldName = gp.Describe(fc).shapeFieldName #  fix missing definition pprint.pprint(fieldNameDict) #Process: Load the fc into a dictionary fcDictionary = {} oidFieldName = gp.describe(fc).oidfieldname searchRows = gp.searchcursor(fc) searchRow = searchRows.next() while searchRow:     oidFieldValue = searchRow.getvalue(oidFieldName)     fcDictionary[oidFieldValue] = []     for field in fieldList:         fcDictionary[oidFieldValue].append(searchRow.getvalue(field.name))     searchRow = searchRows.next() del searchRow,searchRows print print "Size of dictionary",len(fcDictionary) for sample in range(10):     print     idx = random.randint(0,len(fcDictionary))     #Print the values     print fcDictionary[idx]     #Print a name field value     print fcDictionary[idx][fieldNameDict["LCDB1NAME"]]     #Print the extent rectangle     print fcDictionary[idx][fieldNameDict[shapeFieldName]].extent     #Print ALL the coordinate pairs including parts and interior rings     shapeObj = fcDictionary[idx][fieldNameDict[shapeFieldName]]     for p in range(shapeObj.partCount):         part = shapeObj.getPart(p)         pnt = part.next()         while pnt:             print p,pnt.X,pnt.Y             pnt = part.next()             if not pnt:                 pnt = part.next()                 if pnt:                     interiorRing = True                     print "Hole" print "Completed",datetime.datetime.now() - start [12139, <geoprocessing describe geometry object object at 0x16A388D8>, 66.0, u'Pine Forest - Closed Canopy', 66.0, u'Pine Forest - Closed Canopy', u'no', 64.70760971, 182.0, u'Rodney District', 16.0, u'Auckland Region', 11067.655349413637, 646803.66190707707] Pine Forest - Closed Canopy 1811685.23719998 5594544.72269999 1811781.95419998 5594599.59829999 NaN NaN NaN NaN 0 1811685.2372 5594544.9228 0 1811688.5206 5594544.9241 0 1811733.3278 5594561.747 0 1811781.9542 5594599.5983 0 1811780.0872 5594596.3042 0 1811757.3881 5594569.9774 0 1811726.5661 5594549.7439 0 1811689.9306 5594544.7227 0 1811685.2372 5594544.9228 Completed 0:00:14 ... View more

But is it corrupted? This works for me: # create a dictionary of centroids sr = gp.Describe(fcAU).SpatialReference cur = gp.SearchCursor(fcAU,"land = 1") row = cur.next() dAU = {} pt = gp.CreateObject("Point") while row:     pt = row.shape.centroid     au = row.AU_NO     dAU[au] = pt     row = cur.next() del row,cur print "OD records",len(dAU) ... View more

E00 files are ASCII dumps of a coverage. You have to use a special tool to import them that requires an ArcInfo licence. Coverage Tools>Conversion>To Coverage>Import From Interchange File But this will load into a coverage which needs a second export to a geodatabase ArcView 3.x has a simple standalone import utility too. Or use any open source tool. They load them without fuss into shape files. ... View more

If you don't need measures then make sure that they are disabled, or drop them. Sometimes the range of valid measures also defaults to much smaller than you might have by default. When you define new featureclasses in 'high precision' then always define a projection to let ArcGIS calculate suitable limits for xy and m values. Have a look at your coordinate system extents to see if they are unexpectedly small. ... View more

import win32api win32api.MessageBox(0,"msgbox", "title") win32api.MessageBox(0,"ok cancel?", "title",1) win32api.MessageBox(0,"abort retry ignore?", "title",2) win32api.MessageBox(0,"yes no cancel?", "title",3) You can use the Windows API to popup messages and store the response. Have a look in the Pythonwin help for additional tools. You can also use Tkinter (standard library) to create a full dialog. See the example to create a scrolling list to supplement the data driven pages tool on the Resource Centre. Python Selection Script with Tkinter User Interface http://resources.arcgis.com/gallery/file/geoprocessing/details?entryID=F0B54854-1422-2418-A0F6-37042FE94472 ... View more

It's a good question we all ask ourselves every day. I have the "Cup of Coffee Rule". -If any single process takes longer than a cup of coffee, then interrupt the process and Find a Better Way. I don't bother to run any process for days, it will not be likely to produce anything useful, as you have found. Too much is done in memory, then virtual memory, then paging to disk and then there is none left to write out the results. This happens to me when I try to create a large XY event table from a route, but it only takes an hour to crash! I watch the CPU spinning with no I/O. Its obviously a very poor algorithm. I'm not quite sure why you are using a spatial join, presumably using the interactive join tool? Have you calculated the likely size of that output? 250,000 x 250,000 = 6.25E10 ! Yes, it would be a good idea to limit the comparisons to a buffer if you can. It's a worry that you cannot get a process to work on a sample of 35,000 polygons. That should work in minutes, so you are missing something.To just give you a benchmark time I ran a union between 33,000 parcels with a 428,000 polygon landcover dataset as the hardest process I could think of with similar sizes. That gives me a parcel-landcover combined set of  polygons with all attributes for further statistical summary in 5 minutes. Tune your PC. I have a similar machine, but I know that 8 CPU's are irrelevant, only one is used. Do you have heaps of virtual memory on a separate disk in a contiguous partition, or even better a dedicated disk for it? Do you have a local scratch workspace defined as a file geodatabase (not just a folder). Have you got a RAID 0 array running at 10,000RPM? Is there plenty of free disk, unfragmented. Since the process is equivalent to an old clean, you probably still need 13 times the space of the source before starting. Is all your data in local file geodatabases, not on a network or in SDE? Dragging data across to do analysis is likely to be unsuccessful. Tune your data. Have you run a Repair Geometry on the data first? Check for uniform size of polygons. If there are some very large ones, split them up. At 10 there is a dice tool to automate this. Do you have multi-part polygons? These are a great evil for geoprocessing. Explode them first. Be reasonable with your requests. ArcGIS does not often warn you if you have not thought through the task. You might get a hint that a join does not have an attribute index, but not that you do not have a spatial index or simply multiplying the two sets to give you a predicted output size. If I can recast my problem into a simpler one which does not requires polygon splitting then I do. eg Will a point in polygon do? Extract one layer as centroids. Not the same, but maybe good enough. Finding a better way would be to use another Python module or another tool such as FME. ... View more

There is a better way to sort points for labelling. The technique uses a Peano curve to sort the points into a snakes and ladders type of route where the next number is always near the last. This is similar to the way that we used to number labels on polygons when we used to digitise by hand! The Peano Curve is an option on the ArcGIS 10 sort (data management) tool. It basically comes up with a single index for the x,y pairs that is a very useful figure to sort a point cloud so that you can easily find any point if they are sequentially numbered. ... View more

Maybe it is a bug? NIM060284 - The Point class constructor arguments m and z are reversed. fixed in 10 SP1 ... View more

You can directly calculate DMS using the calculate geometry interactive tool on a table in ArcMap. The trick is to create a TEXT field to calculate into and you will be given more output format options. Doesn't work for a script, I know. ... View more

How hard can it be?! The maths to calculate the distance between two points could use good old Pythagoras once you have the point coordinates in a list. The numbers may balloon out because the product of the two sets will often be large. To reduce this, use the Python module itertools  in 2.6 to generate the permutations of the point pairs n*(n-1)/2 so that you don't measure to-from as well as from-to (n*n). To make running through lists repeatedly more efficient, load the ID,x,y for each point into a Python structure such as a dictionary. Numpy has a function to do the maths even faster than Pythagoras implemented in a script. distance = numpy.linalg.norm(x-y) ... View more

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