Latest Contributions by KimOllivier

Try the elementTree module for handling XML. It is much easier to use to extract out elements into Python dictionaries and lists. ... View more

I suspect you are running out of memory. Thirty six fields is not a lot, but what sort of fields are they? You did not say how many records you are joining. If you do the spatial join with a small sample does it work? Most join operations fail for me if there are over a million records as a rule of thumb. If you have loaded tables from Excel each field defaults to 255 text chars, even if it is supposed to be a few chars, and numeric fields are double. Try using schema.ini in conjunction with a CSV dump to control the schema when loading, or load more explicitly using a script. Maybe look at tuning your system. Defrag, move/delete/redefine your swap space to be huge and static by making min size == max size Do you have enough spare scratch space, 30% of the disk is recommended. How much free memory do you have before starting? Overlay processes in workstation used to demand 13 times the size of the source coverages, I suspect the same ratio still applies. Have you made your scratch workspace a file geodatabase? Otherwise you have a 2GB limit as shapefiles. A good test is if any single process is taking longer than a cup of coffee, then you have asked for the impossible, and you are probably invoking endless page faults that will end up crashing hours later with no result. So interrupt it and find a better way. ... View more

That is a bit clumsy, given the time to open ArcMap. Oh well, now I know. I was hoping for a menu dropdown that could refresh the AddIn cache without exiting ArcMap. Like Python does it for modules with reload(module) ... View more

Here is a Keyfile selection tool. Rewritten from Bruce Harold's 9.3 version for 10.0 using sets instead of FIDSet which no longer works. My samples were two polygon featureclasses with 2.4 million records in different file geodatabases. I selected 630,000 in one featureclass and found the corresponding records in the other. Total elapsed time to read each featureclass and create a new selection in ArcMap took 6.5 minutes I must catch up with the latest incarnation of ArcScripts.... Anyway it is short enough to paste in, and I will attempt to add a toolbox with a dialog If you want to detect changes in attributes, or deleted records then the same process can be done on a dictionary instead of just a single key. Create a dictionary indexed by the key of all the attributes that need to be compared and convert the keys to a set as before. When you have the intersection, then you can compare the dictionaries for a similarly fast operation. I use this to create  del/new/change sets of layers. Even geometry can be tested if you are creative about using a proxy for the shape such as a rounded area, length or xy pair of coordinates for poly,line,point featureclasses. # Author: ESRI (#5588) # Date:   April 27th 2009 # Purpose: This script applies an SQLquery to a layer or table view using the "IN" #          operator and a key value set taken from another layer or table view field. #          The output is a new selection in the same layer or table view. #          The process is like the KEYFILE reselection option in ArcInfo Workstation. # edit by Kim to fix selection bugs # 17 June 2009 # 12 August 2009 FIDSet broken at 9.3.1 ?? # 25 April 2012 10.0 completely redesigned using sets, 10.x compatible, faster # Kim Ollivier [email protected] import arcpy import datetime ts0 = datetime.datetime.now() arcpy.env.overwriteOutput = True try:     # raise Exception # uncomment for debugging in PythonWin     inLayer = arcpy.GetParameterAsText(0) # Get the input layer or table view to be subqueried     inField = arcpy.GetParameterAsText(1) # Get the input object subquery field     keyLayer = arcpy.GetParameterAsText(2) # Get the keyfile layer or table view     keyField = arcpy.GetParameterAsText(3) # Get the keyfile field     selOption = arcpy.GetParameterAsText(4) # Get the selection option defaults to NEW_SELECTION except: # debugging     inLayer = "e:/project/BTR/view/Geocoding_business2.lyr" # Get the input layer or table view to be subqueried     inField = "ID" # Get the input object subquery field     keyLayer = "e:/project/BTR/view/Geocoding_business.lyr" # Get the keyfile layer or table view     keyField = "ID" # Get the keyfile field     selOption = "NEW_SELECTION" # Get the selection option defaults to NEW_SELECTION    # Build two sets of keys setInKey = set([row.getValue(inField) for row in arcpy.SearchCursor(inLayer)]) setKey = set([row.getValue(keyField) for row in arcpy.SearchCursor(keyLayer)]) # only uses selected records # make the common set (so fast cannot be measured! setCommon = setInKey.intersection(setKey) print "Common",len(setCommon) arcpy.AddMessage(str(len(setCommon))+ " common records") if len(setCommon) > 0:     # build an SQLquery "[inField] IN (1,2,3....)"     delimitedInField = arcpy.AddFieldDelimiters(inLayer,inField)     subQuery = delimitedInField + " in "+ str(tuple(setCommon)) # creates a perfect SQL style list     # Update the layer with the new selection     arcpy.AddMessage(selOption+" "+inLayer)     result = arcpy.SelectLayerByAttribute_management(inLayer,selOption,subQuery)     outLayer = result.getOutput(0)     arcpy.SetParameter(5,outLayer) else:     arcpy.AddWarning("No common records, so no selection possible") print "Total",datetime.datetime.now() - ts0 keyfile tool ------------ This tool is a replacement for the ARC/INFO workstation SELECT <cover> KEYFILE <cover> <keyitem> It was originally 9.3 proposed by Bruce Harold using FIDSet on the Describe tool but this now fails and seems obsolete. This is because the geoprocessing tools finally work on selected records in layers. A better method is to use the new Python set data structure for comparisons. The set comparison is so fast thst it can hardly be measured with datetime. You still have to create the sets with a SearchCursor, but that may be a lot faster at 10.1 with the new da module. If you have a very large number of records you must have the key field indexed to ensure that the selection draws in ArcMap in a reasonable time because the SQL expression is a list, not an equation. This tool does not check for indexes. There does not seem to be a limit on the number of records to be compared, but the largest size file for testing has been to select 2,500 records from another file of 2.4 million records which took 5 minutes to complete. Selection 630,000 records from a set of 2.4 million and run against another earlier set of 2.4 million tool 6.5 minutes. **Warning** Field Index names in ArcGIS are limited to 16 characters long, so you have to be careful just concatenating _idx to the field name. The nearest equivalent tool is MakeQueryTable. This tool has a number of limitations not in this keyfile tool. The tables being compared can be in different databases and even different types of database. Since only a selection is changed, all layers are in memory for speed. Make sure that you have suitable local memory and scratch space. Both keys have to be the same basic type for the SQL query to be valid. Either both numeric or both text. Limitations ----------- There is no validation or error trapping except if there are no common records. It is recommended that the keys are indexed for large tables. The toolbox is 9.3 format but the script is now written using 10.0 syntax with arcpy Kim Ollivier [email protected] www.ollivier.co.nz 25 April 2012 (ANZAC Day) ... View more

