Latest Contributions by KimOllivier

Forget any idea of iterating over each row. Better analysis tools do that for you. Why not just use an overlay tool between the two featureclasses? If you want to find if a parcel centroid is inside some other polygon featureclass then create a point featureclass of parcel centroids and use Identity_analysis. Performance in my terms means finishing in preferably 30 seconds for something so simple, so you can run it several times until you get it right before morning tea. 😉 import arcpy import datetime start = datetime.datetime.now() arcpy.env.overwriteOutput = True arcpy.env.workspace = "e:/crs/customer/waiheke/mobile.gdb" arcpy.FeatureToPoint_management("parcel", "in_memory/parcel_center","INSIDE") arcpy.management.MakeFeatureLayer("E:/data/census2006/census2006.gdb/nztm/mb2006","jurs_lay","TA06V2 = '007'") arcpy.Identity_analysis("in_memory/parcel_center", "jurs_lay", "in_memory/parcel_jurs") # join back district to parcels arcpy.JoinField_management("parcel", "objectid", "in_memory/parcel_jurs", "ORIG_FID", ["MB06"]) print "Done",datetime.datetime.now() - start, " seconds" # >>> Done 0:00:26.953 seconds ... View more

Set your page size and map size to be the same. Pick the same ratio and a suitable size so that the number of dpi is easy to calculate. Set the percentage margin to zero. Leave on the grid for tests to check that you have got it right. I found it very difficult to get it to the exact pixel size. Unfortunately you cannot export the world file and projection if you are using layout view, only data frame view. If you need those, then use a grid with known dimensions so you can calculate the world file yourself and then use geotiffe.exe to populate the TIFF files afterwards. You may also be able to use PIL python module to do this. ... View more

