Extracting values from rasters at location of points using R

Question

My question is about extracting values from rasters at the location of points. With the function extract this is very easy, and the function gives me a dataframe with the values of all the variables in the points. I want to have in that dataframe the coordinates of each point.

How can I make that happen, perhaps by saying to R that when extracting the values from the raster also add the columns of the location point?

This is my procedure:

presencias=read.table("c:/SDM_R/presencias/P_lentiscus_pres.csv",header=TRUE,sep=";")

lista_variables <-list.files(path="Variables_modelizacion/solo_ascii",pattern='*.asc',full.names=TRUE)
variables <- stack(lista_variables)

variables_presencia<-extract(variables, presencias)

the results are something like these:

> bio1  bio12  bio18  bio2  bio4 
> 90    875    165    95    4886
> 115   1085   158    83    4075
> 135   1153   153    67    3402
> 85    1026   137    99    5203
> 96    667    128    108   5823
> 98    531    109    113   6305
> 132   450    63     123   6598
> 132   569    104    106   5963
> 95    814    196    98    5571
> 146   474    39     114   6603

But I want two more columns with the coordinate data (but could be others columns from my csv extracting location table).

Answer 1

Assuming that presencias and variables share the same projection, this should be an easy task. I recommend you to add these lines of code after your read.table() statement in order to convert presencias dataframe to a SpatialPointsDataFrame object (just refine the names of the columns containing x and y coordinates if they differ from my example).

coordinates(presencias) <- c("x", "y")

To provide a reproducible example, I try to open up the scope of my answer a little more. First of all, download and unzip this ESRI shapefile with more or less important locations in Germany. These will serve as point data later on. You will also need packages dismo, rgdal and raster for this short example, so make sure that these libraries (and all their dependencies) are installed on your local hard drive.

Let's start with loading the required packages.

library(dismo)
library(rgdal)
library(raster)

Next, you should generate a sample RasterLayer. In our case, we will make use of the gmap() function from the dismo package in order to obtain a physical map of Germany.

germany.mrc <- gmap("Germany")

You can now import your point shapefile via readOGR from R's rgdal package. Make sure to adjust the data source name (dsn = ...). The whole projection stuff is obsolete in your particular case. However, it has to be done in our example in order to successfully overlay our point data with the Germany RasterLayer.

# Import SpatialPointsDataFrame
germany.places <- readOGR(dsn = "/path/to/shapefile", layer = "places")
# Define shapefile's current CRS
projection(germany.places) <- CRS("+proj=lonlat +ellps=WGS84")
# Reproject to RasterLayer's CRS
germany.places.mrc <- spTransform(germany.places, CRS(projection(germany.mrc)))

To reduce the huge size of our point data, we will draw a random sample of ten locations in Germany. This should suffice for our purposes.

set.seed(35)
germany.places.mrc.sample <- germany.places.mrc[sample(nrow(germany.places.mrc), 10), ]

Now that the preparation stuff is finished, we could just start to extract the values of those particular pixels our ten randomly sampled points lie within.

data <- data.frame(coordinates(germany.places.mrc.sample),
                   germany.places.mrc.sample$name, 
                   extract(germany.mrc, germany.places.mrc.sample))
names(data) <- c("x", "y", "name", "value")

In order to merge the point coordinates with the extracted pixel values, we just need to set up a dataframe containing the coordinates of our SpatialPointsDataFrame. That's it!

data
           x       y          name value
1  1073490.3 6513446 Veitsteinbach   208
2  1269100.8 6156690   Assenhausen   231
3  1336757.5 6246284    Frauenwahl   195
4   828579.9 6634122      Altenhof   189
5  1571418.1 6662558         Wohla   151
6  1192299.4 6864087     Flechtorf   170
7   976270.0 6362050    Hilsenhain   208
8  1117416.4 6092146      Nestbaum   175
9  1274192.0 6344490 Wappeltshofen   236
10  878488.2 6839843        Leeden   208

Answer 2

Surely you can just do:

variables_presencia$x <- presencias['x']

variables_presencia$y <- presencias['y']

(presuming your coordinate data is in two columns called 'x' and 'y')