计算R中不规则形状的面积
Calculating the area of irregular shape in R
我想计算 R 中的封闭区域，如下图所示：
密码是：
a <- c(0,1,2,2,1,0,-1,-2,-2,-1,0)
b <- c(0,0,0,1,0,0,0,0,1,0,0)
id <- order(a)
AUC <- sum(abs(diff(a[id])*rollmean(b[id],2)))
AUC 的结果是 0.5。为什么不是 1？
当我做同样的事情，但针对不同的向量时，第二个例子：
代码是：
a <- c(0,1,1,0,0,0,1,1,0,0,0)
b <- c(1,1,2,1,0,-1,-1,-2,-1,0,1)
id <- order(a)
AUC <- sum(abs(diff(a[id])*rollmean(b[id],2)))
结果是1，符合我的预期。
最后一个最重要的问题是：
我在数据框中有 2 列 x 和 y（100 行，2 列），当绘制 x 与 y 时，它们形成不规则形状（圆图）的闭合曲线，如下图所示：
我想使用 2 列 x 和 y 计算该形状的面积。有办法吗？
更新：
第三个例子（循环图）的数据是：
structure(c(-0.317857301365341, -0.27254852000745, -0.239116190750992,
-0.212617899743053, -0.188665601384051, -0.164652457277714, -0.516811363687365,
-0.704142203990645, -0.925833008271225, -1.17644987082148, -1.44698060870818,
-1.72501192276032, -1.99567842934217, -2.24329812421979, -2.45336219338653,
-2.61434899714043, -2.71894439436842, -2.76453987881648, -2.75307392975534,
-2.69034804478649, -2.58498022688992, -2.44717942836848, -2.28751639917794,
-2.11583498227994, -1.94040889261091, -1.76740411983085, -1.37034128566268,
-0.945317248194438, -0.603583991912112, -0.245876237842544, 0.12405281597095,
0.500705070183803, 1.07928003293753, 1.3786192588209, 1.93625569819096,
2.32102987210859, 2.67667485884349, 2.9817217810439, 3.21508203608072,
3.35805065942255, 3.39629843099713, 3.32158207702952, 3.13297595252244,
2.83749468084523, 2.44998880428773, 1.99219777869042, 1.34133852767483,
0.905097245217795, 0.332495432968877, 0.000565676359279427, -0.304679355313777,
-0.676576930379378, -0.799504493858939, -0.877037748519715, -0.907917643885644,
-0.895574022744618, -0.847333419315191, -0.773304922267926, -0.68506594358595,
-0.594261154658402, -0.511242949335073, -0.443908244418344, -0.396899995824444,
-0.371311027504522, -2.95318107124334, -2.91494448796198, -2.80017869044755,
-2.61894793180667, -2.38576761259564, -2.11805651326137, -1.41368736155233,
-0.896608829738209, -0.363137342157785, 0.165324041578864, 0.668413502524084,
1.12765173596156, 1.52716559439287, 1.85443533430032, 2.10104197466128,
2.26318484127474, 2.34177274880676, 2.34207035998987, 2.2729829455437,
2.14607786936849, 1.97445291610407, 1.77157880225516, 1.5502403566333,
1.32167995293426, 1.09501806084454, 0.876985405337858, 0.780951474855999,
0.715653776169357, 0.631076472349792, 0.572023400386015, 0.534441560957589,
0.513328184556074, 0.275997813967411, 0.34478353180991, 0.504249489908016,
0.587434743269948, 0.654235885315243, 0.701450060739606, 0.72746455454756,
0.732567090012129, 0.718940924224628, 0.690264488986728, 0.650927827300914,
0.604983707715769, 0.555050419060634, 0.501440460475964, 0.370959009632643,
0.337726723524854, 0.216270960220213, 0.162458598954556, 0.060364386887637,
-0.136702597984821, -0.309488796931729, -0.535968167790145, -0.807754519418657,
-1.11285625385781, -1.43643730574278, -1.76172554792592, -2.0710195962479,
-2.34678141009815, -2.57280751486803, -2.73543188753935, -2.82464965255479,
-2.83501170706251), .Dim = c(64L, 2L), .Dimnames = list(c("309",
"310", "311", "312", "313", "314", "315", "316", "317", "318",
"319", "320", "321", "322", "323", "324", "325", "326", "327",
"328", "329", "330", "331", "332", "333", "334", "335", "336",
"337", "338", "339", "340", "341", "342", "343", "344", "345",
"346", "347", "348", "349", "350", "351", "352", "353", "354",
"355", "356", "357", "358", "359", "360", "361", "362", "363",
"364", "365", "366", "367", "368", "369", "370", "371", "372"
), c("PC1", "PC2")))
不确定您使用的方法是什么以及为什么在这种情况下不起作用。
另一种方法是使用 sf (simple features)，它处理面积计算：
library(sf)
a <- c(0,1,2,2,1,0,-1,-2,-2,-1,0)
b <- c(0,0,0,1,0,0,0,0,1,0,0)
poly <- st_polygon(list(cbind(a,b)))
poly <- st_make_valid(poly)
plot(poly)
st_area(poly)
#> [1] 1
请注意 st_make_valid 的使用，因为 多边形应该是不相交的 ，在您提供的三个示例中并非如此。
在第一个示例中，多边形转换为 2 个多边形和 2 个线串：
st_make_valid(poly)
GEOMETRYCOLLECTION (MULTIPOLYGON (((-2 1, -1 0, -2 0, -2 1)), ((2 0, 1 0, 2 1, 2 0))),
MULTILINESTRING ((0 0, -1 0), (0 0, 1 0)))
对于圆图，如果你想要一个精确的面积，就需要轨迹的交点，而你提供的数据并不是这样。
作为近似值，您可以使用 concaveman:
library(concaveman)
st_polygon(list(concaveman(data)))
plot(st_polygon(list(concaveman(data))), col='grey')
plot(st_linestring(data),add=T,col='red')
st_area(st_polygon(list(concaveman(data))))
[1] 2.693168
我想计算 R 中的封闭区域，如下图所示：
密码是：
a <- c(0,1,2,2,1,0,-1,-2,-2,-1,0)
b <- c(0,0,0,1,0,0,0,0,1,0,0)
id <- order(a)
AUC <- sum(abs(diff(a[id])*rollmean(b[id],2)))
AUC 的结果是 0.5。为什么不是 1？
当我做同样的事情，但针对不同的向量时，第二个例子：
a <- c(0,1,1,0,0,0,1,1,0,0,0)
b <- c(1,1,2,1,0,-1,-1,-2,-1,0,1)
id <- order(a)
AUC <- sum(abs(diff(a[id])*rollmean(b[id],2)))
结果是1，符合我的预期。
最后一个最重要的问题是：
我在数据框中有 2 列 x 和 y（100 行，2 列），当绘制 x 与 y 时，它们形成不规则形状（圆图）的闭合曲线，如下图所示：
更新： 第三个例子（循环图）的数据是：
structure(c(-0.317857301365341, -0.27254852000745, -0.239116190750992,
-0.212617899743053, -0.188665601384051, -0.164652457277714, -0.516811363687365,
-0.704142203990645, -0.925833008271225, -1.17644987082148, -1.44698060870818,
-1.72501192276032, -1.99567842934217, -2.24329812421979, -2.45336219338653,
-2.61434899714043, -2.71894439436842, -2.76453987881648, -2.75307392975534,
-2.69034804478649, -2.58498022688992, -2.44717942836848, -2.28751639917794,
-2.11583498227994, -1.94040889261091, -1.76740411983085, -1.37034128566268,
-0.945317248194438, -0.603583991912112, -0.245876237842544, 0.12405281597095,
0.500705070183803, 1.07928003293753, 1.3786192588209, 1.93625569819096,
2.32102987210859, 2.67667485884349, 2.9817217810439, 3.21508203608072,
3.35805065942255, 3.39629843099713, 3.32158207702952, 3.13297595252244,
2.83749468084523, 2.44998880428773, 1.99219777869042, 1.34133852767483,
0.905097245217795, 0.332495432968877, 0.000565676359279427, -0.304679355313777,
-0.676576930379378, -0.799504493858939, -0.877037748519715, -0.907917643885644,
-0.895574022744618, -0.847333419315191, -0.773304922267926, -0.68506594358595,
-0.594261154658402, -0.511242949335073, -0.443908244418344, -0.396899995824444,
-0.371311027504522, -2.95318107124334, -2.91494448796198, -2.80017869044755,
-2.61894793180667, -2.38576761259564, -2.11805651326137, -1.41368736155233,
-0.896608829738209, -0.363137342157785, 0.165324041578864, 0.668413502524084,
1.12765173596156, 1.52716559439287, 1.85443533430032, 2.10104197466128,
2.26318484127474, 2.34177274880676, 2.34207035998987, 2.2729829455437,
2.14607786936849, 1.97445291610407, 1.77157880225516, 1.5502403566333,
1.32167995293426, 1.09501806084454, 0.876985405337858, 0.780951474855999,
0.715653776169357, 0.631076472349792, 0.572023400386015, 0.534441560957589,
0.513328184556074, 0.275997813967411, 0.34478353180991, 0.504249489908016,
0.587434743269948, 0.654235885315243, 0.701450060739606, 0.72746455454756,
0.732567090012129, 0.718940924224628, 0.690264488986728, 0.650927827300914,
0.604983707715769, 0.555050419060634, 0.501440460475964, 0.370959009632643,
0.337726723524854, 0.216270960220213, 0.162458598954556, 0.060364386887637,
-0.136702597984821, -0.309488796931729, -0.535968167790145, -0.807754519418657,
-1.11285625385781, -1.43643730574278, -1.76172554792592, -2.0710195962479,
-2.34678141009815, -2.57280751486803, -2.73543188753935, -2.82464965255479,
-2.83501170706251), .Dim = c(64L, 2L), .Dimnames = list(c("309",
"310", "311", "312", "313", "314", "315", "316", "317", "318",
"319", "320", "321", "322", "323", "324", "325", "326", "327",
"328", "329", "330", "331", "332", "333", "334", "335", "336",
"337", "338", "339", "340", "341", "342", "343", "344", "345",
"346", "347", "348", "349", "350", "351", "352", "353", "354",
"355", "356", "357", "358", "359", "360", "361", "362", "363",
"364", "365", "366", "367", "368", "369", "370", "371", "372"
), c("PC1", "PC2")))
不确定您使用的方法是什么以及为什么在这种情况下不起作用。
另一种方法是使用 sf (simple features)，它处理面积计算：
library(sf)
a <- c(0,1,2,2,1,0,-1,-2,-2,-1,0)
b <- c(0,0,0,1,0,0,0,0,1,0,0)
poly <- st_polygon(list(cbind(a,b)))
poly <- st_make_valid(poly)
plot(poly)
st_area(poly)
#> [1] 1
请注意 st_make_valid 的使用，因为 多边形应该是不相交的 ，在您提供的三个示例中并非如此。
在第一个示例中，多边形转换为 2 个多边形和 2 个线串：
st_make_valid(poly)
GEOMETRYCOLLECTION (MULTIPOLYGON (((-2 1, -1 0, -2 0, -2 1)), ((2 0, 1 0, 2 1, 2 0))),
MULTILINESTRING ((0 0, -1 0), (0 0, 1 0)))
对于圆图，如果你想要一个精确的面积，就需要轨迹的交点，而你提供的数据并不是这样。
作为近似值，您可以使用 concaveman:
library(concaveman)
st_polygon(list(concaveman(data)))
plot(st_polygon(list(concaveman(data))), col='grey')
plot(st_linestring(data),add=T,col='red')
st_area(st_polygon(list(concaveman(data))))
[1] 2.693168