确定点是否位于重叠多边形内的简单方法。

huangapple
go评论122阅读模式
英文:

simple approach to determine if points intersect overlapping polygons

问题

以下是您提供的内容的翻译：

我有一个包含纬度和经度坐标的简单数据框（我提供了一个子集）。我想确定这些点是否落在一个或两个多边形内（第二个多边形完全覆盖第一个）。

点的子样本

pts <- structure(list(Longitude = c(-26.775, -27.3167, -30.2717, -30.7267, -37.8683, -38.2733,
-44.3033, -44.8233, -57.42, -57.8717, -52.235, -52.54,
-54.2317, -54.6883, -63.6033, -65.415, -66.31, -66.7533),
Latitude = c(62.5183, 62.3, 60.8367, 60.59, 56.6383, 56.3917, 52.8, 52.45,
45.56, 45.46, 47.3917, 47.1317, 46.26, 46.16, 44.29, 43.2633,
43.0983, 43.0267)),
row.names = c(NA, -18L), class = c("tbl_df", "tbl", "data.frame"))

将点表示为空间对象

sp <- SpatialPoints(pts, proj4string = CRS("+proj=longlat +datum=WGS84"))

多边形 A 和 B

pA <- st_polygon(list(as.matrix(data.frame(lon = c(-40, -40, -64, -64, -54, -54, -40),
lat = c(45, 61, 61, 51, 51, 45, 45))))) %>%
st_sfc(crs = st_crs("+proj=longlat +datum=WGS84"))
pB <- st_polygon(list(as.matrix(data.frame(lon = c(-6, -6, -50, -50, -64, -64, -54, -54, -6),
lat = c(45, 66, 66, 70, 70, 51, 51, 45, 45))))) %>%
st_sfc(crs = st_crs("+proj=longlat +datum=WGS84"))

我知道可以使用 st_intersects 来确定点是否在单个多边形内重叠，但我不知道如何扩展到多个重叠的多边形。理想情况下，最终结果将是一个变量，如果点在 A 中，则返回 A，如果点属于不与 A 重叠的较大的 B 多边形，则返回 B，如果点既不在 A 也不在 B 中，则返回 NA。

我发现使用 tidyverse（而不是 datatable）的代码更容易理解，所以我希望有一个使用这些函数的解决方案。

英文:

