# 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)
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?

• I have thought about one solution: Dec 30, 2015 at 3:06

## 1 Answer

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