使用 Python 对网格内的总线长求和

Question

1。简介

一个折线文件，以.shp格式表示某处的路线图。
我想用栅格网络对整个区域进行栅格化。
每个网格应将其内的道路长度相加。

这就是我追求的。（此图是this question中类似项目的结果）

2。我的尝试

2.1 读取折线

目前，我正在考虑使用 shapely 包来应对。

road_file =road = map.readshapefile('roadmap','roadmap',zorder =1,)

路线图显示

2.2 拆分成多行

# x1,y1 --> start point; y1,y2 --> 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])

结果：

2.3 生成网格网络

格子网会遍布整个区域，我会把每个格子都变成一个身材匀称的侍从。

2.4 对网格内每条道路的长度求和

这一步我没有写代码。这是一个简单的例子，一个网格和一条路（实际上是线）。

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)

这样我就可以得到一个格子与一条线的交点长度。

但是当我的数据很大时（100x100 网格和 200 条道路），循环会很慢。必须有一些替代方法来处理这个问题。

如有任何建议，我们将不胜感激！

使用 rpy2 更新替代方法

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 <- readOGR(dsn = ".", layer ="roadmap")
roads_utm <- spTransform(roads, CRS("+init=epsg:21037"))
roads_utm_rst <- raster(extent(roads_utm), ncol = 100,nrow =100,crs = projection(roads_utm))
lengths <- sapply(1:ncell(roads_utm_rst), function(i) {
  tmp_rst <- roads_utm_rst
  tmp_rst[i] <- 1
  tmp_shp <- rasterToPolygons(tmp_rst)

  if (gIntersects(roads_utm, tmp_shp)) {
    roads_utm_crp <- crop(roads_utm, tmp_shp)
    roads_utm_crp_length <- gLength(roads_utm_crp)
    return(roads_utm_crp_length)
  } else {
    return(0)
  }
})
"""

length = r(rscripts)

那么，length是一个列表，里面包含了每一个格子的长度信息，那我就可以直接用python来分析了。

Answer 1

您可以尝试使用 shapely 的 MultiLineString 对象。

然后，您可以只对每个网格单元执行一次交集，而不是让 LineString 列表进行迭代。

roadLengthInCell = roadsMultiLineString.intersection(cellRect).length

使用 Python 对网格内的总线长求和

Sum the total line length inside grid using Python

python

gis

polygon

polyline

shapely

1。简介

2。我的尝试

使用 rpy2 更新替代方法