UseR10085
UseR10085

Reputation: 8136

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

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?

Upvotes: 0

Views: 2016

Answers (2)

Robert Hijmans
Robert Hijmans

Reputation: 47036

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

Upvotes: 6

Elia
Elia

Reputation: 2584

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")

Upvotes: 3

Related Questions