## 设计思想

 1 2 3 4 5  import networkx as nx import math from itertools import tee import shapefile from osgeo import ogr 

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19  def pairwise(iterable): """返回可迭代访问的二值元组 s -> (s0,s1), (s1,s2), (s2, s3), ...""" a, b = tee(iterable) next(b, None) return zip(a, b) def haversine(n0, n1): x1, y1 = n0 x2, y2 = n1 x_dist = math.radians(x1 - x2) y_dist = math.radians(y1 - y2) y1_rad = math.radians(y1) y2_rad = math.radians(y2) a = math.sin(y_dist/2)**2 + math.sin(x_dist/2)**2 \ * math.cos(y1_rad) * math.cos(y2_rad) c = 2 * math.asin(math.sqrt(a)) distance = c * 6371 return distance 

 1 2 3 4 5 6  shp = "road_network.shp" G = nx.DiGraph() r = shapefile.Reader(shp) for s in r.shapes(): for p1, p2 in pairwise(s.points): G.add_edge(tuple(p1), tuple(p2)) 

 1  sg = list(nx.connected_component_subgraphs(G.to_undirected()))[0] 

 1  sg = list(G.to_undirected(c) for c in nx.strongly_connected_components (G))[0] 

 1 2 3  r = shapefile.Reader(r"start_end.shp") start = r.shape(0).points[0] end = r.shape(1).points[0] 

 1 2 3  for n0, n1 in sg.edges(): dist = haversine(n0, n1) sg.edges[n0,n1]["dist"] = dist 

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16  nn_start = None nn_end = None start_delta = float("inf") end_delta = float("inf") for n in sg.nodes(): s_dist = haversine(start, n) e_dist = haversine(end, n) if s_dist < start_delta: nn_start = n start_delta = s_dist if e_dist < end_delta: nn_end = n end_delta = e_dist path = nx.shortest_path(sg, source=nn_start, target=nn_end, weight="dist") #list 

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31  multiline = ogr.Geometry(ogr.wkbMultiLineString) # 构建几何类型:线 line = ogr.Geometry(ogr.wkbLineString) for point in path: line.AddPoint(point[0],point[1]) multiline.AddGeometry(line) wkt = multiline.ExportToWkt() driver = ogr.GetDriverByName("ESRI Shapefile") data_source = driver.CreateDataSource("shortest.shp") srs = osr.SpatialReference() srs.ImportFromEPSG(4326) layer = data_source.CreateLayer("shortest", srs, ogr.wkbLineString ) field_name = ogr.FieldDefn("Name", ogr.OFTString) field_name.SetWidth(14) layer.CreateField(field_name) field_name = ogr.FieldDefn("data", ogr.OFTString) field_name.SetWidth(14) layer.CreateField(field_name) feature = ogr.Feature(layer.GetLayerDefn()) feature.SetField("Name", "shortest") feature.SetField("data", "1.2") line = ogr.CreateGeometryFromWkt(wkt) feature.SetGeometry(line) layer.CreateFeature(feature) feature = None data_source = None 

## 后记

Author by Jerrychoices
Built with Hugo