[SUBCATCHMENTS]
;;Name           Rain Gage        Outlet           Area     %Imperv  Width    %Slope   CurbLen  SnowPack        
;;-------------- ---------------- ---------------- -------- -------- -------- -------- -------- ----------------
SC_5838467073    Raingage         J_5838467073     0.048547 100.0    32.38243786916873 3.253286043802897 0                        
SC_1200996354    Raingage         J_1200996354     0.046269 0.0      31.613519069252817 2.6516504287719727 0                        
SC_5838467075    Raingage         J_5838467075     0.14203  79.875   55.38911621174584 4.760003566741943 0         
以上排水管网的inp文件中关于集水区的部分内容，每个节点都对应有一个集水区。请对下面的代码进行拓展，计算并绘制出Horton阶数和平均集水面积（每个阶数的总集水区面积/该阶数的节点数）的双对数图
def parse_inp_file(filepath):
    with open(filepath, 'r') as file:
        content = file.readlines()

    junctions = []
    conduits = []
    current_section = None

    for line in content:
        # 去除行两端的空白字符
        line = line.strip()

        # 跳过空行
        if not line:
            continue

        # Check for section headers
        if line.startswith("["):
            current_section = line.strip().replace("[", "").replace("]", "").strip()
            continue

        # Parse JUNCTIONS section
        if current_section == "JUNCTIONS" and not line.startswith(";;"):
            parts = line.split()
            junctions.append({
                'Name': parts[0],
                'Elevation': float(parts[1]),
                'MaxDepth': float(parts[2]),
                'InitDepth': float(parts[3]),
                'SurDepth': float(parts[4]),
                'Aponded': float(parts[5])
            })
        # Parse OUTFALLS section
        if current_section == "OUTFALLS" and not line.startswith(";;"):
            parts = line.split()
            junctions.append({
                'Name': parts[0],
                'Elevation': float(parts[1]),
            })

        # Parse CONDUITS section
        elif current_section == "CONDUITS" and not line.startswith(";;"):
            parts = line.split()
            conduits.append({
                'Name': parts[0],
                'FromNode': parts[1],
                'ToNode': parts[2],
                'Length': float(parts[3]),
                'Roughness': float(parts[4]),
                'InOffset': float(parts[5]),
                'OutOffset': float(parts[6]),
                'InitFlow': float(parts[7]),
                'MaxFlow': float(parts[8])
            })

    return junctions, conduits


def create_graph(junctions, conduits):
    G = nx.DiGraph()

    # Add junctions to the graph
    for junction in junctions:
        G.add_node(junction['Name'], elevation=junction['Elevation'])

    # Add conduits to the graph as directed edges
    for conduit in conduits:
        G.add_edge(conduit['FromNode'], conduit['ToNode'], length=conduit['Length'], roughness=conduit['Roughness'])

    return G
def compute_horton_strahler(graph):
    # 定义一个字典来存储每个节点的Horton-Strahler顺序
    hs_order = {node: 1 for node in graph.nodes()}  # 初始所有节点的HS顺序为1
    degrees = graph.in_degree()  # 计算入度

    # 从源到出口遍历通过拓扑排序
    nodes = list(nx.topological_sort(graph))  # 获取拓扑排序的节点

    for node in nodes:
        # 从源节点到当前节点的顺序计算
        if degrees[node] == 0:  # 源节点
            continue

        # 获取所有进入当前节点的前驱节点
        incoming = list(graph.predecessors(node))
        orders = [hs_order[pre] for pre in incoming]

        if len(set(orders)) == 1:  # 相同顺序相遇
            hs_order[node] = orders[0] + 1
        else:  # 不同顺序相遇
            hs_order[node] = max(orders)

    return hs_order