R – Calculating Sum of Line Lengths in R

Question

I have object1 with polygons and object2 with lines. I am looking for a way in R (preferably with the package sf) to calculate the sum of line lengths from object2 for each polygon in object1.

The resulting object should have the same features as object1, but with as additional attribute the length of the lines from object2 across each polygon of object1.

Answer 1

Read some data:

> polys = st_read("./lines_poly.shp", quiet=TRUE)
> lines = st_read("./lines_line.shp", quiet=TRUE)
> plot(polys$geom)
> plot(lines$geom,lwd=2,add=TRUE)

Gives:

Key operation is intersection:

> ints = st_intersection(lines, polys)

Which returns this:

> ints
Simple feature collection with 5 features and 2 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: -3.435776 ymin: 54.40466 xmax: -2.731861 ymax: 54.59513
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
    ID ID.1                       geometry
1    0    0 MULTILINESTRING ((-3.213695...
5    4    0 MULTILINESTRING ((-3.296152...
2    1    1 MULTILINESTRING ((-2.792331...
4    3    1 LINESTRING (-2.733722 54.50...
5.1  4    1 LINESTRING (-2.940045 54.53...

which is all the intersections of the lines with the polygons and ID attributes from the source lines and polygons. Map:

> plot(polys$geom)
> plot(ints$geom,lwd=2,add=TRUE)

So to get the length in each polygon, sum up the line segment lengths by polygon ID, which is ID.0:

> tapply(st_length(ints), ints$ID.1,sum)
       0        1 
61739.89 33265.18 

Polygons that don't have any lines in them are missing from this list, so you can't simply put them in the vector. Also, they might be in a random order...

To put these values in the right place in the polygon object a bit of match and indexing:

> polys$ltot = rep(0,nrow(polys))
> polys$ltot[match(names(ltot),polys$ID)] = ltot


> polys
Simple feature collection with 3 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -3.447643 ymin: 54.30678 xmax: -2.729412 ymax: 54.61172
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  ID                       geometry     ltot
1  0 POLYGON ((-3.428417 54.4994... 61739.89
2  1 POLYGON ((-2.940898 54.4747... 33265.18
3  2 POLYGON ((-3.238902 54.3676...     0.00

Then you can do:

> plot(polys[,"ltot"])