The OBJECTID is a system primary key index of guaranteed unique values for that table that is read only, used by geoprocessing tools, scrolling and other system functions. The OBJECTID_12 name is the result of an operation that combined two tables with an existing non-primary OBJECTID. To avoid the primary key getting named in a nonstandard way, make sure that the source tables do not contain a 'normal' field called objectid before you do the analysis, otherwise it is too late to change it later. If you want to have a user-id that you can control or use as a 'foreign key' then you will have to name it something else. Esri have a policy of not allowing field renaming. Possibly because it breaks internal relationships such as indexes, relates and other complex objects, even if you have none. Some third party tools such as ET Geotools allow it, but it probably just does an add/copy/delete cycle which you could do also. You could try using an alias field name, but this cannot be used in geoprocessing, and it always catches me. Relational tables are specifically not allowed to have implicit order of fields or rows (see Codd's rules) so you have no control of stored order, only output in reports and displays. ... View more

Experiments with the sketch tool have shown up a few things that ArcGIS is doing behind the scenes. You can only use the polygon tool to create multipart polygon features. The box tool automatically finished the sketch prematurely for example. I could not find a way to add more parts after the sketch is finished. The 'parts' in the edit vertices window would be better called 'rings' because some will become parts, and others will be donut rings within a part. After finishing a sketch the rings will be converted into parts and donut rings, depending on whether they are enclosed by a ring or not. If the rings are intersecting, they are cracked and additional rings are created. It appears that the order is also recalculated to ensure clockwise. anticlockwise orders with nested rings. Donut vertices can be coincident at one point with their enclosing ring. (This is slightly different to the OGC rules for simple features) If you move a donut ring outside the outer ring it will turn into a separate part when saved, so you go from a single part to multiple parts. If there is a second ring inside the donut, then it turns into the donut, so we need a way of detecting the inner polygons. ... View more

Shock horror! A bug in my 9.3 script is true. My apologies, I will update the resource file as soon as I have investigated a bit more with some more examples. I am very puzzled that it is possible to have a validated feature and yet it does not appear valid. Surely this is not possible after all these years of polygon features? Maybe we are still missing some valid combinations in our scripts? For example it is valid to have ((a lake in an island) in a lake), and so on... recursively, which I have not tested yet. This would result in a second inner 'outer ring' inside an outer ring in the same part that is not a hole. The rings must also be ordered clockwise for an outer ring and anti-clockwise for a donut. So it would be useful to show the order of the polygons, maybe turn on arrows, to see if we are covering all possible cases. I can see that multi-part to single may destroy these apparent errors. It may be valid to have a preceding null point or two null points that would result in some very strange, but maybe "valid" polygons, without a donut, but an inner ring. ... View more

When you have a geoprocessing tool inside a cursor that loops for every feature you will have performance problems. That is because geoprocessing tools are designed to work most efficiently on a whole featureclass, not an individual feature. It should be a red flag, don't do it. The idea of ArcGIS is that it is geo-relational, for the same reason that SQL queries are 'non procedural'. You can't specify the processing for each tuple in table join so please think of a similar expression for GIS processing as well, instead of "reinventing GIS" in the scripting language. Either recast the problem to process the whole featureclass in one step or use a Python structure to index the attributes, such as a dictionary, which can hold many values in a list for each row. There are several alternative approaches that spring to my mind:     Use an arcpy.management. MakeQueryTable tool that allows a full SQL expression between tables with a one to many relation.     Run an arcpy.analysis.SummaryStatistics Tool (not Frequency which requires an ArcInfo licence) with a case field.     Run an arcpy.SearchCursor() once on a table to collect the counts and place in a python dictionary.     Then run an arcpy.UpdateCursor() once to populate the other table. Note that splitting the two cursors avoids multiplying the number of searches required. If you had 1 million records in each you would then have a billion slow cursor accesses. ... View more

You could use Linear Referencing. Build the polylines for routes, (CreateRoutes) and then add a table with a measure at the chosen interval. This will create a set of points along the route (and will not add one at the end for an odd remainder). You will have to create a dummy table with n x interval records that exceed the likely number. Does that disqualify the method because the easiest way to do that is run a python script to build a dummy table? But you only have to do it once. It can be reused. Trying to think of everything, if you want to add points at a lot of different polylines, each one has to be a separate route ID. def makeEvent(event):     ''' Creates a table of events     using fields plot and mark as a composite key plotmark     all fields hardcoded:         sampling_location         plotmark         plotid         mark         transect         distance_along_transect         fpc_number     '''     try:         curIn = arcpy.InsertCursor(event)         for row in arcpy.SearchCursor(routeFC):             plotmark = row.plotmark             plotid = row.plotid             mark = row.mark             for pt in range(1,traps+1): # one-based counting                 rownew = curIn.newRow()                 rownew.sampling_location = plotid # field Increased                 rownew.transect = mark.strip()                 rownew.trap_number = int(pt)                 rownew.distance_along_transect = int(pt*spacing) # long??                 rownew.plotmark = plotmark                 curIn.insertRow(rownew)         del curIn         del rownew         del row         countrow = str(arcpy.GetCount_management(event))         arcpy.AddMessage(event+" "+countrow)         return None ... View more

Perhaps you can look at your process. Running "network in a loop" rings alarm bells. Perhaps you need to remove previous runs in the loop to release memory? If the first pass works, then you must clean up. Network results are held in memory, perhaps write them out to a file if you need to keep them. ... View more

Create a date field in your geodatabase (not shapefile) and populate it using the datetime object in Python. The datetime module handles all the messy calculations very well and the animation tools need a date field. If you have the source data as  separate strings there are easy conversions to a datetime object using datetime.combine(date,time) The datetime object also knows about time intervals, so you can subtract two dates to get the differences easily, or add useful offsets such as a week, month or year without worrying about days in the month or leap years. eg import datetime print datetime.datetime.now() ... View more

At 10.0 there is a tool to turn a table into points and another to turn points into lines. To save some coding. ... View more

A text file can be loaded as a table, you would need to have a valid header and it helps a lot if you also have a schema.ini file with an entry to define field types, otherwise they all default to char fields 255 wide. ... View more

Have a look at the win32 interface 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) Tkinter provides a more complete dialog designer. ... View more

