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?
Follow the steps:
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:
pointCoordinates=read.csv(“pointfile.csv”)
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
|
