Latest Contributions by XanderBakker

I was looking for a way of extracting elevation values from a DEM and decided to use numpy for this purpose. The idea was to use a polyline, extract points every x m (interval is the raster cell size) and obtain the elevation for each point. At first I tried using the “GetCellValue_management” tool, but when you do that a 100.000 times, it is pretty slow. When you convert a raster to numpy, the numpy array does not know the coordinates of each element of the array. To solve that you have to establish that relationship. What I did was: get the extent of the polyline adapt the polyline extent using the raster extent and cell size create a numpy array using the adapted extent (don’t extract more info than you need) create a list of points from the polyline and using the cell size as interval for each point determine the row and column of the point coordinates extract the elevation value using the row column information. This method turned out to be over 170x faster than using the repetitive use of the tool “GetCellValue_management” and still does not require a license level higher than ArcGIS for Desktop Basic and there is no need for the Spatial Analyst (or any other extension) either. #------------------------------------------------------------------------------- # Name:        dem_numpy.py # Purpose:     extract dem values along a given line using numpy # # Author:      xbakker # # Created:     23/12/2014 #------------------------------------------------------------------------------- import arcpy def main():     import datetime     start_time = datetime.datetime.now()     # inputs     dem = r"D:\Xander\Numpy\fgdb\test.gdb\DEM"     fc = r"D:\Xander\Numpy\fgdb\test.gdb\Line"     dist_from = 25000     dist_to = 125000     # get polyline and extract a part     line = getFirstPolyline(fc)     subline = line.segmentAlongLine(dist_from, dist_to, False) # 10.3 required     cellsize = getRasterCellSize(dem)     # get extent from line and adapt extent to raster     ext1 = getExtentFromPolyline(subline)     ext2 = adaptExtentUsingRaster(ext1, dem, cellsize)     # create numpy array     np = createNumpyArrayUsingExtent(dem, ext2, cellsize)     # loop through points and extract values from raster     d = 0     pnts = getListOfPointsFromPolyline(subline, cellsize)     for pnt in pnts:         row, col = getRowColExtentFromXY(ext2, pnt, cellsize)         z = np.item(row,col)         print "{0}\t{1}\t{2}\t{3}\t{4}\t{5}".format(d, col, row, pnt.X, pnt.Y, z)         d += cellsize     end_time = datetime.datetime.now()     delta = end_time - start_time     print "Elapsed: {0} seconds".format(delta.total_seconds())     # Elapsed: 4.519 seconds def getListOfPointsFromPolyline(line, interval):     pnts = []     d = 0     max_len = line.length     while d < max_len:         pnt = line.positionAlongLine(d, False)         pnts.append(pnt.firstPoint)         d += interval     pnts.append(line.lastPoint)     return pnts def getRowColExtentFromXY(ext, pnt, cellsize):     # r, c start at 0 from upper left     col = int(((pnt.X - ext.XMin) - ((pnt.X - ext.XMin) % cellsize)) / cellsize)     row = int(((ext.YMax - pnt.Y) - ((ext.YMax - pnt.Y) % cellsize)) / cellsize)     return row, col def createNumpyArrayUsingExtent(raster, ext, cellsize):     lowerLeft = arcpy.Point(ext.XMin, ext.YMin)     ncols = int(ext.width / cellsize)     nrows = int(ext.height / cellsize)     return arcpy.RasterToNumPyArray(raster, lowerLeft, ncols, nrows) def adaptExtentUsingRaster(ext, raster, cellsize):     ras_ext = arcpy.Describe(raster).extent     xmin = ext.XMin - ((ext.XMin - ras_ext.XMin) % cellsize)     ymin = ext.YMin - ((ext.YMin - ras_ext.YMin) % cellsize)     xmax = ext.XMax + ((ras_ext.XMax - ext.XMax) % cellsize)     ymax = ext.YMax + ((ras_ext.YMax - ext.YMax) % cellsize)     return arcpy.Extent(xmin, ymin, xmax, ymax) def getRasterCellSize(raster):     desc = arcpy.Describe(raster)     return (desc.meanCellHeight + desc.meanCellWidth) / 2 def getExtentFromPolyline(line):     return line.extent def getExtentFromFC(fc):     return arcpy.Describe(fc).extent def getFirstPolyline(fc):     # return first feature geometry     return arcpy.da.SearchCursor(fc, ("SHAPE@",)).next()[0] if __name__ == '__main__':     main() Please note that the code assumes that the polyline is completely within the raster extent and there is no handling of no data values in the raster. ... View more

