【问题标题】:How to add the extracted values from raster stack to the data.frame of the Spatial object using terra package?如何使用 terra 包将栅格堆栈中提取的值添加到 Spatial 对象的 data.frame?
【发布时间】:2021-09-17 19:58:02
【问题描述】:

我想使用terra 包将光栅堆栈中提取的值添加到空间对象的data.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 字段不存在,这是我们过去使用 raster 包的 sp aurgument 获得的。现在如何在提取的 data.frame 中添加val 字段?

【问题讨论】:

    标签: r r-raster terra


    【解决方案1】:

    只要你是用点提取(或用线和多边形的汇总函数),你可以使用cbindterra

    示例数据

    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
    

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

    【讨论】:

      【解决方案2】:

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

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

      编辑

      提取栅格值的最快选项是exactextract 包中的exact_extract,它有一个参数(称为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")
      

      【讨论】:

      • 感谢您的回答。我知道这是可以做到的。但我想知道如何使用extract 函数本身来实现它,因为它可以在raster::extract 中使用。
      • 我不认为terra::extract 有这样的论点。我想如果你想使用terra 而不是raster,那是因为你有速度问题。在这种情况下,我会选择exactextract 包中的exact_extract,这绝对是最快的栅格提取解决方案,并且有一个参数(称为include_cols)将矢量场附加到由萃取。唯一的问题是它只接受多边形作为输入,但最简单的解决方案是在点周围创建一个小缓冲区并用函数(例如meanmedian)汇总值。
      • 如果您对此解决方案感兴趣,我可以编辑我的答案以包含此类案例
      • 您也可以提供该选项。
      猜你喜欢
      • 2021-09-08
      • 2013-01-18
      • 2016-04-09
      • 1970-01-01
      • 2022-07-13
      • 1970-01-01
      • 2014-11-01
      • 2021-11-30
      • 1970-01-01
      相关资源
      最近更新 更多