8

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.

1

4 Answers 4

10

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

7
  • 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.
    – Mehdi
    May 17, 2016 at 4:58
  • What do you call this raster in mxd table of content?
    – FelixIP
    May 17, 2016 at 5:37
  • I tried DEMH and DEMH_copy as raster names.
    – Mehdi
    May 17, 2016 at 11:22
  • Could you see field demh in shapefile table after crash?
    – FelixIP
    May 17, 2016 at 18:20
  • the field "demh_copy_tif" is added to the attribute table but with NULL values.
    – Mehdi
    May 18, 2016 at 0:03
6

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.

0
4

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

1
  • 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. Mar 30, 2016 at 22:46
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.