Last month a user mentioned the Dijkstra algorithm to determine the shortest path between two nodes in a network. My first reaction was to use Network Analyst instead, but to honor the essence of a previous useless post: Japanese visual multiplication with lines using arcpy and a dash of patriotism (Dijkstra was Dutch) made me look into this a little more. I came across this post; Dijkstra's algorithm for shortest paths « Python recipes « ActiveState Code  and from this post I actually used the code from this one; Rebrained! » Blog Archive » Shortest Path – Python . Remember, basic skills for developing code are search, copy and paste… It is basically recursive code that will ultimately get you the route between two nodes on the network. What you need is a "graph". This is a nested dictionary that lists all the connections between the from and to nodes and the corresponding distances. To be able to use the code I need to convert my "network" into this graph format. What you need is a network (a simple line featureclass), in which all intersections of lines consist of nodes. You will also need a point feature class containing the start and end point. These points should correspond with nodes of the network.This is the result: I came up with the following code: #------------------------------------------------------------------------------- # Name:        dijkstra.py # Purpose:     determine route # # Author:      xbakker # # Created:     02/12/2014 #------------------------------------------------------------------------------- import arcpy import sys def main():     # input featureclasses     fc_red = r"D:\Xander\GeoNet\Dijkstra\fgdb\test.gdb\red"     fc_stops = r"D:\Xander\GeoNet\Dijkstra\fgdb\test.gdb\stops2"     # create line dictionary     dct_lines = createLineDict(fc_red)     # create node and inverted node dictionary     dct_nodes = createNodes(dct_lines)     # update line dictionary with nodes     dct_lines_nodes = combineLinesAndNodes(dct_lines, dct_nodes)     # create graph     dct_graph = createGraph(dct_nodes, dct_lines_nodes)     # determine the node of start and end of route     node_start, node_end = getStartEndNodes(fc_stops, dct_nodes)     # determine route     result = shortestpath(dct_graph, node_start, node_end)     # create line with the result     dist = result[0]     route = result[1]     print "Resulting route    : {0}".format(route)     print "Length of the route: {0}\n".format(dist)     cnt = 0     for node_id in route:         cnt += 1         print "{0}\t{1}\t{2}".format(cnt, dct_nodes[node_id][0], dct_nodes[node_id][1]) def createLineDict(fc):     """Crear dict de lineas con datos relevantes"""     #OBJECTID (line): length, (Xstart, Ystart), (Xend, Yend)     dct_lines = {}     flds = ("OID@", "SHAPE@", "SHAPE@LENGTH")     dct_lines = {r[0]: [r[2], (r[1].firstPoint.X, r[1].firstPoint.Y),                 (r[1].lastPoint.X, r[1].lastPoint.Y)]                 for r in arcpy.da.SearchCursor(fc, flds)}     return dct_lines def createNodes(dct):     """Crear dict de nodos"""     node_format = "node_{0}"     cnt = 0     dct_nodes = {}     dct_nodes_inv = {}     for line in dct.values():         start = line[1]         if not start in dct_nodes.values():             cnt += 1             node_id = node_format.format(cnt)             dct_nodes[node_id] = start         end = line[2]         if not end in dct_nodes.values():             cnt += 1             node_id = node_format.format(cnt)             dct_nodes[node_id] = end     return dct_nodes def combineLinesAndNodes(dct_lines, dct_nodes):     """Crear dict con lineas y codigo de nodos"""     dct_lines_nodes = {}     for oid, line in dct_lines.items():         start = getNodeID(line[1], dct_nodes)         end = getNodeID(line[2], dct_nodes)         dct_lines_nodes[oid] = [line[0], start, end]     return dct_lines_nodes def getNodeID(pnt, dct_nodes):     for node_id, coords in dct_nodes.items():         if coords == pnt:             return node_id             break def createGraph(dct_nodes, dct_lines_nodes):     res = {}     for node_id in dct_nodes.keys():         links = {}         nodes_start = [[line[2], line[0]] for line in dct_lines_nodes.values() if line[1] == node_id]         for node in nodes_start:             links[node[0]] = node[1]         nodes_end = [[line[1], line[0]] for line in dct_lines_nodes.values() if line[2] == node_id]         for node in nodes_end:             links[node[0]] = node[1]         res[node_id] = links     return res def getStartEndNodes(fc, dct_nodes):     cnt = 0     with arcpy.da.SearchCursor(fc, ("SHAPE@")) as curs:         for row in curs:             cnt += 1             if cnt == 1:                 pnt = row[0].firstPoint                 node_start = (pnt.X, pnt.Y)                 node_start = getNodeID(node_start, dct_nodes)             elif cnt == 2:                 pnt = row[0].lastPoint                 node_end = (pnt.X, pnt.Y)                 node_end = getNodeID(node_end, dct_nodes)             else:                 break     return node_start, node_end def shortestpath(graph,start,end,visited=[],distances={},predecessors={}):     """Find the shortest path between start and end nodes in a graph"""     # we've found our end node, now find the path to it, and return     if start==end:         path=[]         while end != None:             path.append(end)             end=predecessors.get(end,None)         return distances[start], path[::-1]     # detect if it's the first time through, set current distance to zero     if not visited:         distances[start]=0     # process neighbors as per algorithm, keep track of predecessors     for neighbor in graph[start]:         if neighbor not in visited:             neighbordist = distances.get(neighbor,sys.maxint)             tentativedist = distances[start] + graph[start][neighbor]             if tentativedist < neighbordist:                 distances[neighbor] = tentativedist                 predecessors[neighbor]=start     # neighbors processed, now mark the current node as visited     visited.append(start)     # finds the closest unvisited node to the start     unvisiteds = dict((k, distances.get(k,sys.maxint)) for k in graph if k not in visited)     closestnode = min(unvisiteds, key=unvisiteds.get)     # now we can take the closest node and recurse, making it current     return shortestpath(graph,closestnode,end,visited,distances,predecessors) if __name__ == '__main__':     main() Have fun! ... View more