用R中的gArea和gIntersection误差计算多边形相交面积

Question

我有两个多边形

coords1[,1]
-15.99719 -16.01133 -16.04000 -16.08308 -16.14041 -16.21175 -16.29681 -16.39525 -16.50668 -16.63064 -16.76664 -16.91413 -17.07252 -17.24116 -17.41939 -17.60647 -17.80166 -18.00417 -18.21319 -18.42788 -18.64736 -18.87076 -19.09718 -19.32570 -19.55541 -19.78538 -20.01468 -20.24240 -20.46760 -20.68940 -20.90688 -21.11919 -21.32546 -21.52486 -21.71659 -21.89988 -22.07399 -22.23822 -22.39191 -22.53443 -22.66522 -22.78374 -22.88953 -22.98215 -23.06124 -23.12646 -23.17757 -23.21435 -23.23666 -23.24441 -23.23756 -23.21615 -23.18025 -23.13002 -23.06566 -22.98742 -22.89562 -22.79063 -22.67288 -22.54283 -22.40101 -22.24800 -22.08440 -21.91088 -21.72813 -21.53690 -21.33795 -21.13208-20.92012 -20.70293 -20.48137 -20.25635 -20.02876 -19.79953 -19.56958 -19.33983 -19.11120 -18.88463 -18.66102 -18.44127 -18.22626 -18.01687 -17.81393 -17.61826 -17.43066 -17.25187 -17.08262 -16.92358 -16.77540 -16.63868 -16.51396 -16.40174 -16.30249 -16.21659 -16.14440 -16.08620 -16.04224 -16.01268 -15.99764 -15.99719


coords1[,2]
6.225418 6.262789 6.300443 6.338229 6.375996 6.413591 6.450863 6.487662 6.523839 6.559250 6.593750 6.627202 6.659471 6.690426 6.719943 6.747904 6.774196 6.798712 6.821354 6.842031 6.860660 6.877166 6.891481 6.903550 6.913322 6.920758 6.925830 6.928515 6.928804 6.926696 6.922198 6.915329 6.906116 6.894597 6.880818 6.864834 6.846711 6.826520 6.804343 6.780269 6.754396 6.726828 6.697675 6.667055 6.635091 6.601912 6.567652 6.532448 6.496442 6.459780 6.422608 6.385077 6.347338 6.309542 6.271842 6.234389 6.197335 6.160829 6.125017 6.090044 6.056051 6.023174 5.991546 5.961295 5.932541 5.905401 5.879985 5.856394 5.834723 5.815060 5.797484 5.782066 5.768867 5.757941 5.749333 5.743075 5.739195 5.737707 5.738617 5.741922 5.747609 5.755654 5.766025 5.778680 5.793569 5.810631 5.829798 5.850993 5.874129 5.899115 5.925849 5.954224 5.984126 6.015434 6.048021 6.081758 6.116508 6.152130 6.188483 6.225418

coords1.sys[,1]
-17.27076 -17.28943 -17.32370 -17.37346 -17.43849 -17.51853 -17.61326 -17.72229 -17.84520 -17.98148 -18.13059 -18.29192 -18.46483 -18.64862 -18.84254 -19.04583 -19.25765 -19.47717 -19.70348 -19.93569 -20.17285 -20.41402 -20.65821 -20.90445 -21.15175 -21.39911 -21.64554 -21.89003 -22.13161 -22.36931 -22.60216 -22.82924 -23.04962 -23.26241 -23.46677 -23.66187 -23.84692 -24.02117 -24.18394 -24.33455 -24.47241 -24.59695 -24.70769 -24.80416 -24.88599 -24.95284 -25.00445 -25.04060 -25.06116 -25.06603 -25.05520 -25.02871 -24.98667 -24.92926 -24.85668 -24.76926 -24.66732 -24.55129 -24.42163 -24.27886 -24.12357 -23.95636 -23.77792 -23.58897 -23.39027 -23.18261 -22.96683 -22.74380 -22.51443 -22.27962 -22.04034 -21.79754 -21.55219 -21.30530 -21.05785 -20.81083 -20.56525 -20.32209 -20.08233 -19.84693 -19.61685 -19.39301 -19.17632 -18.96763 -18.76781 -18.57764 -18.39789 -18.22930 -18.07253 -17.92822
-17.79696 -17.67926 -17.57560 -17.48640 -17.41203 -17.35277 -17.30887 -17.28050 -17.26778 -17.27076

