2

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)

Poly <- readOGR("grille.shp",layer="grille")
plot(Poly)

line <- readOGR("line.shp",layer="line")
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 ?

3

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

poly = read_sf("grille.shp")
line = read_sf("line.shp")
# 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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