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.