# Intersecting lines and polygons and calculating line length in R?

I have shapefiles of polygons and lines at my disposal. I would like to intersect polygons and lines with R and thus calculate the length of the lines present in each polygon. You can find the data here "https://www.dropbox.com/s/b6x9gk8bzhan7zm/test.rar?dl=0" .I know how to do it on Arcgis but I would like to have a code on R to automate this because I have to do the same procedure for a lot of data. The code I put in place on R :

``````setwd("D:/folder/test")

library(sp)
library(maptools)
library(rgeos)
library(rgdal)

plot(Poly)

plot(line)

Inter <- gIntersection(line,Poly,byid=TRUE)
``````

allows me to have an intersection but what I would like to know is how to calculate then the length of the lines present in each polygon ?

I would recommend to use the sf package. You should also consider to use a projected CRS, in your case one, that is true to distance.

``````library("sf")
library("dplyr")

# intersection
int = st_intersection(line, poly)
# find out about the length of each line segment
int\$len = st_length(int)
# use a meaningful id (so far consists only of 0s)
poly\$Id = 1:nrow(poly)
# spatial overlay
join = st_join(poly, int)
# use the ID of the polygon for the aggregation
out = group_by(join, Id.x) %>%
summarize(length = sum(len))
# find out about polygons without line segments
filter(out, is.na(length))
# you can set the length of the polygons without line intersections to 0
# if you want
mutate(out, length = ifelse(is.na(length), 0, length))
``````

Of course, you can do the same in sp but sf already began to replace sp. To find out more about how to use sf, look at its vignettes (`vignette(package = "sf")`) and at our book Geocomputation with R.

• When we make the intersection between "line" and "poly", the polygons containing no line within them are not taken into account. At the end of the operation I would like to have also the (id) of polygons having no line in them but with 0 as length. – Falltok Apr 30 '18 at 16:17
• I have updated the code in accordance with your question but the answer has been already there, hence I am not entirely sure if this is what you meant... – Jannes May 3 '18 at 7:36
• I used the book you showed me above to export the file through : st_write(out,dsn="D:/Outa.shp", driver = "ESRI Shapefile") ,but the exported file have a different CRS than the input files. I used out<-st_transform(out, "+init=epsg:2154") to transform the file before exporting it but I always get a different CRS output. – Falltok May 6 '18 at 17:06
• Not sure what you exactly mean, please provide reproducible examples. `st_write` writes a .prj file containing the expected proj4-string. The shapefile format does not support EPSG codes (see also github.com/r-spatial/sf/issues/380). – Jannes May 7 '18 at 17:54