R使用tmaptools为sf对象中的每一行创建边界框
Posted
技术标签:
【中文标题】R使用tmaptools为sf对象中的每一行创建边界框【英文标题】：R Creating boundary boxes with tmaptools for every row in an sf object 【发布时间】：2021-10-31 03:07:45 【问题描述】：我正在尝试研究如何从 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:
我要做的是从“几何”列中的 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 值作为新几何图形，称为边界框，添加到现有数据框。
我尝试使用 '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”列中的每一行重复相同的值，而不是给出特定的边界框。
我如何让它工作，或者有更好的方法吗？ 非常感谢
一旦每行数据都有一个 bbox，我就需要将边界框转换为空间对象。我使用 'bb_poly' 来做到这一点：
bb_poly('some boundary box data', steps = 1)
【问题讨论】：
圆形会更容易，但您想要每个点的正方形多边形？ 实际上一个圆圈可能会起作用，但椭圆形会更好。我需要能够设置椭圆的宽度和高度。如果我还可以从轴承柱设置圆（或矩形）的方向，那就更好了。 【参考方案1】：您通过将向量存储在数据框中使其变得不必要地复杂。
最好将此数据存储在具有x 和y 列的tibble 中。
见下文。
我在这里使用tribble 和bind_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
现在您只需要一个简单的函数，它以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...
【讨论】：
这很棒。它给了我需要的边界框的角点。但是，一些重要的地理空间信息会从数据结构中丢失。由“bb”函数创建的普通“bbox”将包含 CRS 数据。将 bbox 数据保存为 tibble 会丢失这一点，这会导致我需要做的下一步出现问题 - 应用“bb_poly”。我将在原始问题中添加更多信息。 我已经改变了你的函数，以便将“bbout”存储为一个列表： 我添加了几行代码。希望你现在期待它。buot 可以作为列表返回，但也可以作为tibble 返回。没关系。关键是，您应该得到更容易处理的结果。
效果很好。非常感谢您的帮助。我的下一个挑战是看看我是否可以根据“轴承”列旋转边界框。
我相信你可以！我很高兴能帮上忙。虽然我不知道“基于 'Bearing' 列旋转边界框”是什么。我真的不知道bb_poly 函数是干什么用的。但是，我已经向您展示了您应该遵循的路径。编写一个函数来执行您想要的操作，并通过 tibble 中的每个元素 map。根据你所掌握的知识分析结果。以上是关于R使用tmaptools为sf对象中的每一行创建边界框的主要内容，如果未能解决你的问题，请参考以下文章