## Thursday, October 11, 2012

### Extract Raster Values from Points

The R blog article encourages me to write this solution to extract Raster values from points in R.

In geospatial analysis, extracting the raster value of a point is a common task. If you have few raster files or few points; you can extract the raster value by overlaying a point on the top of the raster using ArcGIS. What will you do, if you have hundreds of raster files and thousands of points? The easy solution is use loop in Python and ArcGIS. Is loop efficient to use? No. Can loop be avoided? Yes.

Then how?

Step 1: Create a Raster stack or Raster brick of your raster files using “raster” package in R.
For example:
rasStack = stack(raster1, raster2,raster3 …rasterN)

Step 2:  Read point data, and convert them into spatial points data frame.
Sample: pointfile.csv
 Point_ID LONGITUDE LATITUDE 1 48.765863 -94.745194 2 48.820552 -122.485318 3 48.233939 -107.857685 4 48.705964 -99.817363
For example:
coordinates(pointCoordinates)= ~ LONGITUDE+ LATITUDE

Step 3: Extract raster value by points
For example:
rasValue=extract(rasStack, pointCoordinates)

Step 4:  Combine raster values with point and save as a CSV file.
For example:
combinePointValue=cbind(pointCoordinates,rasValue)
write.table(combinePointValue,file=“combinedPointValue.csv”, append=FALSE, sep= ",", row.names = FALSE, col.names=TRUE)

Step 5: You should get the results as following table.
Result: combinedPointValue.csv
 Point_ID LONGITUDE LATITUDE raster1 raster2 raster3 1 48.765863 -94.745194 200 500 -100 2 48.820552 -122.485318 178.94 18.90 10.94 3 48.233939 -107.857685 -30.74 -30.74 -0. 4 4 48.705964 -99.817363 0 110 -0.7

1. Thanks, I was have quite the frustrating afternoon trying to get ArcMap's Extract Multi Values to Points, or Sample tool to work without modifying values. I am starting to think I should just totally abandon most of the buggy tools esri has put together and just go straight to the raster package in r.

2. Michal Schnek5/05/2013 11:56 AM

Trying this with one rasterized geotiff image(elevation map) and few long/lat points in a textfile.
It correctly retrieves the elevation information for them, however when I try to create a table with coordinates and respective elevation, instead of coordinates I'll get only "S4 object of class "SpatialPoints"" in each row, followed by elevation.
What should I do to get the nice table with all columns retained from both files?

1. Hi Mike,
aashisDOTlamsal AT sdstateDOTedu? I better look it before giving the answer.

Thanks,
Ashis

3. Hi,

Good Day!

what if i have millions of rows?

i tried to process the same but im getting only 380K as a result. Don't know the status about the rest.

Kindly help on this.

4. Thanks Mike,

Worked like a dream, so quick and easy! Much appreciated.

Helen