JLD475
JLD475

Reputation: 179

Why am I losing columns when extracting raster data with ```terra::extract```?

I created a raster stack with rasters containing different vegetation measurements (i.e., canopy height, veg density). I extracted data from that raster stack to a SpatVector containing GPS points and corresponding data. The output contains the raster data, but does not contain any of the SpatVector data. Sample code below. I am not sure how to add the raster data to the question.

structure(list(Id = c("A1", "A1", "A1", "A1", "A1", "A1", "A1", 
"A1", "A1", "A1"), DateTime_Local = c("2019-06-18 14:00:00", 
"2019-06-18 14:30:00", "2019-06-18 15:00:00", "2019-06-18 15:30:00", 
"2019-06-18 16:00:00", "2019-06-18 16:30:00", "2019-06-18 17:00:00", 
"2019-06-18 17:30:00", "2019-06-18 18:00:00", "2019-06-18 18:30:00"
), Temp_C = c(23.484, 23.388, 23.196, 23.677, 24.738, 24.738, 
24.641, 26.097, 27.37, 28.357), Temp_F = c(74.2712, 74.0984, 
73.7528, 74.6186, 76.5284, 76.5284, 76.3538, 78.9746, 81.266, 
83.0426), Type = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1), Long = c(-97.47462153, 
-97.47462153, -97.47462153, -97.47462153, -97.47462153, -97.47462153, 
-97.47462153, -97.47462153, -97.47462153, -97.47462153), Lat = c(26.58459955, 
26.58459955, 26.58459955, 26.58459955, 26.58459955, 26.58459955, 
26.58459955, 26.58459955, 26.58459955, 26.58459955), Long.1 = c(651903.662642045, 
651903.662642045, 651903.662642045, 651903.662642045, 651903.662642045, 
651903.662642045, 651903.662642045, 651903.662642045, 651903.662642045, 
651903.662642045), Lat.1 = c(2941332.22211244, 2941332.22211244, 
2941332.22211244, 2941332.22211244, 2941332.22211244, 2941332.22211244, 
2941332.22211244, 2941332.22211244, 2941332.22211244, 2941332.22211244
)), row.names = c(NA, -10L), class = "data.frame")
BG_vect <- vect(BG.sf) #SF object containing GPS coordinates and point data

BG.extracted <- terra::extract(veg_stk, BG_vect, fun = mean)
summary(BG.extracted)

Upvotes: 0

Views: 664

Answers (2)

Robert Hijmans
Robert Hijmans

Reputation: 47146

You can combine the original point data with the extracted values using cbind.

Example data

library(terra)
df <- data.frame(id=1:5, Long = 1:5, Lat=1:5, var=letters[1:5])
df
#  id Long Lat var
#1  1    1   1   a
#2  2    2   2   b
#3  3    3   3   c
#4  4    4   4   d
#5  5    5   5   e

r <- rast(xmin=0, xmax=6, ymin=0, ymax=6, nlyr=2, res=.5, names=c("A", "B"))
set.seed(0)
values(r) <- sample(size(r))

It would be efficient to extract the values from the raster directly with the x (longitude) and y (latitude) coordinates in the data.frame

e1 <- extract(r, df[, c("Long", "Lat")])
e1
#  ID   A   B
#1  1 273  46
#2  2 201  18
#3  3  51 238
#4  4 141  27
#5  5 115 106
 

But you can also first create a SpatVector

v <- vect(df, c("Long", "Lat"))
e2 <- extract(r, v)

And in either case, you can cbind the results to the data.frame or the SpatVector.

cbind(df, e1[,-1])
#  id Long Lat var   A   B
#1  1    1   1   a 273  46
#2  2    2   2   b 201  18
#3  3    3   3   c  51 238
#4  4    4   4   d 141  27
#5  5    5   5   e 115 106

