5

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.

1 Answer 1

8

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:

enter image description here

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)

enter image description here

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

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.