如何利用python处理数据满足ce-qual-w2地形文件要求

如何利用python读取地形断面数据，将其处理为CE-QUAL-W2模型可以应用的bth.csv（地形文件)数据格式。

我想要达到的结果

可以通过读取已处理的断面数据，通过python，计算出每一层的宽度。

######断面的数据格式下所示

• rs x z
  1 0 392.07
  1 2.5 391.28
  1 5.1 387.72
  1 6.4 380.54
  1 12.4 377.72
  1 18 375.3
  1 29.8 370.9
  1 37.3 365.9
  1 48.9 363.2
  1 56.8 361.3
  1 71.5 360.9
  1 79.8 358.3
  1 84 358.3
  1 87.3 357.4
  1 95.5 344.7
  1 106.8 342
  1 118.2 340.2
  1 125.1 339.2
  1 130.9 337.7
  1 142.9 324.8
  1 146.9 317.6
  1 154.7 310.1
  1 166.2 310.3
  1 173.7 310.5
  1 174.1 310.5
  1 181.8 310.7
  1 189.7 307.6
  1 197.6 298.9
  1 205.8 294.6
  1 213.8 289.5
  1 221.8 283.7
  1 229.8 275.9
  1 237.8 272.7
  1 245.8 266.6
  1 253.9 263.2
  1 261.7 260.9
  1 269.8 260.2
  1 277.7 260.7
  1 285.7 261.5
  1 293.6 262.2
  1 301.6 262.6
  1 309.6 262.9
  1 317.6 263
  1 325.6 263.1
  1 333.5 262.9
  1 341.6 262.8
  1 349.5 262.8
  1 357.4 262.7
  1 365.3 262.7
  1 373.4 263.3
  1 381.5 263.3
  1 389.4 263.3
  1 397.5 263.2
  1 405.6 263.3
  1 413.5 263.5
  1 421.5 262.8
  1 429.6 262.5
  1 437.5 262.2
  1 445.5 262.2
  1 453.5 262.3
  1 461.5 262.8
  1 469.5 263.6
  1 477.6 264.5
  1 485.5 266.2
  1 493.5 266.9
  1 501.4 269.7
  1 509.4 270.3
  1 517.4 266.5
  1 525.4 266.8
  1 533.4 267.7
  1 541.4 265.5
  1 549.4 265.4
  1 557.3 264.8
  1 565.3 273.8
  1 570.2 286.9
  1 573.5 287.6
  1 581.2 288.7
  1 589.2 290.2
  1 597.2 292.2
  1 605.1 297.1
  1 605.5 301.2
  1 621.2 310
  1 629.2 310.5
  1 635.2 310.8
  1 643.4 311.2
  1 651.3 311
  1 659.7 311.7
  1 667.7 311.8
  1 677.2 312.2
  1 691.4 308.8
  1 701.2 308.2
  1 707 310.3
  1 725 324.7
  1 732.9 326.3
  1 740.8 328.7
  1 749 332.9
  1 757.1 337.6
  1 764.9 339.8
  1 772.9 342.6
  1 780.9 345.8
  1 788.8 348.4
  1 796.9 352.9
  1 804.7 361.2
  1 812.5 366.1
  1 814.9 368.2
  1 826.6 377.72
  1 830.2 382.66
  1 830.7 389.97
  1 831.3 390
  2 0 400.08
  2 10 397.58
  2 16.8 392.92
  2 23 388
  2 26.8 385.41
  2 29.3 383.43
  2 35.9 377.69
  2 39.6 375.3
  2 54.6 362.3
  2 63.1 357.4
  2 72.2 351.1
  2 79.5 348.9
  2 84.3 345.9
  2 91.8 340.6
  2 98.7 338.6
  2 116.2 330.6
  2 124 325.2
  2 132 320.8
  2 140 317.3
  2 148 314.2
  2 155.9 310.7
  2 164 308.6
  2 171.8 307.3
  2 176.7 304
  2 184.7 300.1
  2 200.6 295.1
  2 214.6 291.8
  2 224.7 282.4
  2 232.7 280.9
  2 245.9 276.6
  2 254 274.6
  2 270 268.7
  2 277.9 265.4
  2 285.9 262.6
  2 294 259.7
  2 301.9 254.4
  2 309.9 251.2
  2 317.8 249.7
  2 324.4 249.3
  2 332.3 249.6
  2 340.3 249.7
  2 349.9 256.7
  2 357.9 260.6
  2 365.8 265.4
  2 372.3 269.2
  2 381.8 272.4
  2 389.9 275
  2 396.1 276.3
  2 404.2 277.9

