Calculating area between two SpatialLines using R

Question

How do I calculate the area between two lines using R?

The two lines have the same start and end location, but the number of vertices for each line can vary.

I've tried to create a polygon from the coordinates of the vertices from the two line and then calculate the area by using rgeos::gArea but the resultant SpatialPolygons (sps) is incorrect.

x1 <- c(1,5,4,50)
y1 <- c(1,3,4,50)
line1 <- sp::SpatialLines(list(sp::Lines(sp::Line(cbind(x1,y1)), ID='a')))
x2 <- c(1,5,5,50)
y2 <- c(1,4,6,50)
line2 <- sp::SpatialLines(list(sp::Lines(sp::Line(cbind(x2,y2)), ID='b')))

plot(line1)
plot(line2, add = T, col = "red")

line_pts <- methods::as(line1, "SpatialPoints")
comparison_pts <- methods::as(line2, "SpatialPoints")

line_coords <- sp::coordinates(line_pts)
comparison_coords <- sp::coordinates(comparison_pts)

start <- line_coords[1, ]
end <- line_coords[base::nrow(line_coords), ]

coords <- base::rbind(line_coords, comparison_coords)

p = sp::Polygon(coords)
ps = sp::Polygons(list(p), 1)
sps = sp::SpatialPolygons(list(ps))

plot(sps, add = T, col  ="black") 

Answer 1

You need to reverse one of your lines to make a correctly oriented polygon.

If you modify your coords matrix by reversing the coordinates of the second line:

coords <- base::rbind(line_coords, comparison_coords[4:1,])

and then proceed unchanged:

p = sp::Polygon(coords)
ps = sp::Polygons(list(p), 1)
sps = sp::SpatialPolygons(list(ps))

You get

plot(sps, add = T, col  ="black") 

At this point my concern is that your polygon is twisted, and this can result in area calculations being wrong since twisted parts of polygons can be seen to have negative areas. Using the rgeos package we can compute the apparent area:

> library(rgeos)
> rgeos::gArea(sps)
[1] 23.5

Testing for validity is always a good idea:

> gIsValid(sps)
[1] FALSE
Warning message:
In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
  Self-intersection at or near point 4.4285714285714288 3.5714285714285716

so we should compare this with the polygon reconstructed from the lines that create it, this should undo the self intersection:

> spp2 = gPolygonize(gNode(gBoundary(sps)))
> plot(spp2)

the polygon looks the same but:

> gArea(spp2)
[1] 24.92857

it has grown. If you want to test this is correct, construct an example where you can easily see the overlap and compute by hand the correct area of the polygon.