MODFLOW6-FloPy 中是否有任何自动选项可以根据 DEM 和网格单元高度将拒绝的渗透从 UZF 单元路由到 SFR 河段？
问问题
58 次
1 回答
0
以下所有内容均假设您熟悉 python 和 flopy。
我认为下面的 python 会帮助你，我写它是为了基本上完成你所要求的。您可能需要进行一些编辑才能将输入和输出转换为您的特定应用程序所需的形式。
它需要几个输入：
- elev_arr：这是一个包含地表高程的 2D numpy 数组。对于我在其中使用的模型，我有一个保存的文件，其中包含名为“top1.txt”的所有地表高程。我直接在 MODFLOW 数组中使用了这个文件。如果在您的模型中，您的地表高程出现在不同的图层中，则您需要考虑这一点。
- ibnd：给出高程数组中活动单元的范围
- sfr_dat：是一个元组列表，其中每个元组包含 (lay, row, col)。此列表的顺序很重要！该列表应与出现在 SFR 输入文件的 packagedata 块中的数据的顺序相同。
- sfrlayout_conn_candidate_elev：大于所有活动单元格的最大海拔的值
该函数determine_runoff_conns_4mvr生成并返回一个类似于 UZF1 中使用的 IRUNBND 数组的 2D numpy 数组。我把这个函数放到它自己的python脚本中调用mvr_conns，然后导入它。返回的 numpy 数组的每个位置的值determine_runoff_conns_4mvr对应于拒绝渗透将被路由到的 SFR 范围（由其在 SFR PackageData 块中的索引标识），我认为这就是您想要的。
关于湖泊的说明：如果有的话，下面的脚本需要一些修改。我只有 1 个湖，我以一种非常不稳健的方式处理它（因为我可以侥幸逃脱）。如果您有多个湖泊，我必须将此修复留给您作为练习。
因此，对于 MVR 包，您需要处理由determine_runoff_conns_4mvr. 这是我的做法，但您可能需要根据自己的目的进行更改。基本上，在使用 flopy 实例化 MVR 包时，static_mvrperioddata可以将其中包含的信息传递给参数。perioddata
# This function uses the top elevation array and SFR locations to calculate the irunbnd array from the UZF1 package.
import mvr_conns
# Of course in MF6, MVR now handles these runoff connections
irunbnd = mvr_conns.determine_runoff_conns_4mvr(pth) # at this point, irunbnd is 1-based, compensate below
iuzno = 0
k = 0 # Hard-wire the layer no.
for i in range(0, iuzfbnd.shape[0]):
for j in range(0,iuzfbnd.shape[1]):
if irunbnd[i,j] > 0: # This is a uzf -> sfr connection
static_mvrperioddata.append(['UZF-1', iuzno, 'SFR-1', irunbnd[i, j] - 1, 'FACTOR', 1.])
iuzno += 1
elif irunbnd[i, j] < 0: # This is a uzf -> lak connection
static_mvrperioddata.append(['UZF-1', iuzno, 'RES-1', abs(irunbnd[i, j]) - 1, 'FACTOR', 1.])
iuzno += 1
这是基本函数：def determine_runoff_conns_4mvr(pth)：
elev_arr = np.loadtxt(os.path.join(pth,'dis_support','top1.txt'), dtype=np.float)
ibnd = np.loadtxt(os.path.join(pth,'bas_support','ibnd1.txt'), dtype=np.int)
# Get the sfr information stored in a companion script
sfr_dat = sfrinfo.get_sfrlist()
sfrlayout = np.zeros_like(ibnd)
for i, rchx in enumerate(sfr_dat):
row = rchx[1]
col = rchx[2]
sfrlayout[row - 1, col - 1] = i + 1
sfrlayout_new = sfrlayout.copy()
nrow = 64
ncol = 133
stop_candidate = False
for i in np.arange(0, nrow):
for j in np.arange(0, ncol):
# Check to ensure current cell is active
if ibnd[i, j] == 0:
continue
# Check to make sure it is not a stream cell
if not sfrlayout[i, j] == 0:
continue
# Recursively trace path by steepest decent back to a stream
curr_i = i
curr_j = j
sfrlayout_conn_candidate_elev = 10000.
while True:
direc = 0
min_elev = elev_arr[curr_i, curr_j]
# Look straight left
if curr_j > 0:
if not sfrlayout[curr_i, curr_j - 1] == 0 and not ibnd[curr_i, curr_j - 1] == 0: # Step in if neighbor is a stream cell
if elev_arr[curr_i, curr_j - 1] > 0 and (elev_arr[curr_i, curr_j - 1] < elev_arr[curr_i, curr_j] and
elev_arr[curr_i, curr_j - 1] < sfrlayout_conn_candidate_elev):
sfrlayout_conn_candidate = sfrlayout[curr_i, curr_j - 1]
sfrlayout_conn_candidate_elev = elev_arr[curr_i, curr_j - 1]
stop_candidate = True
elif not elev_arr[curr_i, curr_j - 1] == 0 and not ibnd[curr_i, curr_j - 1] == 0: # Step here if neighbor is not an sfr cell
if elev_arr[curr_i, curr_j - 1] < elev_arr[curr_i, curr_j] and elev_arr[curr_i, curr_j - 1] < min_elev:
elevcm1 = elev_arr[curr_i, curr_j - 1]
min_elev = elevcm1
direc = 2
# Look up and left
if curr_j > 0 and curr_i > 0:
if not sfrlayout[curr_i - 1, curr_j - 1] == 0 and not ibnd[curr_i - 1, curr_j - 1] == 0: # Step in if neighbor is a stream cell
if elev_arr[curr_i - 1, curr_j - 1] > 0 and (elev_arr[curr_i - 1, curr_j - 1] < elev_arr[curr_i, curr_j] and
elev_arr[curr_i - 1, curr_j - 1] < sfrlayout_conn_candidate_elev):
sfrlayout_conn_candidate = sfrlayout[curr_i - 1, curr_j - 1]
sfrlayout_conn_candidate_elev = elev_arr[curr_i - 1, curr_j - 1]
stop_candidate = True
elif not elev_arr[curr_i - 1, curr_j - 1] == 0 and not ibnd[curr_i - 1, curr_j - 1] == 0: # Step here if neighbor is not an sfr cell
if elev_arr[curr_i - 1, curr_j - 1] < elev_arr[curr_i, curr_j] and elev_arr[curr_i - 1, curr_j - 1] < min_elev:
elevrm1cm1 = elev_arr[curr_i - 1, curr_j - 1]
min_elev = elevrm1cm1
direc = 5
# Look straight right
if curr_j < ncol - 1:
if not sfrlayout[curr_i, curr_j + 1] == 0 and not ibnd[curr_i, curr_j + 1] == 0: # Step in if neighbor is a stream cell
if elev_arr[curr_i, curr_j + 1] > 0 and (elev_arr[curr_i, curr_j + 1] < elev_arr[curr_i, curr_j] and
elev_arr[curr_i, curr_j + 1] < sfrlayout_conn_candidate_elev):
sfrlayout_conn_candidate = sfrlayout[curr_i, curr_j + 1]
sfrlayout_conn_candidate_elev = elev_arr[curr_i, curr_j + 1]
stop_candidate = True
elif not elev_arr[curr_i, curr_j + 1] == 0 and not ibnd[curr_i, curr_j + 1] == 0: # Step here if neighbor is not an sfr cell
if elev_arr[curr_i, curr_j + 1] < elev_arr[curr_i, curr_j] and elev_arr[curr_i, curr_j + 1] < min_elev:
elevcm1 = elev_arr[curr_i, curr_j + 1]
min_elev = elevcm1
direc = 4
# Look straight right and down
if curr_i < nrow - 1 and curr_j < ncol - 1:
if not sfrlayout[curr_i + 1, curr_j + 1] == 0 and not ibnd[curr_i + 1, curr_j + 1] == 0: # Step in if neighbor is a stream cell
if elev_arr[curr_i + 1, curr_j + 1] > 0 and (elev_arr[curr_i + 1, curr_j + 1] < elev_arr[curr_i, curr_j] and
elev_arr[curr_i + 1, curr_j + 1] < sfrlayout_conn_candidate_elev):
sfrlayout_conn_candidate = sfrlayout[curr_i + 1, curr_j + 1]
sfrlayout_conn_candidate_elev = elev_arr[curr_i + 1, curr_j + 1]
stop_candidate = True
elif not elev_arr[curr_i + 1, curr_j + 1] == 0 and not ibnd[curr_i + 1, curr_j + 1] == 0: # Step here if neighbor is not an sfr cell
if elev_arr[curr_i + 1, curr_j + 1] < elev_arr[curr_i, curr_j] and elev_arr[curr_i + 1, curr_j + 1] < min_elev:
elevrp1cp1 = elev_arr[curr_i + 1, curr_j + 1]
min_elev = elevrp1cp1
direc = 7
# Look straight up
if curr_i > 0:
if not sfrlayout[curr_i - 1, curr_j] == 0 and not ibnd[curr_i - 1, curr_j] == 0: # Step in if neighbor is a stream cell
if elev_arr[curr_i - 1, curr_j] > 0 and (elev_arr[curr_i - 1, curr_j] < elev_arr[curr_i, curr_j] and
elev_arr[curr_i - 1, curr_j] < sfrlayout_conn_candidate_elev):
sfrlayout_conn_candidate = sfrlayout[curr_i - 1, curr_j]
sfrlayout_conn_candidate_elev = elev_arr[curr_i - 1, curr_j]
stop_candidate = True
elif not elev_arr[curr_i - 1, curr_j] == 0 and not ibnd[curr_i - 1, curr_j] == 0: # Step here if neighbor is not an sfr cell
if elev_arr[curr_i - 1, curr_j] < elev_arr[curr_i, curr_j] and elev_arr[curr_i - 1, curr_j] < min_elev:
elevcm1 = elev_arr[curr_i - 1, curr_j]
min_elev = elevcm1
direc = 3
# Look up and right
if curr_i > 0 and curr_j < ncol - 1:
if not sfrlayout[curr_i - 1, curr_j + 1] == 0 and not ibnd[curr_i - 1, curr_j + 1] == 0: # Step in if neighbor is a stream cell
if elev_arr[curr_i - 1, curr_j + 1] > 0 and (elev_arr[curr_i - 1, curr_j + 1] < elev_arr[curr_i, curr_j] and
elev_arr[curr_i - 1, curr_j + 1] < sfrlayout_conn_candidate_elev):
sfrlayout_conn_candidate = sfrlayout[curr_i - 1, curr_j + 1]
sfrlayout_conn_candidate_elev = elev_arr[curr_i - 1, curr_j + 1]
stop_candidate = True
elif not elev_arr[curr_i - 1, curr_j + 1] == 0 and not ibnd[curr_i - 1, curr_j + 1] == 0: # Step here if neighbor is not an sfr cell
if elev_arr[curr_i - 1, curr_j + 1] < elev_arr[curr_i, curr_j] and elev_arr[curr_i - 1, curr_j + 1] < min_elev:
elevrm1cp1 = elev_arr[curr_i - 1, curr_j + 1]
min_elev = elevrm1cp1
direc = 6
# Look straight down
if curr_i < nrow - 1:
if not sfrlayout[curr_i + 1, curr_j] == 0 and not ibnd[curr_i + 1, curr_j] == 0: # Step in if neighbor is a stream cell
if elev_arr[curr_i + 1, curr_j] > 0 and (elev_arr[curr_i + 1, curr_j] < elev_arr[curr_i, curr_j] and
elev_arr[curr_i + 1, curr_j] < sfrlayout_conn_candidate_elev):
sfrlayout_conn_candidate = sfrlayout[curr_i + 1, curr_j]
sfrlayout_conn_candidate_elev = elev_arr[curr_i + 1, curr_j]
stop_candidate = True
elif not elev_arr[curr_i + 1, curr_j] == 0 and not ibnd[curr_i + 1, curr_j] == 0: # Step here if neighbor is not an sfr cell
if elev_arr[curr_i + 1, curr_j] < elev_arr[curr_i, curr_j] and elev_arr[curr_i + 1, curr_j] < min_elev:
elevrp1 = elev_arr[curr_i + 1, curr_j]
min_elev = elevrp1
direc = 1
# Look down and left
if curr_i < nrow - 1 and curr_j > 0:
if not sfrlayout[curr_i + 1, curr_j - 1] == 0 and not ibnd[curr_i + 1, curr_j - 1] == 0: # Step in if neighbor is a stream cell
if elev_arr[curr_i + 1, curr_j - 1] > 0 and (elev_arr[curr_i + 1, curr_j - 1] < elev_arr[curr_i, curr_j] and
elev_arr[curr_i + 1, curr_j - 1] < sfrlayout_conn_candidate_elev):
sfrlayout_conn_candidate = sfrlayout[curr_i + 1, curr_j - 1]
sfrlayout_conn_candidate_elev = elev_arr[curr_i + 1, curr_j - 1]
stop_candidate = True
elif not elev_arr[curr_i + 1, curr_j - 1] == 0 and not ibnd[curr_i + 1, curr_j - 1] == 0: # Step here if neighbor is not an sfr cell
if elev_arr[curr_i + 1, curr_j - 1] < elev_arr[curr_i, curr_j] and elev_arr[curr_i + 1, curr_j - 1] < min_elev:
elevrp1cm1 = elev_arr[curr_i + 1, curr_j - 1]
min_elev = elevrp1cm1
direc = 8
# if stop candidate found, don't move the cell indices
if not stop_candidate:
# Direc corresponds to:
# |----------------------
# | 5 | 3 | 6 |
# |----------------------
# | 2 | cur_cel | 4 |
# |----------------------
# | 8 | 1 | 7 |
# |----------------------
if direc == 0:
break
elif direc == 1:
curr_i += 1
elif direc == 2:
curr_j -= 1
elif direc == 3:
curr_i -= 1
elif direc == 4:
curr_j += 1
elif direc == 5:
curr_i -= 1
curr_j -= 1
elif direc == 6:
curr_i -= 1
curr_j += 1
elif direc == 7:
curr_i += 1
curr_j += 1
elif direc == 8:
curr_i += 1
curr_j -= 1
if stop_candidate:
sfrlayout_new[i, j] = sfrlayout_conn_candidate
stop_candidate = False
break # Bust out of while loop
elif not stop_candidate:
# Check if encountered ibnd == 0, which may be a lake or boundary that drains out of model
if ibnd[curr_i, curr_j] == 0:
# This condition is dealt with after looping through all cells,
# see comment that starts, "Last step is set..."
break
pass # Commence next downstream cell search
# Last step is set the 0's in the vicinity of the lake equal to the negative of the lake connection
for i in np.arange(0, nrow):
for j in np.arange(0, ncol):
if sfrlayout_new[i, j] == 0 and ibnd[i, j] > 0:
sfrlayout_new[i, j] = -1
# Once all cells are filled, save to array
np.savetxt(os.path.join(pth, 'uzf_support','irunbnd_mf6.txt'), sfrlayout_new, fmt='%5d')
return sfrlayout_new
于 2020-06-24T20:24:34.737 回答