Search code examples
rareacurve

Calculating the area of irregular shape in R

I would like to calculate the closed area in R as in the graph below: example_1

The code is:

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)))

The result for AUC is 0.5. Why it is not 1?

When I do the same, but for different vectors, second example: enter image description here The code is:

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)))

The result is 1, which I have expected.

And the final, most important question is: I have 2 columns x and y in dataframe(100 rows, 2 columns) which form closed curve of irregular shape (cyclogram) when plotted x versus y, Like in this picture: enter image description here I want to calculate area of that shape using 2 columns x and y. Is there a way of doing this?

updated: Data for the 3rd example (cyclogram) is:

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")))
Solution

  • Not sure what is the method you used and why it didn't work in this case.
    An alternative is to use sf (simple features), which handles area calculation :

    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

    Note the use of st_make_valid because polygons should be non intersecting, which is not the case in the three examples you provided.
    In the first example, the polygon is transfomed in 2 polygons and 2 linestrings:

    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)))

    For the cyclogram, the intersection points of the trajectory are needed if you want an exact area, which is not the case in the data you provided.
    As an approximation, you could use 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

    enter image description here