# Calculate sum of line lengths in R

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`.

``````> polys = st_read("./lines_poly.shp", quiet=TRUE)
> plot(polys\$geom)
``````

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

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"])
``````

• Any idea for SpatialLinesDataFrame objetcs? May 3, 2022 at 3:55
• @gonzalez.ivan90 convert to `sf` using `sf::st_as_sf` then proceed as above... May 3, 2022 at 6:53
• Thanks, yes. I used the st_as_sf( ) function! May 3, 2022 at 16:28