Imbedding any geoprocessing tool inside a cursor will always be too slow. They are not designed to be used that way, only on whole datasets at once, any looping is built into the tools and uses indexing and other optimisation. You need to find a tool that does the comparison inside the tool. 1. MakeQueryTable springs to mind, but first you would have to copy the table so they are both in the same database. Nice and simple but not very flexible. You would make the comparison in an SQLquery expression. 2. Find the differences in a python set. This is the fastest by a long way (milliseconds for the whole thing). The catch is that you still have to read the tables once each to make the sets and write out a selected set. But it will still be less than minute or two. ... View more

Kevin,  I see you still don't have any answers for this. Can you provide more information with your question? For example: the whole piece of code you are using, the error you are receiving...  I am not familiar with any of the new features in 10.1, but in 10 you get values from inputs (including, I presume, comboboxes) via   arcpy.GetValueAsText(). Or are you using something like   Tk? If you have not got 10.1 and are not trying to use the new 10.1 only AddIns module for Python then you cannot help. I can't work it out either. I assume we have to understand how namespaces in classes work. We have to find the name of the instance of the class which might be that ID that we fill in at the wizard screen. ... View more

The built-in compare tool seems to be designed to check two featureclasses after a copy and edit. I cannot see the point of stopping after find the first difference without saying which one. What we really need is a compare using a static key (which used to be called the user-id). The fastest way to do this uses the Python Set operations which run in milliseconds. But first you have to create the set for each of the featureclasses which takes some time with a SearchCursor. This only takes a few seconds. Finally you need to create a selection or new featureclass containing the deletes/adds/changes. This is easy enough but a bit slower, taking perhaps minutes. Use an SQL query using IN (list_of_keys) You need to keep track of the primary key for each feature. Do this using a Dictionary keyed with the key and containing all the attributes that you want to compare. Comparing shapes is difficult because it requires them to be in a standalone object so you have to unpack the shape field into a Python string. Instead use the area or length as a proxy for the shape and round it to a suitable precision eg int(round(row.shape.area)*100)) #NAME:  delta_road.py #AUTHOR:  Kevin Bell #EMAIL:  [email protected] #DATE:  20071207 # edited by Kim Ollivier #PURPOSE:  create adds/deletes layer files by comparing 2 point #          feature classes shape and attributes. If the shape has #          not changed, but any of the attributes have, the feature #          will show as a delete, and an add. #NOTE:  buildDict method has hard coded primary key and attribute names. # mods : # use sets for del/add # write out featureclass differences # do not use centroid for lines, they have changed in a year! suggest length # do not use shape for non-points # dont forget to index keys for definition layers def buildDeltaLayers(inFC1, inFC2):     '''build an adds and deletes lyr for a given chrono FC '''     print inFC1, inFC2     d1 = buildDict(inFC1)     d2 = buildDict(inFC2) #   find set differences and intersections to reduce dictionary comparisons     startTime = time.clock()     s1 = set(d1.keys())     s2 = set(d2.keys())     sNew = s2 - s1     print "New",len(sNew)     sDel = s1 - s2     print "Del",len(sDel)     sInt = s1.intersection(s2)     print "Int",len(sInt)     print "sets done"     stopTime = time.clock()     elapsedTime = stopTime - startTime     print "elapsed time = " + str(round(elapsedTime, 1)) + " seconds"     # could weed the dictionaries but will that be slower than just using the intersection set?     compareList = valuesChanged(d1,d2,sInt)     print "Changes",len(compareList)     stopTime = time.clock()     elapsedTime = stopTime - startTime     print "elapsed time = " + str(round(elapsedTime, 1)) + " seconds"     changes = len(compareList)     if changes > 0 and changes < 50000:         print "Changes to be written",changes         makeLYR(inFC2, compareList, "rdchg")             else :         print "no changes, or too many to be believed",changes     # must create delete layer first     # to allow new features to be filtered for ID renumbering     if len(sDel) > 0 :         makeLYR(inFC1, sDel,"rddel")     if len(sNew) > 0 :         makeLYR(inFC2, sNew,"rdnew")     def valuesChanged(dict1, dict2,sBoth):     '''get a list of keys from one dict if a corresponding dict's values are different'''     ##    outList = [key for key in set(dict1.keys() + dict2.keys()) if dict1.get(key) != dict2.get(key)]     outList = [key for key in sBoth if dict1.get(key) != dict2.get(key)]     return outList def buildDict(inputFC): #-----BEWARE OF HARDCODED PRIMARY KEY AND ATTRIBUTES BELOW!!!!!     '''Build a dictionary of the primary key, and its fields'''     startTime = time.clock()     d = {}     cur = gp.SearchCursor(inputFC)     row = cur.Next()     while row:         # only need to check primary keys and shape         # oops row.shape.centroid always fails - why?? because its only an object reference         # you have to get out the coordinates and put it in a real Python object         d[row.GetValue(pk)] = [row.name.upper().strip() ]#, round(row.shape.length)]         # d[row.GetValue(pk)] = [round(row.shape.length)]         row = cur.Next()     del cur     print inputFC,     stopTime = time.clock()     elapsedTime = stopTime - startTime     print "dict created, elapsed time = " + str(round(elapsedTime, 1)) + " seconds"     return d def makeLYR(fc, inList, outLyrName): # BEWARE OF HARDCODED PRIMARY KEY BELOW     '''Given a list, return a LYR file'''         startTime = time.clock()     wc = str(tuple(inList))     print outLyrName,len(inList)     whereclause = pk+" IN " + wc # <----IF DATA ISN'T FILE GDB, YOU MAY NEED QUOTES/BRACKETS     # print whereclause     gp.MakeFeatureLayer_management (fc, outLyrName, whereclause)     print outLyrName,"ORIG layer count",gp.GetCount(outLyrName).GetOutput(0)     # remove changes of just RCL_ID     if outLyrName == 'rdnew':         print "RDDEL layer count",gp.GetCount("rddel").GetOutput(0)         gp.SelectLayerByLocation_management(outLyrName, "WITHIN", "rddel","","REMOVE_FROM_SELECTION")         print "Made NEW layer",outLyrName,round(time.clock() - startTime)," seconds"         print outLyrName,"NEW layer count",gp.GetCount(outLyrName).GetOutput(0)     gp.RefreshCatalog(delta_gdb)     print outLyrName,"layer count",gp.GetCount(outLyrName).GetOutput(0)     layerfn = os.path.dirname(delta_gdb)+"/"+outLyrName +".lyr"     if gp.Exists(layerfn) :         gp.Delete(layerfn)         ##gp.SaveToLayerFile_management(outLyrName, layerfn)     # print "saved layer",round(time.clock() - startTime)     deltaFC = delta_gdb+"/"+outLyrName     # Warning, featureclass MUST be indexed on PK!     if not gp.Exists(fc) :         print fc,"not found"     ## gp.AddIndex_management(fc,pk,fc+"_"+pk+"_idx")     if gp.Exists(deltaFC):         gp.Delete(deltaFC)     gp.CopyFeatures_management(outLyrName,deltaFC)     stopTime = time.clock()     elapsedTime = stopTime - startTime     print outLyrName,"elapsed time = " + str(round(elapsedTime, 1)) + " seconds" #---------------------------------------------------------------------- print "-----------  delta road  ----------------" import arcgisscripting, time,os,sys gp = arcgisscripting.create(9.3) gp.overwriteoutput = True try :     last_gdb = sys.argv[1]     current_gdb = sys.argv[2]     delta_gdb = sys.argv[3] except :     last_gdb    = "E:/road_jan2009.gdb"     current_gdb = "e:/crs/mobile/mobile.gdb"     delta_gdb   = "E:/crs/road_roadname.gdb" if not os.path.exists(delta_gdb) :     gp.CreateFileGDB_management(os.path.dirname(delta_gdb), os.path.basename(delta_gdb)) gp.Workspace = delta_gdb gp.OverwriteOutput = 1 print for laydef in [["road","RCL_ID"]] :     layer = laydef[0]     pk = laydef[1]     startTime = time.clock()     print "Start"     old = last_gdb+"/"+layer     new = current_gdb+"/"+layer     ready = True     if not gp.Exists(old) :         print old,'not found'         ready = False     if not gp.Exists(new) :         print new,'not found'         ready = False     if ready :           buildDeltaLayers(old,new)         msg = "Your del/change/add fc's are in:"         gp.AddMessage(msg)         print str(gp.workspace)         gp.AddMessage(str(gp.workspace))         stopTime = time.clock()         elapsedTime = stopTime - startTime         print layer,"total elapsed time = " + str(round(elapsedTime, 1)) + " seconds"     else :         print "not ready" # del gp print "-------------------------------------------------------------" gp.AddMessage("finished") ... View more

