用 R 中的 gArea 和 gIntersection 错误计算多边形相交区域

如何解决用 R 中的 gArea 和 gIntersection 错误计算多边形相交区域

我有两个多边形

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[,coords1.sys[,col=alpha("black",0.5))

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

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

kt.small.s=data.frame(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[] 中的名称更改为更小的字符串或不带“.”的字符串。我仍然遇到同样的错误。

感谢您的帮助。

解决方法

不确定您追求的百分比是多少，但下面列出了所有领域供您进行数学计算。

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

如果您使用 sp 对象（SpatialPolygonsDataframes 等），通常可以使用 sf 将它们更改为 as(x,'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 日创建}