网络分析(ArcPy)

一.前言

GIS中的网络分析最重要的便是纠正拓扑关系，建立矫正好的网络数据集，再进行网络分析，一般大家都是鼠标在arcgis上点点点，今天说一下Arcpy来解决的方案，对python的要求并不高,具体api参数查询arcgis帮助文档即可。

二.数据资源

在我的资源发布里，下载即可

三.步骤

新建数据库和数据集，并将数据导入数据集中
建立拓扑—导入要素进拓扑—对拓扑添加规则—检查拓扑—将拓扑错误导出
修正拓扑错误，并将选择出的拓扑点导出
调用 arcpy.SelectLayerByAttribute_management 按属性查询点，找到 NERA_DIST 数值相同的两个点，调用 Points To Line 工具将悬挂点连接起来，生成一个新的线图层。
将所有以点生成的线图层与原道路数据进行合并处理
删除所有以点生成的线图层
对拓扑处理后的道路添加字段计算时间成本
基于拓扑修正后的道路网创建网络数据集
创建服务区图层，加载设施点数据，解算设施点服务范围，并保存服务区域数据

四.代码

# --coding:utf-8--
import arcpy
import numpy as np

path = u"C:\\Users\\86152\\Desktop\\各科作业\\GIS算法\\workspace_5"
arcpy.env.workspace = path

arcpy.CreateFileGDB_management(path, "test5.gbd")  # 新建文件地理数据库
prj = u"C:\\Users\\86152\\Desktop\\各科作业\\GIS算法\\data_exp5\\CGCS2000 3 Degree GK CM 114E.prj"  # 坐标系参照

arcpy.CreateFeatureDataset_management("test5.gdb", "test5", prj)  # 新建要素数据集

# in_features = ['street1.shp', 'pharmacy.shp']  # 当数据就在当前工作空间时，可以直接用文件名
in_features = [u"C:\\Users\\86152\\Desktop\\各科作业\\GIS算法\\data_exp5\\street1.shp", u"C:\\Users\\86152\\Desktop"
               u"\\各科作业\\GIS算法\\data_exp5"
               u"\\pharmacy.shp"]
dataset_path = "C:\\Users\\86152\\Desktop\\各科作业\\GIS算法\\workspace_5\\test5.gdb\\test5"
arcpy.FeatureClassToGeodatabase_conversion(in_features, dataset_path)  # 将要素类导入至要素数据集中

arcpy.env.workspace = dataset_path  # 切换工作空间至要素数据集,便于访问
arcpy.CreateTopology_management(dataset_path, "street_Topology")  # 新建拓扑
topo_path = dataset_path + "\\street_Topology"
arcpy.AddFeatureClassToTopology_management(r"street_Topology", r"street1", 1, 1)  # 拓扑中添加要素类
arcpy.AddRuleToTopology_management(topo_path, "Must Not Have Dangles (Line)", "street1", "", "", "")  # 新增拓扑规则
arcpy.ValidateTopology_management(topo_path, "")  # 拓扑验证
arcpy.ExportTopologyErrors_management(topo_path, dataset_path, "F_topo")  # 输出拓扑错误

# 计算 F_topo_point 中每个点到其最近邻居点的距离，并将结果存储在 F_topo_point 的属性表中，生成几个新字段
arcpy.Near_analysis("F_topo_point", "F_topo_point")
# 近邻距离在0到20范围内的点可能是拓扑错误（例如，两个点非常接近但不相连，或者一个点悬挂在边上）
# 从 F_topo_point 中选择出 NEAR_DIST（近邻距离）在0到20之间的点，并将这些点保存到新的要素类 near_point 中
arcpy.Select_analysis("F_topo_point", "near_point", 'NEAR_DIST <20 and NEAR_DIST>0')
# 添加“NEAR”字段到near_point中
arcpy.AddField_management("near_point", "NEAR", "FLOAT", "", 6, "", "", "NULLABLE", "")
# 将NEAR_DIST的值赋值给NEAR字段
arcpy.CalculateField_management("near_point", "NEAR", '[NEAR_DIST]', "")  # 字段赋值最近距离

# 取出最近距离的相同的点
distance_List = []
shprows = arcpy.SearchCursor("near_point", ['NEAR'])
while True:
    shprow = shprows.next()
    if not shprow:
        break
    distance_List.append(shprow.NEAR)
distance_List = np.unique(distance_List)  # 去除重复值，获得唯一的最近距离列表（shp_List）

arcpy.MakeFeatureLayer_management("near_point", "near_point")  # 将要素转为图层

# 对每一个唯一的最近距离值，选择相应的点，并将这些点连接成线，生成新的线图层
line_List = []
i = 0
for p in distance_List:
    i = i + 1
    sql = '"NEAR"=@p'
    sql = sql.replace('@p', str(p))
    print(sql)
    ptl = 'point_to_line@i'
    ptl = ptl.replace('@i', str(i))
    arcpy.SelectLayerByAttribute_management("near_point", "NEW_SELECTION", sql)
    # 根据按属性选择出的点生成线存入ptl图层中
    arcpy.PointsToLine_management("near_point", ptl)
    # 将相应德线图层加入line_List
    line_List.append(ptl)

line_List.append("street1")
arcpy.Merge_management(line_List, "allstreets_temp")  # 将所有以点生成的线图层与原道路数据进行合并处理
arcpy.UnsplitLine_management("allstreets_temp", "allstreets")  # 合并具有重合端点的线

# 删除所有以点生成的线图层
i = 0
for p in distance_List:
    i = i + 1
    ptl = 'point_to_line@i'
    ptl = ptl.replace('@i', str(i))
    arcpy.Delete_management(ptl)

arcpy.AddField_management("allstreets", "Time", "FLOAT", "", 6, "", "", "NULLABLE", "")  # 对allstreets添加字段存储时间成本
arcpy.CalculateField_management("allstreets", "Time", '[Shape_Length]/80', "")  # 计算时间成本

print('接下来交给你了，请你完成网络分析！')

#########此处手动建立网络数据集############

# 根据allstreets创造出相应的网络数据集,完成之后解开下面的代码，注释掉上面的代码，运行即可解决


# dataset_path = "C:\\Users\\86152\\Desktop\\各科作业\\GIS算法\\workspace_5\\test5.gdb\\test5"
# arcpy.env.workspace = dataset_path
# outSAResultObject = arcpy.na.MakeServiceAreaLayer("test5_ND", "Network_analysis", "长度", "", 1000)  # 创建服务区图层
# outNALayer = outSAResultObject.getOutput(0)
# subLayerNames = arcpy.na.GetNAClassNames(outNALayer)
# facilitiesLayerName = subLayerNames["Facilities"]
# polygonsLayerName = subLayerNames["SAPolygons"]
# PolygonsSubLayer = arcpy.mapping.ListLayers(outNALayer, polygonsLayerName)[0]
# arcpy.na.AddLocations(outNALayer, facilitiesLayerName, "pharmacy", "", "")  # 加载pharmacy作为设施点
# arcpy.na.Solve(outNALayer)  # 求解
# arcpy.management.CopyFeatures(PolygonsSubLayer, "result")  # 将求解结果输出为要素