It's much harder than that! I frequently want to adjust stripmap angles after they have been created and it is very difficult to calculate the angle of a polygon. COGO tools are no use because they only handle lines. You would have to split the polygons back to lines and then you would still not have a useful angle. I wrote my own Python script and used math.tan2(y/x) to get bearings, that then had to be normalised. You could use the new tool in 10.x Data Management>Features>Maximum Bounding Geometry and set the option "Add geometry" which includes a nice field MBG_orientation suitable to use, but of course it is now in a new featureclass with all your nice attributes trashed, so you would have to make sure you group by SMAP_NUM , join back and recalculate. Still a bit less work that my Python script. But here it is anyway in case it helps you. # UpdateAngle.py # assumes rectangular polygons # find bearing of longest side # might have trouble with up # 2 June 2011 # Kim Ollivier import arcgisscripting from math import * import datetime gp = arcgisscripting.create(9.3) try:     layer = sys.argv[1] except:     layer = "stripindex"     ws = 'e:/teararoa/nz/current.gdb'     gp.workspace = ws if not gp.Exists(ws):     print "ws",ws     raise Exception rotField = "smap_angle" if  not gp.ListFields(layer,rotField):     gp.AddField(layer,rotField,"FLOAT") def calcbearing(dPt):     '''     calculate bearing of longest side     '''     # points     pt1 = dPt[1]     pt2 = dPt[2]     pt3 = dPt[3]     # lengths     l1 = hypot(pt2[0] - pt1[0],pt2[1] - pt1[1])     l2 = hypot(pt3[0] - pt2[0],pt3[1] - pt2[1])     # bearings     b1 = atan2(pt2[1] - pt1[1],pt2[0] - pt1[0])     b2 = atan2(pt3[1] - pt2[1],pt3[0] - pt2[0])     # print "test",round(l1),round(l2),degrees(b1),degrees(b2)     # print dPt     if l1 > l2 :         degree = degrees(b1)     else:         degree = degrees(b2)     bearing = 90.0 - degree     if  bearing > 90:         bearing = bearing - 180     if  bearing < -90:         bearing = bearing + 180     return bearing cur = gp.UpdateCursor(layer) ##,"ID = 1") row = cur.next() print while row:     feat = row.Shape     for partNum in range(feat.partCount) :         part = feat.getPart(partNum)         dPt = {}         n = 0         for pt in iter(lambda: part.next(),None): # stops at null pt             n+=1             dPt = (pt.X,pt.Y)     rotation = calcbearing(dPt)     row.setValue(rotField,rotation)     row.setValue("LastUpdate",datetime.datetime.now())     # print row.objectid,row.id,row.island,rotation #,dPt     cur.updateRow(row)     row = cur.next() del row,cur Next the challenge will be to write something that will rotate the feature from the attribute. Numpy will be useful here. ... View more

Use a simpler SQL expression like sql expr = "SP10 in ('WS','BS','AP') AND SP11 in ('MM','GA')" If your logic is hard to follow, make subsets and select multiple times using the 'SELECT_FROM...' options. Check each sql query produces the desired values by running it interactively in ArcMap table selections. Also get rid of all the extra double quotes. If you use valid field names you don't ever need to quote them. (ie start with a character, no spaces or punctuation characters, no reserved words) [You do need square brackets for MS Access, but that doesn't count] ... View more

Finally we have got a good description of the problem! The method  only needs 2 vertices from the millions. That simplifies creating the path from the sun angles considerably. I only used the points to find the outer pairs and generate the swept area. I then unioned the original polygon and a moved copy. There would be better ways of finding them instead of by inspection as I did for the proof of concept. That is a red-herring in the end, no need for points. To find the max and min sides there are a number of coordinate geometry solutions. One easy non-algebraic way would be to rotate the polygon(s) by the final angle and find the bounding box to find the two points/vertices, or just the min and max y. These are the only two points that need to be extended by the path to create the swept area. Creating the paths is just a matter of using the table of azimith, dx values to accumulate x,y pairs. This whole solution works because there is no rotation of the polygon. But if you want to calculate a sun shadow from 3D features there is already a tool in 3D Analyst "Feature Shadow Volume" that will take a range of times, for a whole year in necessary. ... View more