coords1.sys[,2]
6.465179 6.538004 6.611223 6.684541 6.757663 6.830295 6.902143 6.972919 7.042338 7.110120 7.175992 7.239689 7.300955 7.359542 7.415215 7.467750 7.516935 7.562573 7.604478 7.642483 7.676435 7.706197 7.731648 7.752688 7.769230 7.781208 7.788574 7.791299 7.789371 7.782798 7.771606 7.755841 7.735567 7.710864 7.681832 7.648589 7.611267 7.570018 7.525007 7.476416 7.424440 7.369288 7.311183 7.250359 7.187060 7.121542 7.054068 6.984910 6.914346 6.842661 6.770143 6.697084 6.623779 6.550522 6.477608 6.405332 6.333984 6.263851 6.195216 6.128355 6.063538 6.001025 5.941069 5.883910 5.829778 5.778892 5.731457 5.687664 5.647688 5.611692 5.579819 5.552199 5.528941 5.510141 5.495874 5.486197 5.481149 5.480751 5.485003 5.493890 5.507375 5.525404 5.547904 5.574784 5.605938 5.641238 5.680543 5.723695 5.770520 5.820829 5.874419 5.931076 5.990571 6.052663 6.117104 6.183634 6.251984 6.321880 6.393041 6.465179


plot(1,1,ylim=c(0,15),xlim=c(-30,-10), t="n", xlab="", ylab="")
polygon(coords1[,1], coords1[,2], border=2,col=alpha("red",0.5))
polygon(coords1.sys[,1], coords1.sys[,2], border=2,col=alpha("black",0.5))

我正在尝试计算重叠百分比，如使用 gArea 和 gIntersection here 所示。

kt.small.g=data.frame(coords1[,1], coords1[,2])
names(kt.small.g)=c("x","y")
summary(kt.small.g)

kt.small.s=data.frame(coords1.sys[,1], coords1.sys[,2])
names(kt.small.s)=c("x","y")
kt.small.g_pol <- Polygons(list(Polygon(kt.small.g)), "small.g")
kt.small.s_pol <- Polygons(list(Polygon(kt.small.s)), "small.s")

kt.small.g_area <- gArea(shape["small.g"])
kt.small.intersections <- gIntersection(shape["small.g"], shape["small.s"])

但是在运行kt.small.g_area <- gArea(shape["small.g"]) 或kt.small.intersections <- gIntersection(shape["small.g"], shape["small.s"]) 之后我得到了错误：

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'obj' in selecting a method for function 'is.projected': NAs not permitted in row index

我已经完成了链接示例中的代码，没有出现错误。我检查了我的 NA 数据，与示例并行结构，并尝试将 shape[] 中的名称更改为较小的字符串或不带“。”的字符串。我仍然得到同样的错误。

感谢您的帮助。

Answer 1

不确定您所追求的百分比，但下面列出了所有区域供您计算。

以下解决方案使用sf 库来处理空间对象。

如果您使用的是sp 对象（SpatialPolygonsDataframes 等），通常可以使用as(x, 'sf') 将它们更改为sf。

library(dplyr)
library(ggplot2)
library(sf)

# Cleaning your data pasted from SO 
poly1_x <- as.numeric(c("-15.99719", "-16.01133", "-16.04000", "-16.08308", "-16.14041", "-16.21175", "-16.29681", "-16.39525", "-16.50668", "-16.63064", "-16.76664", "-16.91413", "-17.07252", "-17.24116", "-17.41939", "-17.60647", "-17.80166", "-18.00417", "-18.21319", "-18.42788", "-18.64736", "-18.87076", "-19.09718", "-19.32570", "-19.55541", "-19.78538", "-20.01468", "-20.24240", "-20.46760", "-20.68940", "-20.90688", "-21.11919", "-21.32546", "-21.52486", "-21.71659", "-21.89988", "-22.07399", "-22.23822", "-22.39191", "-22.53443", "-22.66522", "-22.78374", "-22.88953", "-22.98215", "-23.06124", "-23.12646", "-23.17757", "-23.21435", "-23.23666", "-23.24441", "-23.23756", "-23.21615", "-23.18025", "-23.13002", "-23.06566", "-22.98742", "-22.89562", "-22.79063", "-22.67288", "-22.54283", "-22.40101", "-22.24800", "-22.08440", "-21.91088", "-21.72813", "-21.53690", "-21.33795", "-21.13208","-20.92012", "-20.70293", "-20.48137", "-20.25635", "-20.02876", "-19.79953", "-19.56958", "-19.33983", "-19.11120", "-18.88463", "-18.66102", "-18.44127", "-18.22626", "-18.01687", "-17.81393", "-17.61826", "-17.43066", "-17.25187", "-17.08262", "-16.92358", "-16.77540", "-16.63868", "-16.51396", "-16.40174", "-16.30249", "-16.21659", "-16.14440", "-16.08620", "-16.04224", "-16.01268", "-15.99764", "-15.99719"))
poly1_y <- c(6.225418, 6.262789, 6.300443, 6.338229, 6.375996, 6.413591, 6.450863, 6.487662, 6.523839, 6.55925, 6.59375, 6.627202, 6.659471, 6.690426, 6.719943, 6.747904, 6.774196, 6.798712, 6.821354, 6.842031, 6.86066, 6.877166, 6.891481, 6.90355, 6.913322, 6.920758, 6.92583, 6.928515, 6.928804, 6.926696, 6.922198, 6.915329, 6.906116, 6.894597, 6.880818, 6.864834, 6.846711, 6.82652, 6.804343, 6.780269, 6.754396, 6.726828, 6.697675, 6.667055, 6.635091, 6.601912, 6.567652, 6.532448, 6.496442, 6.45978, 6.422608, 6.385077, 6.347338, 6.309542, 6.271842, 6.234389, 6.197335, 6.160829, 6.125017, 6.090044, 6.056051, 6.023174, 5.991546, 5.961295, 5.932541, 5.905401, 5.879985, 5.856394, 5.834723, 5.81506, 5.797484, 5.782066, 5.768867, 5.757941, 5.749333, 5.743075, 5.739195, 5.737707, 5.738617, 5.741922, 5.747609, 5.755654, 5.766025, 5.77868, 5.793569, 5.810631, 5.829798, 5.850993, 5.874129, 5.899115, 5.925849, 5.954224, 5.984126, 6.015434, 6.048021, 6.081758, 6.116508, 6.15213, 6.188483, 6.225418)

