R 使用 tmaptools 为 sf 对象中的每一行创建边界框

R Creating boundary boxes with tmaptools for every row in an sf object

我正在尝试研究如何根据 sf 对象中每一行的 sf 点几何数据创建边界框。我正在尝试使用 tmaptools 包中的 'bb' 函数,以及 dplyr 和 rowwise()。但是,我得到的输出只是复制到每一行的相同边界框值,而不是根据每一行的点数据计算的特定边界框。

这是数据框 sf 对象的片段:

df1<-structure(list(Altitude = c(65.658, 65.606, 65.562, 65.51, 65.479, 
65.408, 65.342, 65.31, 65.242, 65.17, 65.122), Bearing = c(201.3042, 
201.3042, 201.3042, 201.3042, 201.3042, 201.3042, 201.3042, 201.3042, 
201.3042, 201.3042, 201.3042), TAI = c(0.7535967, 0.7225685, 
0.7142722, 0.7686105, 0.760403, 0.7515627, 0.7905218, 0.6231222, 
0.7246232, 0.6290409, 0.635797), lat_corrd = c(51.28648265, 51.28647067, 
51.28646118, 51.28645183, 51.28644244, 51.28643067, 51.28642109, 
51.28641164, 51.28639984, 51.2863905, 51.28638087), lon_corrd = c(0.866623929, 
0.866616631, 0.866610968, 0.86660517, 0.866598889, 0.866591258, 
0.866585183, 0.866579259, 0.866571906, 0.86656599, 0.866560288
), geometry = structure(list(structure(c(0.866623929, 51.28648265
), class = c("XY", "POINT", "sfg")), structure(c(0.866616631, 
51.28647067), class = c("XY", "POINT", "sfg")), structure(c(0.866610968, 
51.28646118), class = c("XY", "POINT", "sfg")), structure(c(0.86660517, 
51.28645183), class = c("XY", "POINT", "sfg")), structure(c(0.866598889, 
51.28644244), class = c("XY", "POINT", "sfg")), structure(c(0.866591258, 
51.28643067), class = c("XY", "POINT", "sfg")), structure(c(0.866585183, 
51.28642109), class = c("XY", "POINT", "sfg")), structure(c(0.866579259, 
51.28641164), class = c("XY", "POINT", "sfg")), structure(c(0.866571906, 
51.28639984), class = c("XY", "POINT", "sfg")), structure(c(0.86656599, 
51.2863905), class = c("XY", "POINT", "sfg")), structure(c(0.866560288, 
51.28638087), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
"sfc"), precision = 0, bbox = structure(c(xmin = 0.866560288, 
ymin = 51.28638087, xmax = 0.866623929, ymax = 51.28648265), class = "bbox"), crs = structure(list(
    input = "EPSG:4326", wkt = "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"), class = "crs"), n_empty = 0L)), row.names = 5000:5010, class = c("sf", 
"data.frame"), sf_column = "geometry", agr = structure(c(Altitude = NA_integer_, 
Bearing = NA_integer_, TAI = NA_integer_, lat_corrd = NA_integer_, 
lon_corrd = NA_integer_), .Label = c("constant", "aggregate", 
"identity"), class = "factor"))

str(df1)
#Classes ‘sf’ and 'data.frame': 11 obs. of  6 variables:

我想做的是根据 'geometry' 列中的 sfc_point 值创建一个新的边界框,例如:

require(tmaptools)
bb(df1, cx = st_coordinates(df1)[,1], cy = st_coordinates(df1)[,2], width = 0.000012, height = 0.000012, relative = FALSE)
#      xmin       ymin       xmax       ymax 
# 0.8666179 51.2864767  0.8666299 51.2864887

或者更具体地说,是这样的:

bb(df1[i,], cx = st_coordinates(df1)[i,1], cy = st_coordinates(df1)[i,2], width = 0.000012, height = 0.000012, relative = FALSE) 

我想要为每一行计算的结果 xmin、ymin、xmax、ymax 值作为一个新的几何形状,称为 boundary_boxes 添加到现有数据框中。

我尝试使用 'apply' 来做到这一点,但它根本不起作用,而且 'apply' 传递 'geometry' sf 列表类型值的方式似乎是' 正确 'bb'。 然后我尝试使用 'dplyr' 并且效果更好,但仍然不太正确:

df1 %>% rowwise() %>% mutate(boundary_boxes = list(bb(cx = st_coordinates(.)[,1], cy = st_coordinates(.)[,2], width = 0.000012, height = 0.000012, relative = FALSE)))

这几乎可行,但只是为新 'boundary_box' 列中的每一行重复相同的值,而不是给出特定的边界框。

我如何让它工作,或者有更好的方法吗? 非常感谢

一旦我为每行数据创建了一个边界框,我就需要将边界框转换为空间对象。我用 'bb_poly' 来做到这一点:

bb_poly('some boundary box data', steps = 1)

您通过在数据框中存储向量而使它不必要地复杂化。 最好将此数据存储在 tibblexy 列中。 见下文。

我在这里使用 tribblebind_cols 以使其更显眼。

library(tidyverse)
library(tmaptools)

geom = tribble(
  ~x, ~y,
  0.866623929, 51.28648265, 
  0.866616631, 51.28647067,
  0.866610968, 51.28646118, 
  0.86660517, 51.28645183, 
  0.866598889, 51.28644244, 
  0.866591258, 51.28643067, 
  0.866585183, 51.28642109, 
  0.866579259, 51.28641164, 
  0.866571906, 51.28639984, 
  0.86656599, 51.2863905,
  0.866560288, 51.28638087)

df = tibble(
  Altitude = c(65.658, 65.606, 65.562, 65.51, 65.479,65.408, 65.342, 65.31, 65.242, 65.17, 65.122), 
  Bearing = c(201.3042, 201.3042, 201.3042, 201.3042, 201.3042, 201.3042, 201.3042, 201.3042,201.3042, 201.3042, 201.3042), 
  TAI = c(0.7535967, 0.7225685, 0.7142722, 0.7686105, 0.760403, 0.7515627, 0.7905218, 0.6231222, 0.7246232, 0.6290409, 0.635797), 
  lat_corrd = c(51.28648265, 51.28647067, 51.28646118, 51.28645183, 51.28644244, 51.28643067, 51.28642109, 51.28641164, 51.28639984, 51.2863905, 51.28638087), 
  lon_corrd = c(0.866623929, 0.866616631, 0.866610968, 0.86660517, 0.866598889, 0.866591258, 0.866585183, 0.866579259, 0.866571906, 0.86656599, 0.866560288), 
) %>% bind_cols(geom) 
df

输出

# A tibble: 11 x 7
   Altitude Bearing   TAI lat_corrd lon_corrd     x     y
      <dbl>   <dbl> <dbl>     <dbl>     <dbl> <dbl> <dbl>
 1     65.7    201. 0.754      51.3     0.867 0.867  51.3
 2     65.6    201. 0.723      51.3     0.867 0.867  51.3
 3     65.6    201. 0.714      51.3     0.867 0.867  51.3
 4     65.5    201. 0.769      51.3     0.867 0.867  51.3
 5     65.5    201. 0.760      51.3     0.867 0.867  51.3
 6     65.4    201. 0.752      51.3     0.867 0.867  51.3
 7     65.3    201. 0.791      51.3     0.867 0.867  51.3
 8     65.3    201. 0.623      51.3     0.867 0.867  51.3
 9     65.2    201. 0.725      51.3     0.867 0.867  51.3
10     65.2    201. 0.629      51.3     0.867 0.867  51.3
11     65.1    201. 0.636      51.3     0.867 0.867  51.3

现在您只需要一个简单的函数,即 returns 您以 tibble 的形式得到函数 bb 的结果。比如像下面这样。

fbb = function(data){
  bbout = bb(cx=data$x, cy=data$y, width= 0.000012, height= 0.000012, relative= FALSE)
  tibble(
    xmin = bbout["xmin"],
    ymin = bbout["ymin"],
    xmax = bbout["xmax"],
    ymax = bbout["ymax"]
  )
}

下一个再简单不过了。

df %>% 
  nest(data=x:y) %>% 
  mutate(bbout = map(data, fbb)) %>% 
  unnest(c(data, bbout))

输出

# A tibble: 11 x 11
   Altitude Bearing   TAI lat_corrd lon_corrd     x     y  xmin  ymin  xmax  ymax
      <dbl>   <dbl> <dbl>     <dbl>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1     65.7    201. 0.754      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 2     65.6    201. 0.723      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 3     65.6    201. 0.714      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 4     65.5    201. 0.769      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 5     65.5    201. 0.760      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 6     65.4    201. 0.752      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 7     65.3    201. 0.791      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 8     65.3    201. 0.623      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 9     65.2    201. 0.725      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
10     65.2    201. 0.629      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
11     65.1    201. 0.636      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3

或者你说每一行的数据都一样?当然不是。看

df %>%
  nest(data=x:y) %>%
  mutate(bbout = map(data, fbb)) %>%
  unnest(c(data, bbout)) %>% 
  ggplot(aes(x, y)) +
  geom_ribbon(aes(ymin=ymin, ymax=ymax), alpha=0.2)+
  geom_line(aes(x, y))

再添加一项功能fbb_poly,您将拥有盒子!!

fbb_poly = function(data) 
  bb_poly(bb(cx=data$x, cy=data$y, width= 0.000012, 
             height= 0.000012, relative= FALSE), steps = 1)


df %>%
  nest(data=x:y) %>%
  mutate(bbout = map(data, fbb),
         bb_poly = map(data, fbb_poly)) %>%
  unnest(c(data, bb_poly)) 

输出

# A tibble: 11 x 9
   Altitude Bearing   TAI lat_corrd lon_corrd     x     y bbout                                                                                                 bb_poly
      <dbl>   <dbl> <dbl>     <dbl>     <dbl> <dbl> <dbl> <list>                                                                                 <POLYGON [arc_degree]>
 1     65.7    201. 0.754      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8666179 51.28648, 0.8666179 51.28649, 0.8666299 51.28649, 0.8666299 51.28648, 0.86661...
 2     65.6    201. 0.723      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8666106 51.28646, 0.8666106 51.28648, 0.8666226 51.28648, 0.8666226 51.28646, 0.86661...
 3     65.6    201. 0.714      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.866605 51.28646, 0.866605 51.28647, 0.866617 51.28647, 0.866617 51.28646, 0.866605 51...
 4     65.5    201. 0.769      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8665992 51.28645, 0.8665992 51.28646, 0.8666112 51.28646, 0.8666112 51.28645, 0.86659...
 5     65.5    201. 0.760      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8665929 51.28644, 0.8665929 51.28645, 0.8666049 51.28645, 0.8666049 51.28644, 0.86659...
 6     65.4    201. 0.752      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8665853 51.28642, 0.8665853 51.28644, 0.8665973 51.28644, 0.8665973 51.28642, 0.86658...
 7     65.3    201. 0.791      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8665792 51.28642, 0.8665792 51.28643, 0.8665912 51.28643, 0.8665912 51.28642, 0.86657...
 8     65.3    201. 0.623      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8665733 51.28641, 0.8665733 51.28642, 0.8665853 51.28642, 0.8665853 51.28641, 0.86657...
 9     65.2    201. 0.725      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8665659 51.28639, 0.8665659 51.28641, 0.8665779 51.28641, 0.8665779 51.28639, 0.86656...
10     65.2    201. 0.629      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.86656 51.28638, 0.86656 51.2864, 0.866572 51.2864, 0.866572 51.28638, 0.86656 51.28638))
11     65.1    201. 0.636      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8665543 51.28637, 0.8665543 51.28639, 0.8665663 51.28639, 0.8665663 51.28637, 0.86655...