There is only one text source for a Python script, but after an edit it is compressed into a zip bundle *.esriaddin using makeaddin.py I am then supposed to double click on it to install it into the well known location in MyDocuments. I assume this  is different to addins created in java and .NET where the source is in an IDE. The Python wizard is not an integrated environment. So if you are creating a .NET addin you can edit your (non-python) script while ArcMap stays up? It all works just fine if I close and reopen ArcMap, but what a pain! It seems that the addin is loaded into memory when ArcMap starts and no amount of deleting refreshes the image of the addin. ... View more

The Near tool has an option to record the direction angle so you could find the distance to all eastern states. Select the states to the east of the points could be done by comparing the longitude of the point and the centroid of the state, or  maybe the extent of the state. If you make a list of these they could be input into the Near tool. They might need to be in separate featureclasses, one for each state. So while you are  there, project the states into projected units  that you require. At 10.x you could create geometry objects, one for each state, to hold them all in memory instead of 50 featureclasses. I am trying to think of a way of avoid a loop around each point. Another approach might be to create a line eastwards from each point, intersect the whole lot at once and find the first intersection from each point. ... View more

Yes, it is possible, but you will have to write your own tool. I did this for some MrSID image extents way back (2004) with ArcGIS 8.3 using the Python SAX module to parse the metadata XML file. You have to define a handler function for each tag that you want to extract using SAX. I have even found the code, which you could have if you really want to go down that route. I see that Open Office still use this method which is very fast for large XML files typical of word documents. I learnt how to do this from an O'Reilly book: Python & XML It is easier to use the newer ElementTree module. You will need to export the metadata to a file to be able to parse it. ... View more