cbind(v, e2[,-1])
# class       : SpatVector 
# geometry    : points 
# dimensions  : 5, 4  (geometries, attributes)
# extent      : 1, 5, 1, 5  (xmin, xmax, ymin, ymax)
# coord. ref. :  
# names       :    id   var     A     B
# type        : <int> <chr> <int> <int>
# values      :     1     a   273    46
#                   2     b   201    18
#                   3     c    51   238

So to answer your question: you are not loosing columns; it is just that the input data is not replicated in the output. Also, while you could use merge here, that would be rather inefficient.

Upvotes: 1

Sean McKenzie
Sean McKenzie

Reputation: 909

I think all you need to do is merge the resulting data.framefrom terra::extract() back to your SpatVector object. I created a reproducible example with your data. Note, your spatial data only seems to include one point location, so I altered the "Lat" and "Long" columns with some randomly selected locations. I also assumed these data were in the WGS84 Lat/Long coordinate system (EPSG:4269). From that I created a fake raster of canopy height data.

library(terra)
library(sf)

spvect<-structure(list(Id = c("A1", "A1", "A1", "A1", "A1", "A1", "A1", 
                              "A1", "A1", "A1"), DateTime_Local = c("2019-06-18 14:00:00", 
                                                                    "2019-06-18 14:30:00", "2019-06-18 15:00:00", "2019-06-18 15:30:00", 
                                                                    "2019-06-18 16:00:00", "2019-06-18 16:30:00", "2019-06-18 17:00:00", 
                                                                    "2019-06-18 17:30:00", "2019-06-18 18:00:00", "2019-06-18 18:30:00"
                              ), Temp_C = c(23.484, 23.388, 23.196, 23.677, 24.738, 24.738, 
                                            24.641, 26.097, 27.37, 28.357), Temp_F = c(74.2712, 74.0984, 
                                                                                       73.7528, 74.6186, 76.5284, 76.5284, 76.3538, 78.9746, 81.266, 
                                                                                       83.0426), Type = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1), Long = c(-97.47462153, 
                                                                                                                                                  -97.47462153, -97.47462153, -97.47462153, -97.47462153, -97.47462153, 
                                                                                                                                                  -97.47462153, -97.47462153, -97.47462153, -97.47462153), Lat = c(26.58459955, 
                                                                                                                                                                                                                   26.58459955, 26.58459955, 26.58459955, 26.58459955, 26.58459955, 
                                                                                                                                                                                                                   26.58459955, 26.58459955, 26.58459955, 26.58459955), Long.1 = c(651903.662642045, 
                                                                                                                                                                                                                                                                                   651903.662642045, 651903.662642045, 651903.662642045, 651903.662642045, 
                                                                                                                                                                                                                                                                                   651903.662642045, 651903.662642045, 651903.662642045, 651903.662642045, 
                                                                                                                                                                                                                                                                                   651903.662642045), Lat.1 = c(2941332.22211244, 2941332.22211244, 
                                                                                                                                                                                                                                                                                                                2941332.22211244, 2941332.22211244, 2941332.22211244, 2941332.22211244, 
                                                                                                                                                                                                                                                                                                                2941332.22211244, 2941332.22211244, 2941332.22211244, 2941332.22211244
                                                                                                                                                                                                                                                                                   )), row.names = c(NA, -10L), class = "data.frame")
spvect$Long<-runif(nrow(spvect), -97.5, -96.5)
spvect$Lat<-runif(nrow(spvect), 26, 27)
BG.sf<-sf::st_as_sf(spvect, coords=c("Long", "Lat"), crs=4269)
BG.sf[,"Ind"]<-rownames(BG.sf)
BG.vect<-vect(BG.sf)
rst<-rast(extent=ext(BG.vect), nrow=100, ncol=100,  crs=crs(BG.vect))
values(rst)<-rnorm(10000, 100, 12)
names(rst)<-"Canopy Height"
extrctd<-extract(rst,BG.vect)
BG.Final<-terra::merge(BG.vect, extrctd, by.x="Ind", by.y="ID")

Upvotes: 1

Related Questions