在R中沿着边界寻找一组等间距的垂直线

qni6mghb  于 2022-12-20  发布在  其他
关注(0)|答案(1)|浏览(118)

我在R中有一个描述一个州的各个县的 Dataframe ,我想要的结果是一组垂直线(不需要实际的线串,只想知道一个点的平分线的方向),这些垂直线沿着县的边界均匀分布。
我还有几个额外的考虑,一个是我不知道州的外部边界是什么,另一个是我希望这些垂线沿着“唯一”边界,也就是说,如果县A和县B共享一个边界,我只想从它们共享的边界上采样等分线一次。
现在,我正在使用st_boundary(st_combine(counties))来获取跨州的唯一边界线。
为了生成平分线,我尝试了两种方法,但运气不佳。第一种方法是沿着这条线采样点,将每一个连续的点连接成一条线,然后找到垂直于所有这些线的线。这种方法的问题是,随机采样点会从一个县跳到另一个县。所以连接一个随机采样点和下一个随机采样点的线不一定反映县的实际边界。
我的第二个想法是使用st_segmentize(),但我也遇到了麻烦。也就是说,它似乎在分割多线字符串时遇到了问题,这是在县上使用st_combine()的结果。我也很难控制线段的最大长度。
任何方向,意见,或建议将不胜感激!

jucafojl

jucafojl1#

我的方法基于基本几何学:
几个边界而不是县:

bb <- osmdata::getbb("Oborniki Śląskie")
osmdata::set_overpass_url("https://overpass-api.de/api/interpreter")
boundaries <- osmdata::opq(bb, timeout = 25 * 100) |>
  osmdata::add_osm_feature(key = "boundary", value = "administrative") |>
  osmdata::add_osm_feature(key = "admin_level", value = "8") |>
  osmdata::osmdata_sf() |>
  osmdata::unique_osmdata()

boundaries$osm_multipolygons |>
  sf::st_geometry() |>
  plot()

让我们只取两个县,找到它们的公共边界,然后使用st_intersects()来查找具有公共边界的县对。

b <- boundaries$osm_multipolygons |>
  subset(name %in% c("Lubnów", "Jary")) |>
  subset(select = "geometry") |>
  sf::st_transform(crs = "EPSG:2180") |>
  sf::st_intersection()

plot(b$geometry, col = "gray")

由于交集由2个多边形和1个多线串组成,我们只需要选择多线串。我们将它合并到LINESTRING中,并且只选择几何列。

b <- b |>
  subset(sf::st_geometry_type(b) == "MULTILINESTRING") |>
  sf::st_line_merge() |>
  sf::st_geometry()

b |>
  plot(axes = TRUE)

让我们在边界上创建几个采样点。为此,我们将使用带密度参数的st_line_sample()函数。如果您希望有N个点,请查看帮助(计算线串的长度并相应地间隔点)

sample_points <- sf::st_line_sample(b, density = 0.001, type = "regular") |>
  sf::st_cast(to = "POINT")
sample_points |>
  plot(add = TRUE)

现在我们将边界(LINESTRING)转换为对应于节点的POINTs的表。

points_from_boundary <- b |>
  sf::st_coordinates() |>
  as.data.frame() |>
  sf::st_as_sf(coords = c("X", "Y"), crs = "EPSG:2180")

让我们关注第一个采样点:

p1 <- sample_points[1]
plot(p1, pch = 20, col = "red", add = TRUE)

让我们从线串点中找出最近的点

p2 <- points_from_boundary$geometry[sf::st_nearest_feature(sample_points[1], points_from_boundary)]
plot(p2, pch = 20, col = "blue", add = TRUE)

让我们计算红点和蓝点之间的线的Angular :

alpha <- 90 - 180 * atan(
  (sf::st_coordinates(p1)[2] - sf::st_coordinates(p2)[2]) / (sf::st_coordinates(p1)[1] - sf::st_coordinates(p2)[1])
) / pi

让我们创建一个虚拟点P3,距离P1 200米,垂直于alpha旋转。

y1 <- 200 * cos(90 + alpha)
x1 <- 200 * sin(90 + alpha)

p3 <- p1 + c(x1, y1)
sf::st_crs(p3) <- sf::st_crs(b)
plot(p3, pch = 20, col = "green", add = TRUE)

sf::st_cast(sf::st_union(p1, p3), "LINESTRING") |>
  plot(add = TRUE, col = "green")

绿色线垂直于第一个采样点的边界。现在你必须将它包裹起来才能起作用,并应用于所有采样点和所有边界。
问候你格热戈兹

相关问题