Instead of using a codeblock in CalculateField, open an updateCursor in a Python script and put the code in there. You have to extract out the shape as an object, rebuild a new object and update the shape field. A script can be much easier to write, with exception traps.         # POS transect with measures                 posArray = arcpy.Array()         for transit in dPOS.keys():             for mrk in dPOS[transit]:                 ta,tb = trans(mrk[0],mrk[1])                 posArray.add(arcpy.Point(ta,tb,mrk[2],mrk[3]))             rowPOS = curPOS.newRow()             rowPOS.mark = transit.upper()             rowPOS.plotid = plotid             rowPOS.plotmark = plotid + transit.upper()             rowPOS.shape = posArray             curPOS.insertRow(rowPOS)             posArray.removeAll() In this snippet I am looping around a dictionary of coordinates (dPOS), transforming them and building an array of points in posArray for the polyline feature. I am creating the lines, not updating. Note the four values ta,tb,mrk[2],mrk[3] which are x,y,z,m values. In your case you would be extracting the current values from the shape using an update cursor with something like Here is a more complex example showing how you have to iterate over the parts in a polyline to get to the vertices. # fillDonut.py # frequent forum request to fill holes # use an Update Cursor to edit polys with donuts # 22 May 2010 # Kim Ollivier [email protected] # (c) Creative Commons 3.0 (Attribution) # Method # ------ # Polygon features are stored as: # [[pnt,pnt,pnt,pnt,,pnt,pnt,pnt,pnt][pnt,pnt,pnt,pnt]] # Parts are separate lists inside the feature object array # null point in a part indicates donut ring(s) follow # So just step through and truncate at first null pnt per part, #   re-assemble parts into a polygon # Update feature if a donut found, otherwise continue # Input is a layer so can use selection in Arcmap # could also put on menu at ArcMap 10 # can be back-ported to 9.3 # 4 June 2010 # added iter lambda function for part array because isn't interable # cleared polyOuter buffer earlier to fix logic for multiple parts import arcpy,sys def filldonut(inlay) :     '''Edits a layer in-place     fills all features with donuts     requires a layer to honour     selected features     '''     desc = arcpy.Describe(inlay)     shapefield = desc.ShapeFieldName     rows = arcpy.UpdateCursor(inlay)     n = 0     polyGeo = arcpy.Array() # to hold edited shape     polyOuter = arcpy.Array() # to hold outer ring     for row in rows:         feat = row.getValue(shapefield)         qInterior = False         for partNum in range(feat.partCount) :             part = feat.getPart(partNum)             for pt in iter(lambda:part.next(),None) : # iter stops on null pt                 polyOuter.append(pt)             if part.count > len(polyOuter) :                 qInterior = True             polyGeo.append(polyOuter) # reassemble each part             polyOuter.removeAll() # ready for next part         if qInterior : # in any part of this feature, otherwise skip update             n+=1             row.setValue(shapefield,polyGeo)             rows.updateRow(row)         polyGeo.removeAll()     del rows,row     msg = "Features with interior ring filled "+str(n)     arcpy.AddMessage(msg)     arcpy.AddMessage("Refresh view to see donuts disappear!") #---------------------------------------------------- if __name__ == '__main__' :     inlay = sys.argv[1] # arcpy.GetParameterAsText(0)     desc = arcpy.Describe(inlay)     if desc.shapetype != 'Polygon' :         arcpy.AddError(inlay + " Not a polygon class")     else :         filldonut(inlay)         arcpy.RefreshGraphics() (Unfortunately it is possible to fool this code with complex nested imbedded holes in multipart polygons. but it works for simple cases.) ... View more

