0

I want to add the extracted values from raster stack to the data.frame of the Spatial object using terra package

f <- system.file("ex/logo.tif", package="terra")
r <- rast(f)
#Plot the raster
plot(r, 1:3)

#Create a vector file
points <- cbind.data.frame(Longitude = 40, Latitude = 40, val = 0.5)
points <- rbind.data.frame(points-1, points, points+1)
points_vect <- vect(points, geom=c("Longitude", "Latitude"), crs=crs(r, proj=T),
     type = "points")

#Extract the values from raster using the point vector
terra::extract(r, points_vect, xy = T, method = "simple")

#>   ID red green blue  x  y
#> 1  1 255   255  253 30 30
#> 2  2 149   159  186 40 40
#> 3  3 146   156  207 50 50

As you can see from the output val filed is not there which we used to get using sp aurgument of raster package. Now how can I have the val field added in the extracted data.frame?

UseR10085
  • 7,120
  • 3
  • 24
  • 54

2 Answers2

4

As long as you are extracting with points (or use a summarizing function with lines and polygons), you can use bind=TRUE

Example data

f <- system.file("ex/logo.tif", package="terra")
r <- rast(f)
points <- data.frame(Longitude = 40, Latitude = 40, val = 0.5)
points <- rbind(points-1, points, points+1)
points_vect <- vect(points, geom=c("Longitude", "Latitude"), crs=crs(r))
points_vect
# class       : SpatVector 
# geometry    : points 
# dimensions  : 3, 3  (geometries, attributes)
# extent      : 39, 41, 39, 41  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
# names       : Longitude Latitude   val
# type        :     <num>    <num> <num>
# values      :        39       39  -0.5
#                      40       40   0.5
#                      41       41   1.5

Solution

v <- terra::extract(r, points_vect, xy = T, method = "simple", bind=T)
v
# class       : SpatVector 
# geometry    : points 
# dimensions  : 3, 6  (geometries, attributes)
# extent      : 39, 41, 39, 41  (xmin, xmax, ymin, ymax)
# coord. ref. : Cartesian (Meter) 
# names       :   val   red green  blue     x     y
# type        : <num> <int> <int> <int> <num> <num>
# values      :  -0.5   247   254   255  39.5  38.5
#                 0.5   149   159   186  40.5  39.5
#                 1.5   132   141   182  41.5  40.5

Alternative approaches:

First extract

v <- terra::extract(r, points_vect, xy = T, method = "simple")
v
#  ID red green blue  x  y
#1  1 247   254  255 39 39
#2  2 149   159  186 40 40
#3  3 132   141  182 41 41

Then use cbind to combine the extracted values v with points_vect

x <- cbind(points_vect, v)
x     
# class       : SpatVector 
# geometry    : points 
# dimensions  : 3, 9  (geometries, attributes)
# extent      : 39, 41, 39, 41  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
# names       : Longitude Latitude   val    ID   red green  blue     x     y
# type        :     <num>    <num> <num> <num> <num> <num> <num> <num> <num>
# values      :        39       39  -0.5     1   247   254   255    39    39
#                      40       40   0.5     2   149   159   186    40    40
#                      41       41   1.5     3   132   141   182    41    41

Or replace the original values like this

values(points_vect) <- v
points_vect
# class       : SpatVector 
# geometry    : points 
# dimensions  : 3, 6  (geometries, attributes)
# extent      : 39, 41, 39, 41  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
# names       :    ID   red green  blue     x     y
# type        : <num> <num> <num> <num> <num> <num>
# values      :     1   247   254   255    39    39
#                   2   149   159   186    40    40
#                   3   132   141   182    41    41

You can also create a new SpatVector like this

newpts <- vect(v, c("x", "y"), crs=crs(r))
newpts
# class       : SpatVector 
# geometry    : points 
# dimensions  : 3, 6  (geometries, attributes)
# extent      : 39, 41, 39, 41  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
# names       :    ID   red green  blue     x     y
# type        : <num> <num> <num> <num> <num> <num>
# values      :     1   247   254   255    39    39
#                   2   149   159   186    40    40
#                   3   132   141   182    41    41
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
2

Since terra::extract extract values with the order of elements in the vector, you can simply cbind the field to the extracted dataframe or add the needed field as a column with $

ex <- terra::extract(r, points_vect, xy = T, method = "simple")
ex$val <- points_vect$val

EDIT

The fastest option for extracting raster values is exact_extract in the exactextract package, which has an argument (called append_cols) to attach the vector field to the data.frame resulting from the extraction. exact_extract accepts only polygon data (as sf object) as input, but you can get around this by applying a buffer to your points and summarise the values. See the example below

library(exactextractr)
library(sf)

point_sf <- st_as_sf(points_vect)
point_sf <- st_buffer(point_sf,1)
ex <- exactextractr::exact_extract(stack(r),point_sf,fun="mean",append_cols="val")
#you can summarise values using a vector of functions
ex <- exactextractr::exact_extract(stack(r),point_sf,fun=c("mean","median","max","min"),append_cols="val")
Elia
  • 2,210
  • 1
  • 6
  • 18
  • Thanks for your answer. I know that can be done. But I wanted that how it can be implemented using the `extract` function itself as it is available with `raster::extract`. – UseR10085 Jul 07 '21 at 13:39
  • I don't think `terra::extract` has such an argument. I guess if you want to use `terra` instead of `raster`, it is because you have speed issues. In that case, I would opt for `exact_extract` in the `exactextract` package, which is absolutely the fastest solution for raster extraction and has an argument (called `include_cols`) to attach the vector field to the data.frame resulting from the extraction. The only problem is that it only accepts polygon as input, but the simplest solution is to create a small buffer around points and summarise values with a function (such as `mean` or `median`). – Elia Jul 07 '21 at 13:53
  • if you are interested in this solution I can edit my answer to include such a case – Elia Jul 07 '21 at 13:53
  • You can provide that option also. – UseR10085 Jul 07 '21 at 14:15