import arcpy
import os
# 设置工作空间
gdb_path = r"F:\Project\data-2\Transportation.gdb"
arcpy.env.workspace = gdb_path
arcpy.env.overwriteOutput = True
# 定义要素数据集路径
feature_dataset = os.path.join(gdb_path, "Transportation")
# 定义数据源路径
datasets = {
"Roads": r"F:\Project\data-2\道路.SHP",
"GasStations": r"F:\Project\data-2\加油站.SHP",
"Communities": r"F:\Project\data-2\小区.SHP",
"Districts": r"F:\Project\data-2\区划范围.SHP"
}
# 获取源数据投影(使用道路数据)
prj = arcpy.Describe(datasets["Roads"]).spatialReference
# 创建要素数据集
if not arcpy.Exists(feature_dataset):
arcpy.CreateFeatureDataset_management(
out_dataset_path=gdb_path,
out_name="Transportation",
spatial_reference=prj
)
# 导入数据到要素数据集
for name, path in datasets.items():
out_fc = os.path.join(feature_dataset, name)
if not arcpy.Exists(out_fc):
arcpy.FeatureClassToFeatureClass_conversion(
in_features=path,
out_path=feature_dataset,
out_name=name
)
# 2. 道路数据处理(排除高速公路)
valid_roads = os.path.join(feature_dataset, "ValidRoads")
if not arcpy.Exists(valid_roads):
arcpy.analysis.Select(
in_features=os.path.join(feature_dataset, "Roads"),
out_feature_class=valid_roads,
where_clause="类型 NOT IN ('高速公路')"
)
# 验证所有必要要素类存在
required_fcs = ["Roads", "ValidRoads", "GasStations", "Communities", "Districts"]
for fc in required_fcs:
path = os.path.join(feature_dataset, fc)
if not arcpy.Exists(path):
arcpy.AddError(f"错误:缺失必要要素类 {fc}")
exit()
# 3. 创建拓扑
topology_name = "RoadTopology"
topology_path = os.path.join(feature_dataset, topology_name)
# 创建拓扑(确保拓扑不存在)
if not arcpy.Exists(topology_path):
arcpy.management.CreateTopology(
in_dataset=feature_dataset,
out_topology=topology_name,
cluster_tolerance=0.001
)
print(f"拓扑已创建: {topology_path}")
# 获取拓扑参与者的可靠方法
def get_topology_participants(topology_path):
"""获取拓扑中所有参与者的名称列表"""
participants = []
try:
# 方法1: 使用拓扑规则信息获取参与者
desc = arcpy.Describe(topology_path)
if hasattr(desc, "rules") and desc.rules:
for rule in desc.rules:
# 获取规则涉及的主要要素类
if hasattr(rule, "featureClassName") and rule.featureClassName:
if rule.featureClassName not in participants:
participants.append(rule.featureClassName)
# 获取规则涉及的次要要素类(某些规则需要两个要素类)
if hasattr(rule, "featureClass2Name") and rule.featureClass2Name:
if rule.featureClass2Name not in participants:
participants.append(rule.featureClass2Name)
# 方法2: 如果规则方法失败,尝试其他方法
if not participants:
try:
# 使用ListDatasets获取参与者
datasets = arcpy.ListDatasets(topology_path, "FeatureClass")
if datasets:
participants = list(datasets)
except:
pass
return participants
except Exception as e:
print(f"获取拓扑参与者时出错: {str(e)}")
return []
# 获取拓扑状态
topology_desc = arcpy.Describe(topology_path)
# 使用isValid检查拓扑是否有效
if topology_desc.isValid:
print(f"拓扑 {topology_name} 是有效的")
else:
print(f"拓扑 {topology_name} 无效")
# 添加要素类到拓扑
safe_add_featureclass_to_topology(topology_path, valid_roads)
# 添加拓扑规则 - 使用更可靠的方法
add_topology_rule(topology_path, "Must Not Have Dangles", valid_roads, rule_id=2)
# 验证拓扑
print("正在验证拓扑...")
arcpy.management.ValidateTopology(topology_path)
print("拓扑验证完成")
### 添加要素类到拓扑(如果没有参与者)
##if not topology_has_participants(topology_path):
## print("正在添加要素类到拓扑...")
## arcpy.management.AddFeatureClassToTopology(
## in_topology=topology_path,
## in_featureclass=valid_roads,
## xy_rank=1,
## z_rank=1
## )
## print(f"要素类 {valid_roads} 已添加到拓扑")
##else:
## print("拓扑中已有参与者,跳过添加")
# 添加拓扑规则 - 使用更可靠的方法
def add_topology_rule(topology_path, rule_type, feature_class, rule_id=None):
"""添加拓扑规则的健壮方法"""
try:
# 尝试使用规则名称添加
arcpy.management.AddRuleToTopology(
in_topology=topology_path,
rule_type=rule_type,
in_featureclass=feature_class,
subtype="",
in_featureclass2="",
subtype2=""
)
print(f"成功添加规则: {rule_type}")
return True
except arcpy.ExecuteError:
# 如果名称失败,尝试使用规则ID
if rule_id is not None:
try:
arcpy.management.AddRuleToTopology(
in_topology=topology_path,
rule_type=rule_id,
in_featureclass=feature_class,
subtype=""
)
print(f"使用ID {rule_id} 成功添加规则")
return True
except Exception as e:
print(f"使用ID添加规则失败: {str(e)}")
# 获取详细的错误信息
error_msgs = arcpy.GetMessages(2)
print(f"添加拓扑规则失败: {error_msgs}")
return False
# 添加"Must Not Have Dangles"规则
add_topology_rule(topology_path, "Must Not Have Dangles", valid_roads, rule_id=2)
# 验证拓扑
print("正在验证拓扑...")
arcpy.management.ValidateTopology(topology_path)
print("拓扑验证完成")
# 检查拓扑状态
topology_desc = arcpy.Describe(topology_path)
print(f"拓扑状态: {topology_desc.status}")
# ... (后续代码保持不变) ...
# 4. 创建网络数据集
network_dataset_name = "RoadNetwork"
nd_path = os.path.join(feature_dataset, network_dataset_name)
if not arcpy.Exists(nd_path):
try:
# 创建网络数据集
arcpy.na.CreateNetworkDataset(
in_feature_dataset=feature_dataset,
out_name=network_dataset_name,
template="Road"
)
# 获取网络数据集对象
nd_layer = "nd_layer"
arcpy.na.MakeNetworkDatasetLayer(nd_path, nd_layer)
# 添加道路源
arcpy.na.AddSourcesToNetworkDataset(
in_network_dataset=nd_layer,
sources_table=[[
valid_roads, # 要素类路径
"SHAPE", # 几何字段
"类型", # 道路类型字段
"Shape_Leng", # 长度字段
"", "", "", "",
"Roads" # 源名称
]]
)
# 添加网络属性:成本字段
arcpy.na.AddNetworkAttributes(
in_network_dataset=nd_layer,
attributes=[["Cost", "Length", "DOUBLE", "Shape_Leng"]]
)
# 设置高速公路限制
arcpy.na.SetRestrictionUsage(
in_network_dataset=nd_layer,
restriction_attribute_name="RoadType",
restriction_usage="PROHIBITED",
where_clause="类型 = '高速公路'"
)
# 保存并构建网络数据集
arcpy.na.SaveNetworkDatasetLayer(nd_layer)
arcpy.na.BuildNetwork(nd_path)
print("网络数据集创建成功")
except Exception as e:
arcpy.AddError(f"网络数据集创建失败: {str(e)}")
# 获取详细的错误信息
msgs = arcpy.GetMessages(2)
arcpy.AddError(msgs)
exit()
# 5. 空间关联分析
gas_district = os.path.join(feature_dataset, "GasStations_District")
if not arcpy.Exists(gas_district):
arcpy.analysis.SpatialJoin(
target_features=os.path.join(feature_dataset, "GasStations"),
join_features=os.path.join(feature_dataset, "Districts"),
out_feature_class=gas_district,
join_operation="JOIN_ONE_TO_ONE",
match_option="WITHIN"
)
# 6. 服务区分析
service_areas = os.path.join(feature_dataset, "ServiceAreas")
if not arcpy.Exists(service_areas):
arcpy.na.GenerateServiceAreas(
in_network_dataset=nd_path,
out_feature_class=service_areas,
facilities=gas_district,
break_values=[2000, 5000], # 2km和5km
travel_direction="TRAVEL_FROM",
polygon_overlap="OVERLAP",
merge_by_break_value="MERGE_BY_BREAK_VALUE"
)
# 7. 小区可达性分析
community_access = os.path.join(feature_dataset, "CommunityAccessibility")
if not arcpy.Exists(community_access):
arcpy.analysis.SpatialJoin(
target_features=os.path.join(feature_dataset, "Communities"),
join_features=service_areas,
out_feature_class=community_access,
match_option="WITHIN"
)
# 8. 添加可达性分类字段
if not arcpy.ListFields(community_access, "AccessLevel"):
arcpy.management.AddField(
in_table=community_access,
field_name="AccessLevel",
field_type="TEXT",
field_length=20
)
code_block = """def classify(dist):
if dist <= 2000:
return '高可达性'
elif dist <= 5000:
return '中可达性'
else:
return '低可达性'"""
arcpy.management.CalculateField(
in_table=community_access,
field="AccessLevel",
expression="classify(!NEAR_DIST!)",
expression_type="PYTHON3",
code_block=code_block
)
# 9. 筛选有效小区
valid_communities = os.path.join(feature_dataset, "ValidCommunities")
if not arcpy.Exists(valid_communities):
arcpy.analysis.Select(
in_features=community_access,
out_feature_class=valid_communities,
where_clause="AccessLevel IN ('高可达性', '中可达性')"
)
print("处理完成!结果保存在 ValidCommunities 要素类中")
修正代码中的错误