Extracting Values to Points without ArcGIS Spatial Analyst?

Question

I have a raster DEM & multiple points. I want to determine the DEM elevation at each point. Other than manually examining the DEM at each point, (very time consuming with lots of points) is there any way to do this without access to the Spatial Analyst Extract Values to Points tool? I do have access to 3d Analyst if that helps.

• Aaron: This may be helpful: gis.stackexchange.com/questions/46893/…

Answer 1

Script:

import arcpy, traceback, sys, os
pntFile=arcpy.GetParameterAsText(0)
rasters=arcpy.GetParameterAsText(1)
rasters=rasters.split(';')
theFields=[x.name.lower() for x in arcpy.ListFields(pntFile)]
result=arcpy.GetCount_management(pntFile)
nF=int(result.getOutput(0))
p=arcpy.Point()
try:
 def showPyMessage():
 arcpy.AddMessage(str(time.ctime()) + " - " + message)
 for raster in rasters:
 desc=arcpy.Describe(raster)
 theFLD=raster.lower()
 arcpy.AddMessage("Sampling "+theFLD)
 if not(theFLD in theFields):
 try:arcpy.AddField_management(pntFile, theFLD, "FLOAT")
 except:
 aF=raster.split(os.sep)[-1]
 theFLD=aF.split(".")[0]
 theFLD="F"+theFLD
 arcpy.AddField_management(pntFile, theFLD, "FLOAT")
 arcpy.SetProgressor("step", "", 0, nF)
 with arcpy.da.UpdateCursor(pntFile,("SHAPE@XY",theFLD)) as rows:
 for row in rows:
 p.X,p.Y=row[0]
 myArray = arcpy.RasterToNumPyArray(raster,p,1,1,-9999)
 row[1]=myArray[0,0]
 rows.updateRow(row)
 arcpy.SetProgressorPosition()
 del row,rows
except:
 message = "\n*** PYTHON ERRORS *** "; showPyMessage()
 message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
 message = "Python Error Info: " + str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()

Interface:

enter image description here

Parameters:

enter image description here

to work from mxd. Give raster(s) short unique name in table of content. Script will create same name field and attempt to populate it with cell values.

It is slow, thus I use it with less than 1000 points

• Mehdi: Hi, thanks for the script. I tried to run it for a couple of times but it gives this error: "line 27, in <module> for row in rows: Tue May 17 12:56:55 2016 - Python Error Info: <type 'exceptions.RuntimeError'>: A column was specified that does not exist." any ideas? thanks.
• FelixIP: What do you call this raster in mxd table of content?
• Mehdi: I tried DEMH and DEMH_copy as raster names.
• FelixIP: Could you see field demh in shapefile table after crash?
• Mehdi: the field "demh_copy_tif" is added to the attribute table but with NULL values.

Answer 2

3D analyst's Add Surface Information will add a Z field to your vector data with the data value from the overlapping raster layer:

Interpolates surface elevation properties for point, multipoint, and polyline features.

That's for v10, I didn't catch which version of ArcGIS you were using.

Answer 3

You can do this with Hawthorne Beyer's free Geospatial Modelling Environment, (GME, formerly known as Hawth's Tools). There is a tool in there, Intersect Points With Raster, which as its name implies, acts like the Intersect tool in ArcGIS but allows you to intersect a point layer with a raster, like the Extract Values to Points tool. You can also apply an SQL query to your point layer to only extract raster values for a subset of points. You will need to install or upgrade R on your PC, and GME as well (it relies on ArcGIS but it can run independently of an ArcGIS session).

• Jeffrey Evans: At this point you may as well use R. Including adding libraries and reading data, it only takes about 4 lines of code to do this and you can query the data on the fly using an index syntax.

Answer 4

Another option to extract raster values to point shapefile is to use QGIS which is a free open source software. From plugin manager in QGIS download Plugin: Point sampling tool, which can extract multiple raster values to point shapefile similar to Extract Multi Values to Points which requires Spatial Analyst extension in ArcGIS.

The tool works perfectly if the projection of the raster(s) data and the shapefile are same. Working with different projections for the raster(s) and shapefile data will create a shapefile with empty columns of raster values.