# Extract Raster Values using Numpy

7134
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 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 arcpydef 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 secondsdef 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 pntsdef 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, coldef 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) / 2def getExtentFromPolyline(line):    return line.extentdef getExtentFromFC(fc):    return arcpy.Describe(fc).extentdef 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.

Comments 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
Revision #:
1 of 1
Last update:
‎12-23-2014 04:17 PM
Updated by: Contributors