Sum the total line length inside grid using Python

1. Introduction

• A polylines file which in .shp format represent the roadmap of somewhere.
• I want to raster the whole area with grid network.
• Each grid should sum the road length inside it.

This is what I pursued.(This figure is a result of similar project which shown in this question)

2. My Attempt

2.1 read the polylines

For now, I was thinking about using shapely package to cope with.

road_file =road = map.readshapefile('roadmap','roadmap',zorder =1,)

The roadmap shows like

2.2 Split it into multi-lines

# x1,y1 --> start point; y1,y2 --> end point
road = {"name":[],"x1":[],"x2":[],"y1":[],"y2":[],"distance":[]}
l = []
ll = []
i = 0
for info, shape in zip(roadmap.info, roadmap):
    i+=1
    road["name"].append("road" + str(i))
    l = info
    ll.append(info)
    road["x1"].append(shape[0][0])
    road["x2"].append(shape[-1][0])
    road["y1"].append(shape[0][1])
    road["y2"].append(shape[-1][1])

The result:

2.3 Generate the grid network

The grid network will be all over the area, and I'll transform each grid into a shapely squire.

2.4 Sum the length of each road inside of the grid

This step I didn't write in code. And here is an simple example for one grid with one road(line, indeed).

import shapely
polygon = [(1.0, 1.0), (2.0, 1.0), (2.0, 2.0), (1,2), (1, 1)]
shapely_poly = shapely.geometry.Polygon(polygon)
line = [(1, 2), (2,1)]
shapely_line = shapely.geometry.LineString(line)
intersection_line = list(shapely_poly.intersection(shapely_line).coords)

In this way, I can get the intersection length of one grid with one line.

But when my data was vast(100x100 grids & 200 roads), the loop would be so slow. There must be some alternative way to deal with this issue.

Any advices would be appreciate!

Update an alternative method using rpy2

import rpy2
from rpy2.robjects import r
from rpy2.robjects.packages import importr
rgdal  = importr('rgdal')
raster = importr("raster")
rgeos  = importr("rgeos") 

rscript = """
##Reference the question I mentioned above
roads <- readOGR(dsn = ".", layer ="roadmap")
roads_utm <- spTransform(roads, CRS("+init=epsg:21037"))
roads_utm_rst <- raster(extent(roads_utm), ncol = 100,nrow =100,crs = projection(roads_utm))
lengths <- sapply(1:ncell(roads_utm_rst), function(i) {
  tmp_rst <- roads_utm_rst
  tmp_rst[i] <- 1
  tmp_shp <- rasterToPolygons(tmp_rst)

  if (gIntersects(roads_utm, tmp_shp)) {
    roads_utm_crp <- crop(roads_utm, tmp_shp)
    roads_utm_crp_length <- gLength(roads_utm_crp)
    return(roads_utm_crp_length)
  } else {
    return(0)
  }
})
"""

length = r(rscripts)  

Then, length is a list contain length information of every grid, then I can use python directly to analyse it.