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]])
    ####  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.


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 1


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.