title: Sum the total line length inside grid using Python
body: <h3>1. Introduction</h3>
  
  <ul>
  <li>A polylines file which in <em>.shp</em> format represent the roadmap of somewhere.    </li>
  <li>I want to raster the whole area with grid network.   </li>
  <li>Each grid should sum the road length inside it.      </li>
  </ul>
  
  <p>This is what I pursued.(This figure is a result of similar project which shown in <a href="https://gis.stackexchange.com/questions/119993/convert-line-shapefile-to-raster-value-total-length-of-lines-within-cell">this question</a>)<br>
  <a href="https://i.stack.imgur.com/0RvoC.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/0RvoC.png" alt="enter image description here"></a></p>
  
  <h3>2. My Attempt</h3>
  
  <p>2.1 read the polylines</p>
  
  <p>For now, I was thinking about using <em>shapely</em> package to cope with.     </p>
  
  <pre><code>road_file =road = map.readshapefile('roadmap','roadmap',zorder =1,)
  </code></pre>
  
  <p>The roadmap shows like<br>
  <a href="https://i.stack.imgur.com/YMCME.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/YMCME.png" alt="enter image description here"></a></p>
  
  <p>2.2 Split it into multi-lines   </p>
  
  <pre><code># x1,y1 --&gt; start point; y1,y2 --&gt; end point
  road = {"name":[],"x1":[],"x2":[],"y1":[],"y2":[],"distance":[]}
  l = []
  ll = []
  i = 0
  for info, shape in zip(roadmap.info, roadmap):
      i+=1
      road["name"].append("road" + str(i))
      l = info
      ll.append(info)
      road["x1"].append(shape[0][0])
      road["x2"].append(shape[-1][0])
      road["y1"].append(shape[0][1])
      road["y2"].append(shape[-1][1])
  </code></pre>
  
  <p>The result:<br>
  <a href="https://i.stack.imgur.com/xY0I0.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/xY0I0.png" alt="enter image description here"></a>         </p>
  
  <p>2.3 Generate the grid network </p>
  
  <p>The grid network will be all over the area, and I'll transform each grid into a shapely squire. </p>
  
  <p><a href="https://i.stack.imgur.com/6P68b.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/6P68b.png" alt="enter image description here"></a></p>
  
  <p>2.4 Sum the length of each road inside of the grid   </p>
  
  <p>This step I didn't write in code. And here is an simple example for one grid with one road(line, indeed).    </p>
  
  <pre><code>import shapely
  polygon = [(1.0, 1.0), (2.0, 1.0), (2.0, 2.0), (1,2), (1, 1)]
  shapely_poly = shapely.geometry.Polygon(polygon)
  line = [(1, 2), (2,1)]
  shapely_line = shapely.geometry.LineString(line)
  intersection_line = list(shapely_poly.intersection(shapely_line).coords)
  </code></pre>
  
  <p>In this way, I can get the intersection length of one grid with one line.    </p>
  
  <p>But when my data was vast(100x100 grids &amp; 200 roads), the loop would be so slow. There must be some alternative way to deal with this issue.    </p>
  
  <p>Any advices would be appreciate!   </p>
  
  <h3>Update an alternative method using rpy2</h3>
  
  <pre><code>import rpy2
  from rpy2.robjects import r
  from rpy2.robjects.packages import importr
  rgdal  = importr('rgdal')
  raster = importr("raster")
  rgeos  = importr("rgeos") 
  
  rscript = """
  ##Reference the question I mentioned above
  roads &lt;- readOGR(dsn = ".", layer ="roadmap")
  roads_utm &lt;- spTransform(roads, CRS("+init=epsg:21037"))
  roads_utm_rst &lt;- raster(extent(roads_utm), ncol = 100,nrow =100,crs = projection(roads_utm))
  lengths &lt;- sapply(1:ncell(roads_utm_rst), function(i) {
    tmp_rst &lt;- roads_utm_rst
    tmp_rst[i] &lt;- 1
    tmp_shp &lt;- rasterToPolygons(tmp_rst)
  
    if (gIntersects(roads_utm, tmp_shp)) {
      roads_utm_crp &lt;- crop(roads_utm, tmp_shp)
      roads_utm_crp_length &lt;- gLength(roads_utm_crp)
      return(roads_utm_crp_length)
    } else {
      return(0)
    }
  })
  """
  
  length = r(rscripts)  
  </code></pre>
  
  <p>Then, <em>length</em> is a list contain length information of every grid, then I can use python directly to analyse it.</p>