3

I'm relatively new to arcpy and python.

I have a script that creates triangular polygons given two point shapefiles of centroid points and vertexes around them. I need to run this over ~170k points, which would take more than a week at my available processing power, especially if Murphy's Law strikes during the script execution. I haven't had much luck digging through the Arcpy docs and coming up with anything that could give me the same results in less time.

Can anyone point me in the direction of a different tool chain for this if they can think of one?

Or notice anything in the code that would hamper performance?

with centroidCursor as cecursor:
    for record in cecursor:
        whereClause = "IN_FID = " + str(record[0])
        selection = arcpy.SelectLayerByAttribute_management(nearTableView, "NEW_SELECTION", whereClause) # get the all the records that have that centroid ID from the near table

    with arcpy.da.SearchCursor(selection, "NEAR_FID") as selcursor: # then append the vertex IDs to a list
        tempList = [row[0]] # create temporary list
        for selrow in selcursor:
            tempList.append(selrow[0]) # iterate though rows appending vertex IDs to list

        idNum = record[0]
        whereClause = "OBJECTID = {0} OR OBJECTID = {1} OR OBJECTID = {2}".format(tempList[0],tempList[1],tempList[2])
        selected = arcpy.SelectLayerByAttribute_management(cornerpoints, "NEW_SELECTION", whereClause) # select the vertexes from the list
        lines = arcpy.PointsToLine_management(selected, "in_memory", Close_Line="CLOSE") # convert points to a temporary lines
        poly = arcpy.FeatureToPolygon_management(lines, workspace + "/mesh_polys/polygons{0}".format(idNum)) # then convert lines to polygons stored in a dataset

I'm going to first try going straight from points to polygons with the Geometry class. From the discussion in the comments the follow up will be one of two methods:

A - Using a layer definition query inside the nested searchCursors (I'm unfamiliar with this and it's throwing a bunch of different errors as I'm working with it after FelixIP pointed me in it's direction).

B - creating a custom table with some kind of relevant information tied to the cursor_ID (vertex geometry, perhaps?).

I now understand what the definition queries do and they are absolutely faster.

 whereClause = "OBJECTID IN ({0} , {1} , {2})".format(tempList[0],tempList[1],tempList[2])
 cornerpoints.definitionQuery = whereClause
        selection = cornerpoints
        lines = arcpy.PointsToLine_management(selection, workspace + "/mesh_lines/line{0}".format(row[0]), Close_Line="CLOSE")# convert points to a temporary lines

However, I'm now getting this error when I try to run the script:

ExecuteError: IOError: "GPL0" does not exist

There's very little about what can cause this error that I can find. Anyone know what I'm missing?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
4xle
  • 123
  • 6
  • You may be restricted by available license level, but if not I'm assuming interpolating your vertex points to a TIN or creating Thiessen polygons around your centroids isn't an option? TIN would definitely be triangles, but Thiessen wouldn't necessarily be (probably wouldn't). Another thing that comes to mind is spatially joining the vertices to the nearest three points to get a polygon ID attribute and then using that to restrict vertices used in polygon creation. – Chris W Apr 18 '15 at 22:04
  • Tried both already. TIN got really close to creating the right network, but I don't believe there is a way to make it respect the proximity of the vertexes to the centroids. And the resulting Thiessen polygons were triangles with oddly rounded corners that didn't align with the respective vertexes, but were otherwise generally in the right places. I can understand using a spatial join to effectively get a polygon ID (I'm using the centroid ID to label it now) - are you referring to the Polygon class in arcpy as the creation method? – 4xle Apr 18 '15 at 22:45
  • Ah, so the polygons might have gaps between? TIN wouldn't pay attention to centroids at all, it would just look at the vertices. I'm not an experienced scripter (understatement) and what I suggest may be pretty much what you are already doing with the Near table except the spatial join would theoretically generate three points for every current vertex (no idea on performance or time). I wasn't seeing/paying attention to any current connection between vertices and centroids in attribute form. I think what I'm suggesting eliminates the search cursors and just feeds attributed points to PtL. – Chris W Apr 18 '15 at 23:03
  • Replace select by attributes to layer definition query. B) replace object id = by objectid in () and make it definition query. – FelixIP Apr 19 '15 at 04:35
  • Thanks Felix, I will give LDQ a shot. I've never used it before so I'll have to check the docs on it's implementation. – 4xle Apr 19 '15 at 06:52

1 Answers1

3

Having looked at your code where is row[] being created, I guess we are not seeing the full picture?

I've seen this "GPL0" error before. It's indicating that the first parameter is invalid. If you look at the syntax section of the help for this tool it states the input features are of Data Type FeatureLayer. As we are not seeing the full code we have no idea how you are creating cornerpoints, is it a FeatureLayer?

I know this is easier said than done, but if it were me I would have abandon the approach you have taken and built the geometries manually. You want to avoid running intensive processing where you are constantly calling a tool within a loop. As I've understood it there is always a cost when calling a tool. Have a look at this page it shows how to build geometries using arcpy.

The first part of your code that builds tempList could be jigged around to read once into a dictionary which would take about 2-3 seconds, then use that to get the points that build the triangles. Again you might want to consider reading those into a dictionary first and then its just reading from dictionaries which will be blisteringly fast!

Hornbydd
  • 43,380
  • 5
  • 41
  • 81
  • +1 for dictionary, where key=[in_fid], values=list of near_fids. I'd add +1 for geometries. Also replace iteration through 1st table NearTableView by TableToNumpyArrayy – FelixIP Apr 19 '15 at 20:40
  • I knew that there had to be a better way than what I was doing, a dictionary is so simple that I forgot that it could fit what I would need perfectly. Assuming I don't run into any weirdness with my files (which I don't think I will), this will probably solve all my problems. Now I'll just have to build it. – 4xle Apr 19 '15 at 22:14
  • Don't forget to "solved" next to @Hornbydd answer please – FelixIP Apr 20 '15 at 00:07
  • Of course! And thank you for your input as well Felix, it's much appreciated. For the record, 'blisteringly fast' is the way to describe the results, it takes about three minutes to generate the polygons. Cheers! – 4xle Apr 20 '15 at 22:08
  • So you've turned a week of crunching into just 3 minutes... Awesome! Definitely worth treating yourself to a pint! :) – Hornbydd Apr 20 '15 at 22:59