Latest Contributions by KimOllivier

Ok, I didn't think I would get away with that brief specification. A sample script takes 3 seconds to iterate 100 times. I have shown one case creating a layer of the random selection. You could put this back into a table of contents. # randomselection.py # 19 May 2011 # kimo@ollivier.co.nz # # Input: parcel polygon featureclass #        item to sum (area) #        threshold (check min, mean, max values for reasonable figure) # Assumes values always less than threshold or you get a null list # Output: Total sum and layer with selected parcels # set up to allow iteration # uses ArcGIS 10 but could use 9.3 import arcpy import sys import os import random import datetime try:     fc = sys.argv[1]     item = sys.argv[2]     threshold = sys.argv[3] except:     # debugging and standalone operation     fc = "c:/project/forum/mobile.gdb/parcel"     item = "area"     threshold = 100000.0 # square metres def getdict(fc,item):     ''' Make a dictionary     once for the sampler     '''     dParcel = {}     for row in arcpy.SearchCursor(fc,"featcode = 'parcel'"):         if item.lower() == 'area':             dParcel[row.objectid] = row.shape.area         else :             dParcel[row.objectid] = row.GetValue(item)     return dParcel def getsample(dParcel,threshold):     '''Given a dict of values     return a random selection     less than the threshold.         Shuffles the list and pops off     the ids until threshold reached     ie sampling without replacement     '''     total = 0     lstObj = []     lstKey = dParcel.keys()     random.shuffle(lstKey)     while total < threshold:         nObj = lstKey.pop()         nVal = dParcel[nObj]         total = nVal         if total > threshold:             break # do not exceed threshold         else :             lstObj.append(nObj)             total += nVal     return (total,lstObj) # --------------- main --------------- print start = datetime.datetime.now() print start # # make a dictionary dParcel = getdict(fc,item) print "records",len(dParcel) # sample the dictionary total,lstObj = getsample(dParcel,threshold) print "Total: ",total print "No of parcels: ",len(lstObj) print "objectid: value" for ob in lstObj:     print ob,":",dParcel[ob],",  ", print # make up an sql query to select the features sqlQuery = "OBJECTID in (" + ",".join([str(id) for id in lstObj]) + ")" print sqlQuery arcpy.env.overwriteOutput = True arcpy.MakeFeatureLayer_management(fc,"fc_lay",sqlQuery) lyrName = arcpy.CreateUniqueName("xx.lyr",os.path.dirname(os.path.dirname(fc))) print lyrName arcpy.SaveToLayerFile_management("fc_lay",lyrName) # Or sample the dictionary 100 times for it in range(100):     total,lstObj = getsample(dParcel,threshold)     print total,len(lstObj)     # finished elapsed = datetime.datetime.now() - start print "\n\nElapsed",elapsed, "Well done!" # positive feedback 2011-05-19 20:40:43.640000 records 7170 Total:  145316.242757 No of parcels:  21 objectid: value 726 : 809.726863177 ,   4530 : 809.419546424 ,   4752 : 1012.70775115 ,   6120 : 809.397052955 ,   3022 : 809.418725814 ,   253 : 986.240771943 ,   4445 : 1132.16238039 ,   7161 : 928.75695316 ,   5747 : 964.13873029 ,   2188 : 864.834770195 ,   4164 : 20.9192428737 ,   187 : 839.573627519 ,   2745 : 811.268061104 ,   2572 : 1002.32513709 ,   1445 : 1381.14243957 ,   2674 : 812.661986759 ,   6925 : 25126.7099296 ,   5479 : 1395.56544711 ,   4791 : 810.812432473 ,   5406 : 848.317697436 ,   6543 : 72658.1213786 ,  OBJECTID in (726,4530,4752,6120,3022,253,4445,7161,5747,2188,4164,187,2745,2572,1445,2674,6925,5479,4791,5406,6543) c:/project/forum\xx0.lyr 185921.392647 4 145871.629063 27 154201.63183 1 ... 140524.256299 0 ... 190822.946912 9 193668.483773 23 334292.358346 16 141967.626185 23 15109126.1898 56 ... 413536.709991 2 157130.891363 6 Elapsed 0:00:02.917000 Well done! ... View more

Note to ArcView users (myself included), you won't be able to use the solution found in the link, as this tool references two other toolbox tools (Point to Feature and Frequency) that are not available at the ArcView license level. How and why this functionality was essentially left out is maddening to me, ha. /end rant I keep getting caught by this too. I build a tool for ArcView with a higher licence, then it won't work on my laptop. You don' need to use Frequency. Frequency_analysis is a hangover from the ArcInfo tools. There is an equivalent Arcview level tool called Summarize_analysis that includes all the functionality of Frequency and more! I am sure that Esri did not intend us to have to upgrade just to count features. http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00080000001z000000.htm FeatureToPoint is an ArcInfo tool, but for goodness sake, if you insist on using a viewer for analysis, then get cracking with your own Python script to replace this tool. It isn't  very hard to run a cursor and get out the centroid of each feature to export as a point layer, it's built into the basic toolkit for Arcview. I can't see why it is a higher licence myself, maybe it is to encourage us to learn Python (or keep a third party addin such as the excellent ETGeotools in business). Everything that was in DS Mapbook can be done with data driven pages and a couple of free scripts that are much easier to modify than the DS Visual Basic. After all DS (Developer Sample) was an add-in itself, it was not part of the base. It is no reason not to upgrade IMHO. ... View more

I have a mapbook with 1,440 pages that label neighboring tiles and I do not see a good reason to upgrade to 10 based on what I have tried and have read. Also what is a good way to export random pages, I do not use page numbers but a name instead. jim There is a fine example Python based tool that creates a dialog with map sheet names instead of numbers on the resource pages. I know it is so hard to find these, no arcscripts search now, so here is the whole script unzipped: from Tkinter import * import arcpy # SelectPopup.py sample Python script using Tkinter for a user interface # User selects feature from list, code selects and zooms to feature # Works with file geodatabases and sdc datasets # Use in ArcMap by attaching script to a button for best results # In ArcGIS 10, you can add scripts and models to buttons in the UI without using VBA. # See "Adding a custom tool to a menu or toolbar" in this help topic: # http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//002400000005000000.htm # Replace the values for the 4 variables below to match your data # if you change the variable name, be sure to modify the code # below to match. # Layer on which the selection will be performed Layer = "states" # Table or feature class with field containing unique list of feature names to use for selections # Can be the data source for the layer above or a separate table. Use a separate table containing # only the valid names if there are large numbers of blank or non-unique names or a definition query # on the layer. # SelTable = "D:\Work\Demos\Done\Arkansas10\GeoStor.gdb\NRHP_YrBlt" SelTable = "D:\Work\ESRIData&Maps\mexico\data\states.sdc\states" # Name of field in SelTable and Layer that contains feature names SelField = "ST_NAME" # Title to Appear on Dialog Title = "Select A State" # End of section for entering your data parameters class Application(Frame):                def __init__(self, master=None):         Frame.__init__(self, master)         self.grid()              self.createWidgets(master)     def createWidgets(self, master):         self.yScroll  =  Scrollbar ( self, orient=VERTICAL )         self.yScroll.grid ( row=0, column=1, sticky=N+S )         self.stList = Listbox (self, yscrollcommand=self.yScroll.set)         self.stList.grid( row=0, column=0 )         self.yScroll["command"] = self.stList.yview         #populate list with choices, sorted ascending         rows = arcpy.SearchCursor(SelTable, "", "", SelField, SelField + " A")         for row in rows:             self.stList.insert( END, row.getValue(SelField) )         #add selection and quit button         self.selButton = Button (self, text='Select', command=self.selectFeat)         self.selButton.grid( row=0, column=2 )         self.quitButton = Button ( self, text='Quit', command=master.destroy )            self.quitButton.grid( row=1, column=2 )     def selectFeat(self):         sel = self.stList.curselection()         myState = self.stList.get(sel[0])         arcpy.SelectLayerByAttribute_management(Layer, "NEW_SELECTION", "\"" + SelField + "\" = '" + myState + "'")         mxd = arcpy.mapping.MapDocument("CURRENT")         for df in arcpy.mapping.ListDataFrames(mxd):             if df.name == mxd.activeView:                 df.zoomToSelectedFeatures()                 arcpy.RefreshActiveView()         root = Tk() app = Application(master=root)                   app.master.title(Title) app.mainloop() ... View more

Run a gp.Describe(layer) and you can get out the Layer properties FIDSet. This is a clumsy string of all the FIDs that are in the selection. But it enables you to get a list of which features are selected. ... View more

ArcGIS 10 unfortunately does not install to be used as a development environment. Pythonwin is not installed. You need to install this IDE (or another of your choice). Before you do this I recommend that you fix the Python install. It can work properly if you understand that Esri installs another copy of Python in a non-standard place. This will result in other Python module installs not being seen because the paths are not common between the two versions. My solution is to keep the two items needed to point to ArcPy (Desktop10.pth) and the numpy folder in the non standard path c:/Python26/ArcGIS10.0/lib/site-packages, delete the whole Esri second copy, reinstall Python and Pythonwin  in the standard place (C:/Python26) and restore the Esri pointer file and numpy into C:/Python26/Lib/site-packages. Then Pythonwin will see the Esri modules and everyone else's modules, the Python help will work and things will be sweet. You set the editor to Pythonwin in ArcMap geoprocessing properties to allow you to right-click on the script tool to open in Pythonwin. After you install ArcGIS 10 SP2 you will have to repeat this exercise! It will add back an empty path C:/Python26/Desktop10.0/Lib/site-packages but there will be no second copy of Python so no geoprocessing tools will work. I can see a possible reason is to isolate the python version that is hard-coded into arcgisscripting or Arcpy from any other patches. However I have had no trouble with any 2.6.x release. The disadvantages outweigh the fear IMHO. ... View more

Yes, Python has good standard modules for random selections, permutations and other sampling techniques. You need to run a cursor first to get all the details of the parcels into a dictionary/list in memory so that you can accumulate random selections until you get to your target value. This will be very fast. Then create an SQL query with the selected parcel ids to process those parcels with a cursor. If you don't want to actually update the featureclass, even this would be unnecessary. The process to get the parcel details can run once outside your loop. Then when you run the process 100 times, it will make a selection, run the SQL query to make a layer and then process the layer. The SQL query would be in the form "parcel_id in (23,43,65,876,...)" Make sure that the parcel_id has an attribute index defined, it is not by default except for the objectid. ... View more

CalculateField is best kept for ModelBuilder. It simply wraps a cursor around the expression, so why not use a cursor if you are writing a script? It has a lot more capability: You can project data on the fly for calculation The field list can restrict the fields retrieved for efficiency You can use an sql expression to select a subset of records You can return records sorted It is possible to trap unexpected data It is much easier to program complicated functions All the other script variables are available It is just as fast as calculatefield When you add a spatial reference to the cursor the attributes are projected on the fly for returned variables such as the area calculation. There is no need to project your source. (You might also need a transformation to be set.) I personally do not hold any data in decimal degrees if I am going to run any geoprocessing tool or model. This is because I live in the temperate zone (37S - 47S) where the distortion between latitude and longitude scales is 50% and unacceptable for everything using distance, bearings, area or tolerances. [The only tool to properly account for geodetic calculations is the point buffer tool. Everything else assumes planar coordinates.] So what shape is your convex hull boundary? Does it have straight lines between points, or is it curved to trace the longitude paths? Before you reproject it, has it been densified to allow for non-linear curves? So I agree with Dan that it might be better to project the source points before you even run the convex hull tool. # ----- snip ----- projectedRef = gp.CreateObject("SpatialReference") projectedRef.CreateFromFile("MySr.prj") factor = 123223324 # the factor from map units to acres # in metric societies this is simply a power of 10 from metres to hectares cur = gp.UpdateCursor("checkhull","",projectedRef) row = cur.next() while row:     row.hullarea = row.shape.area * factor     cur.updateRow(row)     row = cur.next() del row,cur # must close to flush file buffers Very similar code for ArcPy at 10.x ... View more

My solution to merely updating the date used a Python module called elementTree, which was written some time ago and may be superceded by better tools, but it does the job at 9.3 well for me. http://docs.python.org/library/xml.etree.elementtree.html http://effbot.org/zone/element-index.htm I export the metadata for a given featureclass into XML. Unfortunately you cannot access the metadata in a geodatabase directly otherwise. There is no tool to unload metadata, you have to do this manually in ArcCatalog at 9.3, fixed at 10.0. If the featureclass was a shapefile it would already be an XML file. Then I load this template into a Python structure with elementTree that makes it easy to find each component for dates and change it. ElementTree then writes out the edited structure faithfully to a revised XML file. (While I was there I deleted the history) Finally there is a tool that will load metadata into a geodatabase. # 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 LINZ 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         pass     ## 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/customer/conservation/metadata") edDate = str(datetime.datetime.now().date()).replace("-","") edDate = publishDate # '20100914' gp.AddMessage(edDate+" edit date") gp.AddMessage(publishDate+" publish date") gp.AddMessage(createDate+" create date") ws = "e:/crs/customer/conservation/corax.gdb" metasrc = "e:/crs/customer/conservation/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"         gp.AddMessage(fc+" updated")     else :         print fcxml,"not found"         gp.AddError(fcxml+" 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 :         print tabxml,"not found"         gp.AddError(tabxml+" not found") # geodatabase metadata etxml = alter(metasrc+"/corax.xml",edDate,publishDate,createDate) gp.AddWarning("Load corax_et.xml to geodatabase by hand") There are a lot more tools at 10 to handle metadata. ... View more

There is a datatype in the tool properties that allows you to sketch a polygon (or point or line) graphic that is passed into the script. You can create a featureclass from it. The data type is a featureset. You need a template defined in a layer file to use it. You can sketch any type of polygon as in the draw tool, edit it and add attributes to fields defined in the template. ... View more

There is now a right-click specific tool to calculate centroids in the Calculate Geometry... into a field. You do not need to use the Field Calculator, but if you did, then you can use a Python expression instead of VBA which is now not supported. ... View more

If you are using Data Driven Pages in ArcGIS 10 you can set the dataframe rotation in a Python script before exporting each page. See the help example scripts. To get the rotation required I have in the distant past written an AML that found the diagonal corners of the page and calculated the average angle between the two points in two different projections to set the cogo bearing adjustment. You could do the same in Python and maybe store it in the page index for easy retrieval during printing. If the angle is not predictable because the parcel layout is the old flat earth projection, then maybe find a line in the view of each sheet (say a road centreline) and extract the angle from that to derive the angle required. You can find the bearing angle easily in Python using math.atan2(y,x) where these are the deltas between the line ends. Then subtract the bearing from 90 or 180 degrees to get the angle. ... View more

The Python code generated by Modelbuilder does not look like an analyst would code the problem. There is huge redundancy and repetition in it that is very inefficient scripting, but is Ok in ModelBuilder. (If you wrap a <CODE> script... </CODE> pair around your script (the # button on the toolbar) when posting it retains the indentation, otherwise we can't read it.) You say that you are using 9.3 but I see that the script is using ArcGIS 9.2 .create(), not .create(9.3). Things have moved on with Python now able to run inside the ArcGIS process at 9.3 and 10.0, avoiding creating a separate process for each tool. Recommendation: Upgrade to 9.3 or 10.0, and run 'in process'. Run it from a tool and check the box 'Run Python script in process' on the Source tab.You may not even need to refactor the code. But if you do... System toolboxes do not need to be added, only custom toolboxes not on the search path. Recommendation: Remove references to system toolboxes to save time and clutter. CalculateField is designed for ModelBuilder because cursors are not available. Cursors are just as fast and can have a lot more logic built in to do more than one thing on a pass. Recommendation: Replace multiple CalculateField steps with a single cursor and internal Python logic. Temporary featureclasses do not need to be written to disk, you can use "in_memory" workspace. Workspaces should always be a file geodatabase for speed and better indexing, and no size limits. Recommendation: Don't use dBase tables or shapefiles for intermediate featureclasses. MakeFeatureLayer does not need to spell out all the fieldmap if you are not modifying it. Recommendation: Leave the defaults empty if you are not changing it. Plain try/except blocks are not useful, they hide the error reporting. Remove them and replace them with specific error traps. Never have a plain except: block because it also traps keyboard interrupts. Specify what error you are trying to trap, eg arcgisscripting.ExecuteError, or Exception or something. Recommendation: Don't use try/except blocks in scripts unless you are actively trapping an error that you want to skip. In a script just let it fail, you will get enough feedback to identify the line and cause. Python errors are rare or even nonexistent since the script is entirely geoprocessing calls. Python has much faster structures such as lists, dictionaries and sets for finding duplicates. I find it is faster to load selected featureclass attributes once into these python structures, do the comparison and then write out the results, or create an SQL expression to make a layer and write out. ... View more

You would have to trap for a null record when there is a null value for cityname. while row :     # Process: Copy Features     try:               outFeature = os.path.join(env.workspace, row.CITY_NAME)         print outFeature         arcpy.CopyFeatures_management(row, outFeature)     except arcpy.ExecuteError:         print "failed to write"     row = searchCursor.next() [You do not have to add toolboxes for system tools. This comes from generating a script from Modelbuilder where this redundant and non-portable line is added.] Geoprocessing tools are designed to work on a whole featureclass or layer, not inside a cursor. Your script will only write one feature to a new featureclass, overwriting any previous featureclass. Is this what you intend? If you want more than one feature written out you would have to rewrite the process to select all the rows that meet an expression as a layer and then copy the layer to a new featureclass. You would need a list of citynames from either running a cursor to collect them or maybe Frequency and then open a cursor to get the list. There is no need to set the output path if you have set the workspace. arcpy.env.overwriteOutput = True arcpy.Frequency_management(cities_mexico,"in_memory/citList","CITY_NAME") for row in arcpy.SearchCursor("in_memory/citList"):   arcpy.MakeFeatureLayer(cities_mexico,"temp_lay","CITY_NAME = '"+row.CITY_NAME+"'")   arcpy.CopyFeatures("temp_lay",city) (I have not run this in a Python window so there may be syntax and case errors in this) ... View more

I didn't confuse the commands, I should have made it clearer that I meant any tool that 'merges rasters'. The Mosaic help does not warn you not to add a lot of tiles at once, and I expect that the mosiac dataset builder has similar limits, probably using the same underlying calls. It does not seem that the new command has solved all the space problems as we might have expected. Have you tried building the mosaic dataset a few tiles at a time? ... View more

Import win32print and set the printer properties. Open up the help for PythonWin and you will see a module win32print which has all the functions that appear in the Windows printer settings dialogs. However it would only set printer attributes, not custom driver attributes such as map size. import win32print win32print.GetDefaultPrinter() \\\\PUMA\\HPColourLaser ... View more