poly2_x <- c("-17.27076 -17.28943 -17.32370 -17.37346 -17.43849 -17.51853 -17.61326 -17.72229 -17.84520 -17.98148 -18.13059 -18.29192 -18.46483 -18.64862 -18.84254 -19.04583 -19.25765 -19.47717 -19.70348 -19.93569 -20.17285 -20.41402 -20.65821 -20.90445 -21.15175 -21.39911 -21.64554 -21.89003 -22.13161 -22.36931 -22.60216 -22.82924 -23.04962 -23.26241 -23.46677 -23.66187 -23.84692 -24.02117 -24.18394 -24.33455 -24.47241 -24.59695 -24.70769 -24.80416 -24.88599 -24.95284 -25.00445 -25.04060 -25.06116 -25.06603 -25.05520 -25.02871 -24.98667 -24.92926 -24.85668 -24.76926 -24.66732 -24.55129 -24.42163 -24.27886 -24.12357 -23.95636 -23.77792 -23.58897 -23.39027 -23.18261 -22.96683 -22.74380 -22.51443 -22.27962 -22.04034 -21.79754 -21.55219 -21.30530 -21.05785 -20.81083 -20.56525 -20.32209 -20.08233 -19.84693 -19.61685 -19.39301 -19.17632 -18.96763 -18.76781 -18.57764 -18.39789 -18.22930 -18.07253 -17.92822 -17.79696 -17.67926 -17.57560 -17.48640 -17.41203 -17.35277 -17.30887 -17.28050 -17.26778 -17.27076")
poly2_x <- stringr::str_split(poly2_x, pattern = ' ') %>% unlist()
poly2_y <- c(6.465179, 6.538004, 6.611223, 6.684541, 6.757663, 6.830295, 6.902143, 6.972919, 7.042338, 7.11012, 7.175992, 7.239689, 7.300955, 7.359542, 7.415215, 7.46775, 7.516935, 7.562573, 7.604478, 7.642483, 7.676435, 7.706197, 7.731648, 7.752688, 7.76923, 7.781208, 7.788574, 7.791299, 7.789371, 7.782798, 7.771606, 7.755841, 7.735567, 7.710864, 7.681832, 7.648589, 7.611267, 7.570018, 7.525007, 7.476416, 7.42444, 7.369288, 7.311183, 7.250359, 7.18706, 7.121542, 7.054068, 6.98491, 6.914346, 6.842661, 6.770143, 6.697084, 6.623779, 6.550522, 6.477608, 6.405332, 6.333984, 6.263851, 6.195216, 6.128355, 6.063538, 6.001025, 5.941069, 5.88391, 5.829778, 5.778892, 5.731457, 5.687664, 5.647688, 5.611692, 5.579819, 5.552199, 5.528941, 5.510141, 5.495874, 5.486197, 5.481149, 5.480751, 5.485003, 5.49389, 5.507375, 5.525404, 5.547904, 5.574784, 5.605938, 5.641238, 5.680543, 5.723695, 5.77052, 5.820829, 5.874419, 5.931076, 5.990571, 6.052663, 6.117104, 6.183634, 6.251984, 6.32188, 6.393041, 6.465179)

# Make polygons as sf objects
poly1 <- data.frame(x = poly1_x, y = poly1_y) %>% 
  st_as_sf(coords = c('x', 'y')) %>%
  st_combine() %>%
  st_cast('POLYGON')

poly2 <- data.frame(x = poly2_x, y = poly2_y) %>% 
  st_as_sf(coords = c('x', 'y')) %>%
  st_combine() %>%
  st_cast('POLYGON')

# Plot
ggplot() + 
  geom_sf(data = poly1, fill = 'red', alpha = .5) + 
  geom_sf(data = poly2, fill = 'blue', alpha = .5)

# Area of polygon1: (red)
st_area(poly1)
#> [1] 6.62603

# Area of polygon2: (blue)
st_area(poly2)
#> [1] 13.88645

# Area covered only by BOTH polygons (purple)
st_intersection(poly1, poly2) %>%
  st_area
#> [1] 5.678321

^{由reprex package (v0.3.0) 于 2020 年 12 月 30 日创建}