I have a simple dataframe with latitude and longitude coordinates (I've provided a subset below). I'd like to determine if any of the points fall within one or both polygons (the second polygon completely overlaps the first).

  1. # subsample of points
  2. pts &lt;- structure(list(Longitude = c(-26.775, -27.3167, -30.2717, -30.7267, -37.8683, -38.2733, 
  3.                                     -44.3033, -44.8233, -57.42, -57.8717, -52.235, -52.54, 
  4.                                     -54.2317, -54.6883, -63.6033, -65.415, -66.31, -66.7533), 
  5.                       Latitude = c(62.5183, 62.3, 60.8367, 60.59, 56.6383, 56.3917, 52.8, 52.45, 
  6.                                    45.56, 45.46, 47.3917, 47.1317, 46.26, 46.16, 44.29, 43.2633, 
  7.                                    43.0983, 43.0267)), 
  8.                  row.names = c(NA, -18L), class = c(&quot;tbl_df&quot;, &quot;tbl&quot;, &quot;data.frame&quot;))
  10. # points as a spatial object
  11. sp &lt;- SpatialPoints(pts, proj4string = CRS(&quot;+proj=longlat +datum=WGS84&quot;))
  13. # polygon A &amp; B
  14. pA &lt;- st_polygon(list(as.matrix(data.frame(lon = c(-40, -40, -64, -64, -54, -54, -40),
  15.                                               lat = c(45, 61, 61, 51, 51, 45, 45))))) %&gt;%
  16.   st_sfc(crs = st_crs(&quot;+proj=longlat +datum=WGS84&quot;))
  17. pB&lt;- st_polygon(list(as.matrix(data.frame(lon = c(-6, -6, -50, -50, -64, -64, -54, -54, -6),
  18.                                                  lat = c(45, 66, 66, 70, 70, 51, 51, 45, 45))))) %&gt;%
  19.   st_sfc(crs = st_crs(&quot;+proj=longlat +datum=WGS84&quot;))

I know I can use st_intersects to determine if points overlap in a single polygon, but I don't know how to expand this if when you have multiple overlapping polygons. Ideally, the end result would be a variable that returns A if the point is just in A, B if the point is part of the larger B polygon that does not overlap A, and NA if the point is not in either polygon.

I find it much easier to follow code that uses tidyverse (rather than datatable), so I'd prefer a solution with those functions.

答案1

得分: 2

If you have collected your polygons into sf object and find point intersects with that, st_intersects() will return a "sparse geometry binary predicate list",它包含每个点的交叉几何索引的列表:

  1. library(sf)
  2. library(sp)
  3. library(dplyr)
  4. library(ggplot2)
  6. ### 移动到 sf 对象，假设没有部分重叠，按大小排列多边形
  7. ab_sf <- st_sf(geometry = c(pA, pB), row.names = c("pA", "pB")) %>%
  8.   tibble::rownames_to_column("poly") %>%
  9.   mutate(size_rank = rank(desc(st_area(geometry)))) %>%
  10.   arrange(size_rank)
  11. ab_sf
  12. #> Simple feature collection with 2 features and 2 fields
  13. #> Geometry type: POLYGON
  14. #> Dimension:     XY
  15. #> Bounding box:  xmin: -64 ymin: 45 xmax: -6 ymax: 70
  16. #> Geodetic CRS:  +proj=longlat +datum=WGS84
  17. #>   poly size_rank                       geometry
  18. #> 1   pB         1 POLYGON ((-6 45, -6 66, -50...
  19. #> 2   pA         2 POLYGON ((-40 45, -40 61, -...
  21. points_sf <- st_as_sf(sp_)
  23. ### 查找交叉点，返回的列表将包括交叉索引的有序列表
  24. ### 例如仅为 pA 的 (1)；为 pA 和 pB 的 (1,2)；(空) 表示无交叉，
  25. ### 从每个中提取最后一项，以获取最小交叉多边形的索引
  26. isect_last_idx <- st_intersects(points_sf, ab_sf) %>% sapply(last)
  27. isect_last_idx
  28. #>  [1]  1  1  1  1  1  1  2  2 NA NA  2  2 NA NA NA NA NA NA
  30. ## 将索引映射到多边形名称
  31. points_sf$poly_name <- ab_sf[["poly"]][isect_last_idx]
  32. points_sf
  33. #> Simple feature collection with 18 features and 1 field
  34. #> Geometry type: POINT
  35. #> Dimension:     XY
  36. #> Bounding box:  xmin: -66.7533 ymin: 43.0267 xmax: -26.775 ymax: 62.5183
  37. #> Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs
  38. #> First 10 features:
  39. #>                    geometry poly_name
  40. #> 1   POINT (-26.775 62.5183)        pB
  41. #> 2     POINT (-27.3167 62.3)        pB
  42. #> 3  POINT (-30.2717 60.8367)        pB
  43. #> 4    POINT (-30.7267 60.59)        pB
  44. #> 5  POINT (-37.8683 56.6383)        pB
  45. #> 6  POINT (-38.2733 56.3917)        pB
  46. #> 7     POINT (-44.3033 52.8)        pA
  47. #> 8    POINT (-44.8233 52.45)        pA
  48. #> 9      POINT (-57.42 45.56)      <NA>
  49. #> 10   POINT (-57.8717 45.46)      <NA>
  51. ab_sf %>%
  52.   arrange(size_rank) %>%
  53.   ggplot() +
  54.   geom_sf(aes(fill = poly, linetype = poly, linewidth = poly), alpha = .8) +
  55.   scale_linewidth_manual(values = c(0,2)) +
  56.   geom_sf(data = points_sf,  aes(shape = poly_name), size =  5) +
  57.   scale_shape_manual(na.value = 1, values = c(9,13)) +
  58.   theme_minimal()

确定点是否位于重叠多边形内的简单方法。

使用的数据:

  2. pts <- structure(list(Longitude = c(-26.775, -27.3167, -30.2717, -30.7267, -37.8683, -38.2733, 
  3. -44.3033, -44.8233, -57.42, -57.8717, -52.235, -52.54, 
  4. -54.2317, -54.6883, -63.6033, -65.415, -66.31, -66.7533), 
  5. Latitude = c(62.5183, 62.3, 60.8367, 60.59, 56.6383, 56.3917, 52.8, 52.45, 
  6. 45.56, 45.46, 47.3917, 47.1317, 46.26, 46.16, 44.29, 43.2633, 
  7. 43.0983, 43.0267)), 
  8. row.names = c(NA, -18L), class = c("tbl_df", "tbl", "data.frame"))
  9. # points as a spatial object
  10. sp_ <- SpatialPoints(pts, proj4string = CRS("+proj=longlat +datum=WGS84"))
  11. # polygon A & B
  12. pA <- st_polygon(list(as.matrix(data.frame(lon = c(-40, -40, -64, -64, -54, -54, -40),
  13. lat = c(45, 61, 61, 51, 51, 45, 45))))) %>%
  14. st_sfc(crs = st_crs("+proj=longlat +datum=WGS84"))
  15. pB<- st_polygon(list(as.matrix(data.frame(lon = c(-6, -6, -50, -50, -64, -64, -54, -54, -6),
  16. lat = c(45, 66, 66, 70, 70, 51, 51, 45, 45))))) %>%
  17. st_sfc(crs = st_crs("+proj=longlat +datum=WGS84"))

创建于 2023-05-25，使用 reprex v2.0.2

英文:

If you have collected your polygons into sf object and find point intersects with that, st_intersects() will return a "sparse geometry binary predicate list", it includes a list of intersecting geometry indices for each point:

  1. Sparse geometry binary predicate list of length 18, where the predicate was `intersects&#39;
  2. first 10 elements:
  3. 1: 1
  4. 2: 1
  5. 3: 1
  6. 4: 1
  7. 5: 1
  8. 6: 1
  9. 7: 1, 2
  10. 8: 1, 2
  11. 9: (empty)
  12. 10: (empty)

If your polygons are ordered correctly in sf, last item for every point will be the polygon with highest row number. Well, if you just want know if point intersects with 1, 2, (,n) or no polygons, all you really need is to get the lengths() of that list.

  1. library(sf)
  2. library(sp)
  3. library(dplyr)
  4. library(ggplot2)
  6. ### move to sf objects, assume that there are no partial overlaps,
  7. ### arrange polygons by their size
  8. ab_sf &lt;- st_sf(geometry = c(pA, pB), row.names = c(&quot;pA&quot;, &quot;pB&quot;)) %&gt;% 
  9.   tibble::rownames_to_column(&quot;poly&quot;) %&gt;% 
  10.   mutate(size_rank = rank(desc(st_area(geometry)))) %&gt;%  
  11.   arrange(size_rank)
  12. ab_sf
  13. #&gt; Simple feature collection with 2 features and 2 fields
  14. #&gt; Geometry type: POLYGON
  15. #&gt; Dimension:     XY
  16. #&gt; Bounding box:  xmin: -64 ymin: 45 xmax: -6 ymax: 70
  17. #&gt; Geodetic CRS:  +proj=longlat +datum=WGS84
  18. #&gt;   poly size_rank                       geometry
  19. #&gt; 1   pB         1 POLYGON ((-6 45, -6 66, -50...
  20. #&gt; 2   pA         2 POLYGON ((-40 45, -40 61, -...
  22. points_sf &lt;- st_as_sf(sp_)
  24. ### find intersects, returned list will include ordered list of intersecting indices
  25. ### i.e (1) for only pA; (1,2) for pA and pB; (empty) for no intersect,
  26. ### extract last item from each to get index of smallest intersecting polygon
  27. isect_last_idx &lt;- st_intersects(points_sf, ab_sf) %&gt;% sapply(last)
  28. isect_last_idx
  29. #&gt;  [1]  1  1  1  1  1  1  2  2 NA NA  2  2 NA NA NA NA NA NA
  31. ## map indeces to polygon names
  32. points_sf$poly_name &lt;- ab_sf[[&quot;poly&quot;]][isect_last_idx]
  33. points_sf
  34. #&gt; Simple feature collection with 18 features and 1 field
  35. #&gt; Geometry type: POINT
  36. #&gt; Dimension:     XY
  37. #&gt; Bounding box:  xmin: -66.7533 ymin: 43.0267 xmax: -26.775 ymax: 62.5183
  38. #&gt; Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs
  39. #&gt; First 10 features:
  40. #&gt;                    geometry poly_name
  41. #&gt; 1   POINT (-26.775 62.5183)        pB
  42. #&gt; 2     POINT (-27.3167 62.3)        pB
  43. #&gt; 3  POINT (-30.2717 60.8367)        pB
  44. #&gt; 4    POINT (-30.7267 60.59)        pB
  45. #&gt; 5  POINT (-37.8683 56.6383)        pB
  46. #&gt; 6  POINT (-38.2733 56.3917)        pB
  47. #&gt; 7     POINT (-44.3033 52.8)        pA
  48. #&gt; 8    POINT (-44.8233 52.45)        pA
  49. #&gt; 9      POINT (-57.42 45.56)      &lt;NA&gt;
  50. #&gt; 10   POINT (-57.8717 45.46)      &lt;NA&gt;
  52. ab_sf %&gt;% 
  53.   arrange(size_rank) %&gt;% 
  54.   ggplot() +
  55.   geom_sf(aes(fill = poly, linetype = poly, linewidth = poly), alpha = .8) +
  56.   scale_linewidth_manual(values = c(0,2)) +
  57.   geom_sf(data = points_sf,  aes(shape = poly_name), size =  5) +
  58.   scale_shape_manual(na.value = 1, values = c(9,13)) +
  59.   theme_minimal()

确定点是否位于重叠多边形内的简单方法。<!-- -->

Used data:

  2. pts &lt;- structure(list(Longitude = c(-26.775, -27.3167, -30.2717, -30.7267, -37.8683, -38.2733, 
  3. -44.3033, -44.8233, -57.42, -57.8717, -52.235, -52.54, 
  4. -54.2317, -54.6883, -63.6033, -65.415, -66.31, -66.7533), 
  5. Latitude = c(62.5183, 62.3, 60.8367, 60.59, 56.6383, 56.3917, 52.8, 52.45, 
  6. 45.56, 45.46, 47.3917, 47.1317, 46.26, 46.16, 44.29, 43.2633, 
  7. 43.0983, 43.0267)), 
  8. row.names = c(NA, -18L), class = c(&quot;tbl_df&quot;, &quot;tbl&quot;, &quot;data.frame&quot;))
  9. # points as a spatial object
  10. sp_ &lt;- SpatialPoints(pts, proj4string = CRS(&quot;+proj=longlat +datum=WGS84&quot;))
  11. # polygon A &amp; B
  12. pA &lt;- st_polygon(list(as.matrix(data.frame(lon = c(-40, -40, -64, -64, -54, -54, -40),
  13. lat = c(45, 61, 61, 51, 51, 45, 45))))) %&gt;%
  14. st_sfc(crs = st_crs(&quot;+proj=longlat +datum=WGS84&quot;))
  15. pB&lt;- st_polygon(list(as.matrix(data.frame(lon = c(-6, -6, -50, -50, -64, -64, -54, -54, -6),
  16. lat = c(45, 66, 66, 70, 70, 51, 51, 45, 45))))) %&gt;%
  17. st_sfc(crs = st_crs(&quot;+proj=longlat +datum=WGS84&quot;))

<sup>Created on 2023-05-25 with reprex v2.0.2</sup>

答案2

得分: 2

  1. # 将两个多边形组合成一个 sf 对象
  2. combined_poly <- st_sf(geometry = c(pA, pB), row.names = c("pA", "pB")) 
  3. # 将点转换为 sf 对象
  4. points_as_sf <- st_as_sf(sp)
  5. st_intersects(points_as_sf, combined_poly, sparse = F)
  6. #       [,1]  [,2]
  7. #  [1,] FALSE  TRUE
  8. #  [2,] FALSE  TRUE
  9. #  [3,] FALSE  TRUE
  10. #  [4,] FALSE  TRUE
  11. #  [5,] FALSE  TRUE
  12. #  [6,] FALSE  TRUE
  13. #  [7,]  TRUE  TRUE
  14. #  [8,]  TRUE  TRUE
  15. #  [9,] FALSE FALSE
  16. # [10,] FALSE FALSE
  17. # [11,]  TRUE  TRUE
  18. # [12,]  TRUE  TRUE
  19. # [13,] FALSE FALSE
  20. # [14,] FALSE FALSE
  21. # [15,] FALSE FALSE
  22. # [16,] FALSE FALSE
  23. # [17,] FALSE FALSE
  24. # [18,] FALSE FALSE
  25. case_when(st_intersects(points_as_sf, combined_poly, sparse = F)[, 1] ~ "A",
  26. st_intersects(points_as_sf, combined_poly, sparse = F)[, 2] ~ "B",
  27. T ~ NA)
  28. # [1] "B" "B" "B" "B" "B" "B" "A" "A" NA  NA  "A" "A" NA  NA  NA  NA  NA  NA
英文:

In addition to the answer by @margusl building on the (default) sparse matrix form of sf::st_intersects() you can consider the alternative non-sparse output = a dense logical matrix.

While it is not as computationally efficient as sparse output I have found it easier to work with in many use cases.

Consider this piece of code; what it does is:

  • combines your polygons to a single sf object
  • casts your sp points to sf format
  • performs intersection of your two objects; the output is a logical matrix / points as rows, polygons as columns

The logical output produced by sf::st_intersects() is then passed to dplyr::case_when() which creates the variable - A if intersection with small polygon is true, B if intersection with big polygon is true and NA as a fallback default in all other cases.

This output can be easily passed to a vector variable or a new column in the points data frame.

  1. # two polygons to one sf object
  2. combined_poly &lt;- st_sf(geometry = c(pA, pB), row.names = c(&quot;pA&quot;, &quot;pB&quot;)) 
  3. # points as sf object
  4. points_as_sf &lt;- st_as_sf(sp)
  5. st_intersects(points_as_sf, combined_poly, sparse = F)
  6. #       [,1]  [,2]
  7. #  [1,] FALSE  TRUE
  8. #  [2,] FALSE  TRUE
  9. #  [3,] FALSE  TRUE
  10. #  [4,] FALSE  TRUE
  11. #  [5,] FALSE  TRUE
  12. #  [6,] FALSE  TRUE
  13. #  [7,]  TRUE  TRUE
  14. #  [8,]  TRUE  TRUE
  15. #  [9,] FALSE FALSE
  16. # [10,] FALSE FALSE
  17. # [11,]  TRUE  TRUE
  18. # [12,]  TRUE  TRUE
  19. # [13,] FALSE FALSE
  20. # [14,] FALSE FALSE
  21. # [15,] FALSE FALSE
  22. # [16,] FALSE FALSE
  23. # [17,] FALSE FALSE
  24. # [18,] FALSE FALSE
  25. case_when(st_intersects(points_as_sf, combined_poly, sparse = F)[, 1] ~ &quot;A&quot;,
  26. st_intersects(points_as_sf, combined_poly, sparse = F)[, 2] ~ &quot;B&quot;,
  27. T ~ NA)
  28. # [1] &quot;B&quot; &quot;B&quot; &quot;B&quot; &quot;B&quot; &quot;B&quot; &quot;B&quot; &quot;A&quot; &quot;A&quot; NA  NA  &quot;A&quot; &quot;A&quot; NA  NA  NA  NA  NA  NA

huangapple
  • 本文由 发表于 2023年5月26日 01:12:32
  • 转载请务必保留本文链接：https://go.coder-hub.com/76334776.html
  • geospatial
  • r
  • r-sf
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定