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.
Good job...We should start Numpy Snippets ...
