1

enter image description hereI have specific points in shp file and list of rasters in the same coord system and want to extract values from rasters for each point and place these values in column in Attribute Table of shp file of points.

Rasters are separate files.

Is it possible to go from vector to raster?

Shp file pf points:

infc = 'samplepoints1.shp'

code just a start:

import arcpy

infile = r"C:\Users\anna_pavlenko\Desktop\SILVER\biscuit_dnbr.tif"
#create raster object
ftc = arcpy.sa.Raster(infile)
#rows and columns
ftc.height 
ftc.width 
#cell size
cellSize = ftc.meanCellHeight
print 'The cell  size is ', cellSize, 'm'

rasterList = r"C:\Users\biscuit_dnbr.tif"
rstArray = arcpy.RasterToNumPyArray(rasterList)

rows, cols = rstArray.shape
print "number of rows: " ,  rows
print "number of columns: ", cols

for rowNum in xrange(rows):
    for colNum in xrange(cols):
        value = rstArray.item(rowNum, colNum)

If I use

inRasterList = [["biscuit_dnbr.tif", "foothills_dnbr.tif", "dNBR1"] ,     ["biscuit_rdnbr.tif", "foothills_rdnbr.tif", "RdNBR1"],["silver_dnbr.tif", "pony_complex_dnbr.tif", "dNBR2"], ["silver_rdnbr.tif", "pony_complex_rdnbr.tif", "RdNBR2"]] 
#  ["foothills_dnbr.tif","dNBR1"], "foothills_rdnbr.tif", "RdNBR1"],
#  ["pony_complex_dnbr.tif", "dNBR2"], "pony_complex_rdnbr.tif", "RdNBR2"]]

ExtractMultiValuesToPoints(inPointFeatures, inRasterList, "NONE")

This function create values for 1 place (polygon) and extract value for point inside this polygon only and place them in columns: "dNBR1", "dNBR2", "RdNBR1", "RdNBR2".
For another set of rasters it create 0 -zeros for first points and add new columns : "dNBR1_1", "dNBR2_1", "RdNBR1_1", "RdNBR2_" each time.

1
  • There is a tool for this task sample sa
    – FelixIP
    Commented Apr 18, 2016 at 5:53

3 Answers 3

3

If you don't have the spatial analyst extension or want to use Numpy, you can convert your x and y point coordinates to the row and column that is close to that point. The Extract Values tools have many more options, and produce slightly different results than those in the code below. I'm not advocating one method over another. I also don't know if this would be any faster.

import numpy as np

Get information about your raster file:

rast = arcpy.Raster(r"somePath")
desc = arcpy.Describe(rast)

Get the upper left corner of the raster and spatial reference system (you can ignore the spatial reference if all your points and rasters are the same).

ulx = desc.Extent.XMin
uly = desc.Extent.YMax
#Spatialrefernce
sr = desc.spatialReference

Convert the raster to a numpy array

rstArray  = arcpy.RasterToNumPyArray(rast)

Access your Feature Class using an update cursor to add the value. I'm using the projectAs tool to make sure the point is in the same coordinate system as the raster. Calculate the row and col based on how far the point is from the upper left corner. This can very depending on rounding. Then pass this into your numpy array.

with arcpy.da.UpdateCursor("NHDPoint",["SHAPE@","value"]) as uc:
    for row in uc:
        pnt = row[0].projectAs(sr)
        #assuming every point falls to the left and below uperleft corner
        deltaX = pnt.centroid.X - ulx
        deltaY = lly- pnt.centroid.Y
        arow = int(deltaY/rast.meanCellHeight)
        acol = int(deltaX/rast.meanCellWidth)
        row[1] = rstArray[arow,acol]
        uc.updateRow(row)

This code assumes that the points will fall to the right of upper left corner, and below it. It also just used the default rounding when converting the float to an integer.

Another option, which might be closer to what the Extract tools is to use a distance to the center of the cell from your points to determine which value to use.

1
  • I was able to get this to work by adding lly = desc.Extent.YMin as it was failing on deltaY = lly- pnt.centroid.Y
    – tpcolson
    Commented Nov 18, 2018 at 16:09
2

In addition to the tool pointed out by HDunn, you can use Extract Multi Values to Points if you have many raster data with the same point shapefile. Here is another snippet from ArcGIS help.

import arcpy
from arcpy.sa import *
from arcpy import env 
env.workspace = "c:/sapyexamples/data"
ExtractMultiValuesToPoints("observers.shp", [["elevation", "ELEV"], 
                           ["costraster", "COST"], ["flowdir", "DIR"]], "NONE")
2
  • This function create values for
    – Adams
    Commented Apr 18, 2016 at 6:35
  • The above function will update the original file and add fields for of raster names used, which means this tool modifies the input data. Make sure to have a copy of your original input shapefile as backup.
    – ahmadhanb
    Commented Apr 18, 2016 at 6:48
1

If you're already using ArcPy, You might as well use the Extract Values to Points tool (assuming you have Spatial Analyst extension)

Here's the code snippet from ESRI's site

ExtractValuesToPoints("rec_sites.shp", "elevation",
                  "C:/sapyexamples/output/outValPnts","INTERPOLATE",
                  "VALUE_ONLY")

Iterate over your raster layers like you did and extract values to your point layer

1
  • This function place output value in separate file - the values in the attribute table of an output feature class. I want in the same file of point infc = 'samplepoints1.shp' just to add new column in Attribute Table.
    – Adams
    Commented Apr 18, 2016 at 6:31

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.