More

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

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">


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


Watch the video: Clip Raster in ArcMap Basic processing in GIS