0

我想使用包将光栅堆栈中提取的值添加到 Spatial 对象的terradata.frame

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

正如您从输出val文件中看到的那样，我们过去使用sp包的参数来获得它不存在raster。现在如何val在提取的 data.frame 中添加字段？

4

2 回答 2

2

由于terra::extract使用向量中元素的顺序提取值，您可以简单地cbind将字段添加到提取的数据框中，或者将所需字段添加为列$

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

编辑

提取栅格值的最快选项是exact_extract在exactextract包中，它有一个参数（称为append_cols）将矢量场附加到提取后的 data.frame。 exact_extract仅接受多边形数据（作为sf对象）作为输入，但您可以通过将缓冲区应用于您的点并汇总值来解决此问题。请参阅下面的示例

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")
于 2021-07-07T13:34:50.083 回答
2

只要您使用点进行提取（或对线和多边形使用汇总函数），就可以使用cbindwithterra

示例数据

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

利用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

并将提取的值v与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

或者像这样替换原始值

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

SpatVector或者像这样创建一个新的

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
于 2021-07-07T16:10:36.253 回答