I am using dictionaries to update tables instead of a join more as well. I tried refactoring my clumsy lines to use the oneline list comprehension but it turned out to be marginally slower. 222565 dictionary count 0:00:46.594000 222565 dictionary count 0:00:48.125000 I note that you do not bother to specify a subset of fields when opening the cursor. If you have a lot of fields it apparently helps a lot to only list the relevant fields for the calculations. Not so easy to generalise I suppose, but it may help with memory management too. Has anyone done some tests on the 10.1 da module that has  rewritten cursors? Maybe we will not need dictionaries after all. ... View more

But your data does contain some "unexpected data". The error message says so: AttributeError: <unprintable AttributeError object> <type 'exceptions.UnicodeEncodeError'>: 'ascii' codec can't encode characters in position 29-37: ordinal not in range(128) It even tells you where the characters are. It is very easy these days to enter unicode by mistake from the keyboard or a web page cut and paste. Not all Python functions handle unicode as if they were ascii. If you have tried to convert a unicode string that contains a character other than the ascii set then you need to encode/decode the value to an equivalent ascii value. Just doing >> str(u'\u0100') Traceback (most recent call last):   File "<interactive input>", line 1, in <module> UnicodeEncodeError: 'ascii' codec can't encode character u'\u0100' in position 0: ordinal not in range(128) will upset Python >>> str(u'\u0100'.encode('utf-8')) '\xc4\x80' And this will probably upset the reader, but it will explain the problem. Here is what I had to work out to add macrons to a database # unicode.py This Python file uses the following encoding: utf-8 import encodings.utf_8_sig comment = """ <HEAD> <META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=iso-8859-1"> To make your pages display properly in Unicode, you need to change this to: <HEAD> <META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=UTF-8"> import codecs f = codecs.open("file","w",encoding='utf-8') f.write(stuff) f.close() """ ##&#256 = A- ##&#274 = E- ##&#298 = I- ##&#332 = O- ##&#362 = U- ##&#257 = a- ##&#275 = e- ##&#299 = i- ##&#333 = o- ##&#363 = u- dMacron = { "A" : u'\u0100', "a" : u'\u0101', "E" : u'\u0112', "e" : u'\u0113', "I" : u'\u012A', "i" : u'\u012B', "O" : u'\u014C', "o" : u'\u014D', "U" : u'\u016A', "u" : u'\u016B' } f1 = open("e:/project/training/unicode/unicode_sample.txt","w") print # read-only f1.encoding # = "utf8" # print >> f1,dMacron # cannot do this to ascii encoder mO = u"\N{LATIN CAPITAL LETTER O WITH MACRON}" word = "Maori".replace("o",dMacron["o"]) print >> f1,word.encode('utf-8') # PyUnicode.__print__ word.encode('utf-8') lstKey = dMacron.keys() lstKey.sort() for k in lstKey : print >> f1,dMacron .encode('utf-8'), print dMacron , ##mA = u'\u0100' ##ma = u'\u0101' ##mE = u'\u0112' ##me = u'\u0113' ##mI = u'\u012A' ##mi = u'\u012B' ##mO = u'\u014C' ##mo = u'\u014D' ##mU = u'\u016A' ##mu = u'\u016B' # word = "Maori".replace("o",dMacron["o"]) mO = u"\N{LATIN CAPITAL LETTER O WITH MACRON}" word = "MAORI".replace("O",mO) print >> f1,word.encode('utf-8') print word f1.close() print "success" I once reinstalled the entire operating system on a Sun 3, including printer drivers when the new postscript printer stopped working. It had stopped because the printer had run out of paper and the Sun driver did not understand postscript error messages. 🙂 I don't ever reinstall a software package when I hit a bug now, life is too short. ... View more