下图是想要的得到的计算效果，计算出每一层的宽度，

目前自己用于数据处理的代码

PS：程序可以设置最低的断面高程高度，如果其他断面高于这个最低值 ，就会用0宽度来填充这层，所以后边代码有个计算0的排列

import numpy as np
from scipy.interpolate import interp1d
import pandas as pd
import openpyxl


#计算 0的排列
# def calc_numb_zero():
#     workbook = openpyxl.Workbook()
#     sheet = workbook.active
#     for x in range(1,len(list_zmin_diff)+1):
#         for i in range(1,list_zmin_diff[x-1]+1):
#             cell = sheet.cell(row=i, column=x)
#             cell.value = 0
#     #保存excel文件
#     workbook.save('填充zero1.xlsx')
# calc_numb_zero()



df = pd.read_excel('cross-section.xlsx')

#断面最大最小值
z_min = df['z'].min()
z_max = df['z'].max()
#各断面最低点与所有断面的最低点的差列表
list_zmin_diff = []

workbook = openpyxl.Workbook()
acive_sheet = workbook.active

for i in range(160,0,-1):
    #读取 RS值不同的断面
    df_rs = df[df['rs']==i]
    #各断面最低点与所有断面的最低点的差
    z_diff = np.ceil(df_rs['z'].min()-z_min).astype(int)
    #找到当前断面的最低点值的索引及对应最低点的x坐标
    min_z_loc = df_rs['z'].idxmin()
    x_min = df_rs.loc[min_z_loc, 'x']
    # x值最小值左侧列表
    list_x_l=np.array(df_rs.loc[:min_z_loc,'x'])
    # # x值最小值右侧列表
    list_x_r=np.array(df_rs.loc[min_z_loc+1:,'x'])
    # # z值最小值左侧列表
    list_z_l=np.array(df_rs.loc[:min_z_loc,'z'])
    # # z值最小值右侧列表
    list_z_r=np.array(df_rs.loc[min_z_loc+1:,'z'])
    # list_zmin_diff.append(z_diff)

    #左右端点索引
    left_point = df_rs.index[0]
    right_point = df_rs.index[-1]
    # 每个断面的活动层数
    active_seg = 191-z_diff
    print(active_seg)

    for j in range(1,active_seg+1):
        #初始两个 x点 用于后续计算 width
        width = 0
        x1 = 0
        x2 = 0
        #创建插值对象
        x1_interp = interp1d(list_z_l, list_x_l)
        x2_interp = interp1d(list_z_r, list_x_r)
        #判断是否超出左端点
        try:
            x1 = x1_interp(df_rs['z'].min() + j)
            if x1 >= df_rs.loc[left_point, 'x']:
                x1 = df_rs.loc[left_point, 'x']
        except ValueError:
            x1 = df_rs.loc[left_point, 'x']

        #判断是否超出右端点
        try:
            x2 = x2_interp(df_rs['z'].min() + j)
            if x2 >= df_rs.loc[right_point, 'x']:
                x2 = df_rs.loc[right_point, 'x']
        except ValueError:
            x2 = df_rs.loc[right_point, 'x']
        width = x2 - x1
        cell = acive_sheet.cell(row=j,column=160-i+1)
        cell.value=width
#
# workbook.save('width.xlsx')

运行结果及详细报错内容

我现在的代码如上所示，但是部分Layer(SEGMENT的层)会出现下层宽度大于上层宽度。并且还会出现部分的最底层几层宽度按为断面的左右端点之差的绝对值，不符合最初的计算想法。

