Using Python to convert line shapefile to raster, value=total length of lines within cell

Question

I have already read the similar question Convert line shapefile to raster, value=total length of lines within cell using R.
And I want to realize the function using Python.

===== My attempt ====

1. generate the grid network

    ### (1) 140 rows x 187 columns (1km x 1km grid)
    lon_x = np.linspace(xc1,xc2,187)
    lat_y = np.linspace(yc1,yc2,140)
   
    ### (2) empty array for saving the length data within the grid
    sh = (140*187,2)
    grids = np.zeros(140*187*2).reshape(*sh)
    grids_road = np.zeros(140*187)
   
    ### (3) generate the networks and k-d tree structure
    xx = lon_x
    yy = lat_y
   
    k = 0
    for j in range(0,yy.shape[0],1):
        for i in range(0,xx.shape[0],1):
            grids[k] = np.array([xx[i],yy[j]])
            k+=1
   
   ####  the k-d tree algorithm I used can be applied for finding the start/end grid agree with one road
   T = spatial.KDTree(grids)
   

2. Loop all my road polyline => intersection with grid => saving the length

   ### define the search radius
   x_delta = (lon_x[24] - lon_x[23])
   y_delta = (lat_y[24] - lat_y[23])
   R = np.sqrt(x_delta**2 + y_delta**2)
   
   ### loop the road(I save the polyline as pandas dataframe which contain the start point x1,y1 and end point x2, y2 and length) 
   
   for i in range(0,len(road),1):              
   
       sta = sorted(T.query_ball_point([road.x1.iloc[i],road.y1.iloc[i]],r=R))[0]
       end = sorted(T.query_ball_point([road.x2.iloc[i],road.y2.iloc[i]],r=R))[0]
   
       dt = (grids[end][0] - grids[sta][0])/x_delta
       dx = round((grids[end][0] - grids[sta][0])/x_delta)  
       dy = round((grids[end][1] - grids[sta][1])/y_delta)  
   
       # using shapely to change the road into line shapefile
       line = [(road.x1.iloc[i], road.y1.iloc[i]), (road.x2.iloc[i],road.y2.iloc[i])]
       shapely_line = shapely.geometry.LineString(line)
   
       # if the road extend just in one grid, there'll be no intersection
       if road.distance.iloc[i] < 1.0:
            index = sta
            grids_road[index] = road.distance.iloc[i]
       if road.distance.iloc[i] >= 1.0:
          if (dx > 0) & (dy == 0):
             for j in range(0, int(dx),1):
                 k = 0
                 x1,x2,x3,x4 = grids[sta][0] + j*x_delta, grids[sta][0] + (j+1)*x_delta,\
                          grids[sta][0] + (j+1)*x_delta, grids[sta][0] + j*x_delta
                 y1,y2,y3,y4 = grids[sta][1] + k*y_delta, grids[sta][1] + k*y_delta,\
                          grids[sta][1] + (k+1)*y_delta, grids[sta][1] + (k+1)*y_delta,
                 polygon = [(x1,y1),(x2,y2), (x3,y3), (x4,y4), (x1,y1)]
                 shapely_poly = shapely.geometry.Polygon(polygon)
                 intersection = list(shapely_poly.intersection(shapely_line).coords)
                 ds1,ds2 = intersection[0],intersection[1]
                 index = sta + j   
                 grids_road[index] = vincenty(ds1,ds2).miles
   

Then I classify the road into different scene based on the trends(dx, dy can represent this property). But I found this solution is too rigid, and I want to achieve the target with more efficient code.

Wish for your guide!

The figure shows my early progress which contain some vertical and horizontal road.
http://i4.tietuku.com/028bdd8ef1207557.png

Add

This method was too primary comparing with the method using R.
And, when the road wasn't parallel to Lon or Lat, the code using this method was too hard to write for different situations. So, I was think whether python has some similar functions or not?

Answer 1

Maybe you will find this post useful: http://www.machinalis.com/blog/python-for-geospatial-data-processing/