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 ====
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,1): for i in range(0,xx.shape,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)
Loop all my road polyline => intersection with grid => saving the length
### define the search radius x_delta = (lon_x - lon_x) y_delta = (lat_y - lat_y) 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)) end = sorted(T.query_ball_point([road.x2.iloc[i],road.y2.iloc[i]],r=R)) dt = (grids[end] - grids[sta])/x_delta dx = round((grids[end] - grids[sta])/x_delta) dy = round((grids[end] - grids[sta])/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] + j*x_delta, grids[sta] + (j+1)*x_delta,\ grids[sta] + (j+1)*x_delta, grids[sta] + j*x_delta y1,y2,y3,y4 = grids[sta] + k*y_delta, grids[sta] + k*y_delta,\ grids[sta] + (k+1)*y_delta, grids[sta] + (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,intersection index = sta + j grids_road[index] = vincenty(ds1,ds2).miles
Wish for your guide!
The figure shows my early progress which contain some vertical and horizontal road.
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?