如何用FloPy在MODFLOW-NWT中提取各应力期水头并计算空间均值

一、优化水头读取逻辑 #

原代码重复调用hds.get_times()会浪费资源，先一次性提取所有应力期时间点：

import numpy as np
import flopy.utils.binaryfile as bf

hds_file = "hist.bhd"
hds = bf.HeadFile(hds_file, text='head', precision='single')
times = hds.get_times()  # 一次性获取所有应力期时间点

二、保存各应力期的水头 #

提供两种实用保存方式：

1. 单个文件保存单应力期水头 #

适合单独查看某一应力期结果，用np.save保留数组维度信息：

for kper in range(len(times)):
    head = hds.get_data(totim=times[kper])
    # 文件名包含应力期编号，方便识别
    np.save(f"head_kper_{kper+1}.npy", head)
    # 若需可读文本格式（小模型适用），可替换为：
    # np.savetxt(f"head_kper_{kper+1}.txt", head.flatten())

2. 单个文件打包所有应力期水头 #

用np.savez打包存储，节省空间且便于批量读取：

all_heads = {}
for kper in range(len(times)):
    head = hds.get_data(totim=times[kper])
    all_heads[f"kper_{kper+1}"] = head

# 保存为npz格式压缩文件
np.savez("all_heads.npz", **all_heads)

# 后续读取方式示例：
# loaded_data = np.load("all_heads.npz")
# head_kper1 = loaded_data["kper_1"]

三、计算空间/时间均值 #

MODFLOW输出的无效水头默认值为-9999.0，计算时需先过滤：

1. 每个应力期的空间平均水头 #

计算单应力期内所有有效单元的水头均值：

spatial_means = []

for kper in range(len(times)):
    head = hds.get_data(totim=times[kper])
    valid_heads = head[head != -9999.0]
    # 避免空数组报错，无有效单元时设为NaN
    mean_head = np.mean(valid_heads) if len(valid_heads) > 0 else np.nan
    spatial_means.append(mean_head)
    print(f"应力期 {kper+1} 空间平均水头: {mean_head:.2f}")

# 保存均值到文本文件
np.savetxt("spatial_mean_heads.txt", spatial_means, fmt="%.2f", header="Stress Period Spatial Mean Head")

2. 每个模型节点的时间平均水头 #

计算单个空间单元在所有应力期的水头均值：

# 读取所有应力期水头，转为四维数组（应力期数×层数×行数×列数）
all_heads_array = np.array([hds.get_data(totim=t) for t in times])

# 创建数组存储每个节点的时间均值，初始值设为NaN
time_mean_per_node = np.full_like(all_heads_array[0], np.nan)

# 遍历所有节点计算均值
for layer in range(time_mean_per_node.shape[0]):
    for row in range(time_mean_per_node.shape[1]):
        for col in range(time_mean_per_node.shape[2]):
            node_heads = all_heads_array[:, layer, row, col]
            valid_mask = node_heads != -9999.0
            if np.any(valid_mask):
                time_mean_per_node[layer, row, col] = np.mean(node_heads[valid_mask])

# 保存节点时间均值
np.save("time_mean_per_node.npy", time_mean_per_node)

完整代码示例 #

import numpy as np
import flopy.utils.binaryfile as bf

hds_file = "hist.bhd"
hds = bf.HeadFile(hds_file, text='head', precision='single')
times = hds.get_times()

# 保存所有应力期水头+计算空间均值
all_heads = {}
spatial_means = []

for kper in range(len(times)):
    head = hds.get_data(totim=times[kper])
    all_heads[f"kper_{kper+1}"] = head
    
    # 计算空间均值
    valid_heads = head[head != -9999.0]
    mean_head = np.mean(valid_heads) if len(valid_heads) > 0 else np.nan
    spatial_means.append(mean_head)
    print(f"应力期 {kper+1}: 空间平均水头 = {mean_head:.2f}")

# 保存结果文件
np.savez("all_heads.npz", **all_heads)
np.savetxt("spatial_mean_heads.txt", spatial_means, fmt="%.2f", header="Stress Period, Spatial Mean Head")

# 计算每个节点的时间均值（可选）
all_heads_array = np.array(list(all_heads.values()))
time_mean_per_node = np.full_like(all_heads_array[0], np.nan)

for layer in range(time_mean_per_node.shape[0]):
    for row in range(time_mean_per_node.shape[1]):
        for col in range(time_mean_per_node.shape[2]):
            node_heads = all_heads_array[:, layer, row, col]
            valid_mask = node_heads != -9999.0
            if np.any(valid_mask):
                time_mean_per_node[layer, row, col] = np.mean(node_heads[valid_mask])

np.save("time_mean_per_node.npy", time_mean_per_node)

内容的提问来源于stack exchange，提问作者Dima