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

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)

### (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

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)

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
shapely_line = shapely.geometry.LineString(line)

# if the road extend just in one grid, there'll be no intersection
index = sta
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

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.

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