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

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