在使用Python读取地形断面数据并将其转换为CE-QUAL-W2模odel的bth.csv格式时，首先需要了解bth.csv的格式要求。CE-QUAL-W2模型的bth.csv文件应该具有以下格式：

Branch, Distance, Layer, Elevation, Area

这里，Branch 是分支或流段的编号，Distance 是从分支的开始处到给定断面的距离，Layer 是从湖底向上的层次编号（湖底是1，然后递增），Elevation 是每层的底部海拔高度，Area 是每层在该断面的面积。

以下是一个Python脚本的例子，用于读取你给出的断面数据并将其转换为bth.csv格式。这个脚本假设断面数据已经按照从湖/河流的上游到下游的顺序进行了排序。由于没有给出每一层的面积，这个脚本假设每个断面的面积是1平方米。

请注意，这个脚本可能需要根据你的具体需求进行调整。

import pandas as pd

# 读取地形断面数据
df = pd.read_csv('断面数据.txt', sep=' ', header=None, names=['rs', 'x', 'z'])

# 创建一个新的DataFrame来保存bth数据
bth = pd.DataFrame(columns=['Branch', 'Distance', 'Layer', 'Elevation', 'Area'])

# 对于每一个不同的分支
for branch in df['rs'].unique():
    # 提取该分支的数据
    branch_df = df[df['rs'] == branch].copy()
    branch_df.sort_values('x', inplace=True)
    # 对于每一个断面
    for distance in branch_df['x'].unique():
        # 提取该断面的数据并按深度进行排序
        segment_df = branch_df[branch_df['x'] == distance].copy()
        segment_df.sort_values('z', inplace=True)
        # 从湖底开始，为每一层添加一行bth数据
        for i, row in enumerate(segment_df.itertuples()):
            bth = bth.append({
                'Branch': branch,
                'Distance': distance,
                'Layer': i + 1,
                'Elevation': row.z,
                'Area': 1  # 假设面积为1
            }, ignore_index=True)

# 保存bth数据到CSV文件
bth.to_csv('bth.csv', index=False)

这个脚本会为每个断面的每一层创建一行bth数据，其中Layer从1开始并逐层递增，Elevation是该层的底部高度，Area被假设为1平方米。所有的bth数据最后被保存为一个CSV文件。

不知道你这个问题是否已经解决, 如果还没有解决的话:

• 以下回答来自chatgpt:

  我将为你提供处理地形断面数据并转化为CE-QUAL-W2模型可用的bth.csv数据格式的解决方案。

  首先，你需要读取地形断面数据，可以使用Python的pandas库来读取csv文件。请确保你已经安装了pandas库。

  import pandas as pd
  
  data = pd.read_csv("地形断面数据.csv")
  

  接下来，我们可以从读取的数据中提取我们需要的列，并对这些数据进行处理。

  x = data["x"]  # 提取x列
  z = data["z"]  # 提取z列
  

  然后，我们需要计算每一层的宽度。我们可以通过计算每一层上下两个z值的差值得到宽度，同时可以应用我们提到的最低断面高程高度的功能。

  widths = []
  for i in range(len(z)-1):
      width = z[i] - z[i+1]
      if width > 0:
          widths.append(width)
      else:
          widths.append(0)
  

  接下来，我们需要将计算得到的宽度写入到bth.csv文件中。我们可以使用Python的csv库来实现。

  import csv
  
  with open("bth.csv", "w", newline="") as file:
      writer = csv.writer(file)
      writer.writerow(["width"])
      for width in widths:
          writer.writerow([width])
  

  以上就是整个处理过程的一个简单示例。你可以根据你的实际需求进行修改和优化代码。请注意，这里只提供了一种解决方案，可能还会有其他更好的方法来处理地形断面数据。

  如果你有任何其他问题或者对以上解决方案有疑问，请随时向我提问。

如果你已经解决了该问题, 非常希望你能够分享一下解决方案, 写成博客, 将相关链接放在评论区, 以帮助更多的人 ^-^