# Extract Raster Values using Numpy

7716
5
12-23-2014 04:17 PM

## Extract Raster Values using Numpy

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 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)

# 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)

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()

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. Nice Xander  !!! fix def spelling line 32 and 66....Arrat... to ...Array..

There is lots you can do with geometry objects and numpy without licenses ie your defs beginning lines 49,60 and 72 Good job...We should start Numpy Snippets ... Lead the way, and I'll follow! ... or do we first need to get some python editor with numpy and arcpy support on the iPad? edit: thanks for the catch! I will get back on the PC tonight...during the interim, I do have a python implementation for the iThingy but no external modules are available  Same here on my Android devices... would be cool though... ID = [ 0,1,2,3,4,5,6,7,8,9]

Ys = [9,7,6,3,4,1,5,8,2,0]

Xs =[0, 2, 1, 3, 5, 4, 7, 6, 9, 8]

XY_array = np.array(zip(ID,Xs,Ys), dtype= [('ID','<i4'),('X','<f8'),('Y','<f8')])

print XY_array

array([(0, 0.0, 9.0), (1, 2.0, 7.0), (2, 1.0, 6.0), (3, 3.0, 3.0),

(4, 5.0, 4.0), (5, 4.0, 1.0), (6, 7.0, 5.0), (7, 6.0, 8.0),

(8, 9.0, 2.0), (9, 8.0, 0.0)],

dtype=[('ID', '<i4'), ('X', '<f8'), ('Y', '<f8')])

Just bought Pythonista from the iStore for \$6.99, comes with numpy, matplotlib and a few others I haven't looked at yet ( Pythonista  )

No syntax highlighting on the iThingy however....baby steps I suppose

Version history
Last update:
‎12-12-2021 03:37 AM
Updated by: 