我有一个输出的工作模型。我在试着利用弗罗比来询问结果。具体而言,我希望能够为特定的单元格生成(动态)区域预算信息。当然,我可以麻烦地构建我感兴趣的区域文件,但我希望能够提取结果,因为我需要它们,而不必经过这一中间步骤。
我正在寻找一个工作流,我可以在Python中使用该工作流从CBB (已经生成)中提取特定单元格的净流信息(希望避免生成区域文件、导入和提取区域文件的遗留风格)。
编辑
试图让这个实用程序开始工作时遇到了问题。这是最新的尝试。不知道fl皮是否被CBB本身窒息了,或者我提供的区域是否给它带来了问题。
该模型为7层196列241列。我试着提取Layer=3,Row=58,Col=30。我为区域创建的list对象如下:
zon_lst = []
for lay in range(7):
for row in range(196):
for col in range(241):
if lay == 2 and row == 57 and col == 29:
zon_lst.append(1)
else:
zon_lst.append(0)然后,我使用以下内容将其写入一个文件:
def chunk_list(alist, n):
for i in range(0, len(alist), n):
yield alist[i:i + n]
def zon_gen(mylst, rows, cols, file_cols):
# Assumes integers for zonation file
frmt = '{:>3}'
list_by_lay = chunk_list(mylst, rows * cols)
astr = '\n'.join(['\n'.join([''.join([frmt.format(i) for i in seq]) for seq in chunk_list(layer, file_cols)]) for layer in list_by_lay])
return astr
zon_str = zon_gen(zon_lst, 196, 241, 10)
with open('t_26.zon', 'w+') as file:
file.write('# For Scoping Work\n')
file.write('1\n')
file.write('t_26\n')
file.write('INTERNAL 1 (10I3) 1\n')
file.write(zon_str)然后,我为区域预算类/方法构建了modflow模型:
import flopy
mf = flopy.modflow.Modflow(modelname='t_26_scope', version='mf2k')
zon = flopy.modflow.ModflowZon.load(r"t_26.zon",mf,nrow=196, ncol=241)
zb = flopy.utils.zonbud.ZoneBudget(r"P2Rv8.2_1000yr.cbb", zon)所有这些都很好地运行到最后一个命令,其中我得到了以下错误:
Traceback (most recent call last):
File "<input>", line 1, in <module>
File "C:\ProgramData\Miniconda3\envs\ca_cie\lib\site-packages\flopy\utils\zonbud.py", line 53, in __init__
self.cbc = CellBudgetFile(cbc_file)
File "C:\ProgramData\Miniconda3\envs\ca_cie\lib\site-packages\flopy\utils\binaryfile.py", line 618, in __init__
self._build_index()
File "C:\ProgramData\Miniconda3\envs\ca_cie\lib\site-packages\flopy\utils\binaryfile.py", line 708, in _build_index
self._skip_record(header)
File "C:\ProgramData\Miniconda3\envs\ca_cie\lib\site-packages\flopy\utils\binaryfile.py", line 768, in _skip_record
self.file.seek(nbytes, 1)
OSError: [Errno 22] Invalid argument请记住,我仍然有兴趣跳过文件写入部分。我把它包括进来是为了展示我的工作到目前为止
发布于 2019-10-29 20:14:42
ZoneBudget类所需的矩阵必须具有以下维度(泛型术语):
ndarray_dim = lays * rows * cols在我这里提出的示例问题中,相应的维度是:7* 196 * 241,结果是:
len(my_array)
>>> 7
len(my_array[0])
>>> 196
len(my_array[0][0])
>>> 241正如前面指出的,我需要找出为什么逐单元预算二进制文件的二进制读取器不能读取文件的原因,因此我将准备好这一点,作为以后的另一次讨论。谢谢你幽默这个问题“杰森·贝利诺!
另外,区域文件生成器不会生成块样式的文件,就像flopy所期望的那样。我在这里添加了更新的代码:
def chunk_list(alist, n):
for i in range(0, len(alist), n):
yield alist[i:i + n]
def block_gen(mylst, rows, cols, file_cols, field_len):
# mylst is a 1D array whose length matches that of the number of cells in the model
# rows, cols: These are the total rows, columns; respectively
# Assumes integers for zonation file
frmt = '{:>%d}' % field_len
zon_str = ''
for lay in chunk_list(mylst, (rows * cols)):
zon_str += 'INTERNAL ({:}I{})\n'.format(file_cols, field_len)
for block in chunk_list(lay, cols):
for line in chunk_list(block, file_cols):
zon_str += ''.join([frmt.format(cell) for cell in line]) + '\n'
return zon_str
def write_zb_zonefile(filepath, lays, rows, cols, zon_str):
with open(filepath, 'w+') as file:
file.write('{:<6}{:<6}{:<6}\n'.format(lays, rows, cols))
file.write(zon_str)
return如果将生成的文件与flopy.utils.zonbud.read_zbarray(<filepath>)一起使用,您将得到一个可用的zone_dict,以便在flopy.utils.zonbud.ZoneBudget类中应用。
由于我希望避免每次需要一个新的区域预算时都要编写一个文件,所以我修改了block_gen()函数,并创建了一个可以在ZoneBudget类中直接使用的文件:
from numpy import array
def arr_block_gen(mylst, rows, cols):
# mylst: 1D array containing zone definition for every cell in the model
# rows: total number of rows in the model
# cols: total number of columns in the model
zon_arr = []
lay_idx = 0
for lay in chunk_list(mylst, (rows * cols)):
row_idx = 0
zon_arr.append([])
for block in chunk_list(lay, cols):
zon_arr[lay_idx].append(block)
row_idx += 1
lay_idx += 1
zon_arr = array(zon_arr)
return zon_arrUPDATE最终证明了该方法失败的原因是该方法的默认行为是将CBB文件路径传递给二进制文件读取器类,而二进制文件读取器类的默认操作是以单精度读取CBB文件。因此,我使用二进制文件读取器(传递"double“作为精度)创建CBB对象,然后显式地传递CBB对象,而不是像我在问题中显示的那样传递文件。
发布于 2019-10-29 13:35:53
可以使用ZoneBudget类flopy提取预算信息。它接受一个numpy数组作为输入;为了向后兼容,还包含了read_zbarray和write_zbarray配套实用程序。如果你遇到一个问题,在这里张贴一个最低限度的工作例子,我会看看我是否能帮助你。
https://stackoverflow.com/questions/58599117
复制相似问题