# -*- coding: utf-8 -*- """ Created on Wed Mar 20 14:44:36 2024 @author: ZMK """ import numpy as np import flopy from openpyxl import load_workbook import os import Base as base import CalHead #河流的分段下标 riv_seg={0:[454,479],1:[480,505],2:[506,527],3:[528,562]} riv_seg_celles=[26,26,22,35] #边界的分段下标 ibound_seg={0:[1,86],1:[87,111],2:[112,142],3:[143,170],4:[171,240], 5:[241,282],6:[283,354],7:[355,393],8:[394,436],9:[437,453]} def base_excel(model_name,types): if types =='降雨': return rch_excel(model_name) if types =='河流': return river_excel(model_name) if types =='开采量': return well_excel(model_name) if types =='边界': return ibound_excel(model_name) return [] #降雨cxcel def rch_excel(model_name): paths = base.model_dir + model_name +"\\降雨.xlsx" data=[] if not os.path.exists(paths): return data wb = load_workbook(filename = paths) ws = wb[wb.sheetnames[0]] for row in ws.iter_rows(): tmp =[] for cell in row: tmp.append(cell.value) data.append(tmp) wb.close() #分析降水并且更新 rch_analysis(data,model_name) return data def rch_analysis(data_array,model_name): periods =CalHead.get_model_period(model_name) #周期数 periods_len= len(periods) array = np.asarray(data_array) fid= array[1:17,2] params = array[1:17,3] float_params = np.asarray(params,dtype=float) #取出来 1-16行 4-col_last列的数据 col_last = 4 + periods_len data = array[1:17,4:col_last] float_data = np.asarray(data,dtype=float) for i in range(0,len(float_data)): for j in range(0,len(float_data[i])): float_data[i][j] = float_data[i][j] * float_params[i] /30/100 rch_update(float_data,fid,model_name,periods_len) #更新模型rch数据 def rch_update(float_data,fids,model_name,periods_len): if model_name == base.not_allowed_model or model_name in base.archive_models: return "本模型为验证模型,不允许修改！" model_ws = base.model_dir + model_name mx = flopy.modflow.Modflow.load("modflow.nam", model_ws = model_ws, exe_name="mf2005", verbose=True, version="mf2005", check=False) areas= base.getAreas() for per in range(periods_len): #rch 一个周期的数据 item = mx.rch.rech.__getitem__(kper=per) array2d = item.get_value() per_data = float_data[:,per] print(per_data) for i in range(0,len(per_data)): data =round(float(per_data[i]),8) fid = fids[i] tuples = areas[fid] for entity in tuples: x = entity[0] y = entity[1] array2d[x][y]= data mx.rch.rech.__setitem__(key=per, value=array2d) rch = flopy.modflow.ModflowRch(mx,nrchop=mx.rch.nrchop, ipakcb=mx.rch.ipakcb, rech=mx.rch.rech, irch =mx.rch.irch) rch.write_file(check=False) #河流数据解析 def river_excel(model_name): paths = base.model_dir + model_name +"\\河流.xlsx" data=[] if not os.path.exists(paths): return data wb = load_workbook(filename = paths) ws = wb[wb.sheetnames[0]] for row in ws.iter_rows(): tmp =[] for cell in row: tmp.append(cell.value) data.append(tmp) wb.close() result =[] for i in range(1,len(data)): result.append(data[i]) river_analysis(result,model_name) return result #分析河流数据 def river_analysis(data_array,model_name): periods =CalHead.get_model_period(model_name) periods_len= len(periods) row_last= 1 + periods_len array = np.asarray(data_array) data = array[1:row_last,1:2] data2 = array[1:row_last,1:4] merge_array = np.concatenate((data,data2),axis=1) params = array[1:row_last,5:9] float_data = np.asarray(merge_array,dtype=float) float_params = np.asarray(params,dtype=float) for i in range(0,len(float_data)): for j in range(0,len(float_data[i])): temp = round(float_data[i][j] * float_params[i][j] *30*86400 ,4) last = round(temp/riv_seg_celles[j]/30,4) float_data[i][j]=last print(float_data) river_update(float_data,model_name,periods_len) #更新河流数据 def river_update(float_data,model_name,periods_len): if model_name == base.not_allowed_model or model_name in base.archive_models: return "本模型为验证模型,不允许修改！" model_ws = base.model_dir + model_name ml = flopy.modflow.Modflow.load("modflow.nam", model_ws = model_ws, exe_name="mf2005", verbose=True, version="mf2005", check=False) #循环设置wel文件，更新数据 # lrcq = {0:[[2, 3, 4, -100.]], 1:[[2, 3, 4, -100.]]} lrcq = {} for per in range(periods_len): wel = [] wel = ml.wel.stress_period_data.__getitem__(kper=per) #河流分段数据 segment0 = float_data[per][0] segment1 = float_data[per][1] segment2 = float_data[per][2] segment3 = float_data[per][3] #存储每个应力期的数据 array2d = [] count = 1 for Layer, Row, Column, Q in wel: array = [] if count>= riv_seg[0][0] and count <= riv_seg[0][1]: array = [Layer, Row, Column, segment0] elif count>= riv_seg[1][0] and count <= riv_seg[1][1]: array = [Layer, Row, Column, segment1] elif count>= riv_seg[2][0] and count <= riv_seg[2][1]: array = [Layer, Row, Column, segment2] elif count>= riv_seg[3][0] and count <= riv_seg[3][1]: array = [Layer, Row, Column, segment3] else: array = [Layer, Row, Column, Q] array2d.append(array) count +=1 lrcq[per] = array2d flopy.modflow.ModflowWel(ml,ipakcb= ml.wel.ipakcb, dtype=ml.wel.dtype, options=ml.wel.options, stress_period_data=lrcq) ml.write_input() #解析开采量的数据 def well_excel(model_name): paths = base.model_dir + model_name +"\\开采量.xlsx" data=[] if not os.path.exists(paths): return data wb = load_workbook(filename = paths) ws = wb[wb.sheetnames[0]] for row in ws.iter_rows(): tmp =[] for cell in row: tmp.append(cell.value) data.append(tmp) wb.close() result =[] for i in range(1,len(data)): result.append(data[i]) well_analysis(result,model_name) return result #开采量的分析 def well_analysis(data_array,model_name): periods =CalHead.get_model_period(model_name) periods_len= len(periods) row_last= 1 + periods_len array = np.asarray(data_array) data = array[1:row_last,1:5] float_data= np.asarray(data,dtype=float) #转置数组 zz_array = float_data.transpose() zz_data = [] for i in range(50): zz_data.append(zz_array[0]) for i in range(49): zz_data.append(zz_array[1]) for i in range(18): zz_data.append(zz_array[2]) for i in range(12): zz_data.append(zz_array[3]) zz_data.append(zz_array[0]) np_data = np.asarray(zz_data,dtype=float) well_scale = np.loadtxt(base.well_scale_path, dtype=str) float_scale= np.asarray(well_scale,dtype=float) ##抽水量 = 各城区月开采量*街道比例*10000/(街道井数*30) for i in range(0,len(np_data)): for j in range(0,len(np_data[i])): tmp = np_data[i][j] * float_scale[i][0]*10000/( float_scale[i][1] *30) np_data[i][j] = round(tmp,4) well_update(np_data,model_name,periods_len) #更新开采量 def well_update(np_data,model_name,periods_len): if model_name == base.not_allowed_model or model_name in base.archive_models: return "本模型为验证模型,不允许修改！" model_ws = base.model_dir + model_name ml = flopy.modflow.Modflow.load("modflow.nam", model_ws = model_ws, exe_name="mf2005", verbose=True, version="mf2005", check=False) #循环设置wel文件，更新数据 # lrcq = {0:[[2, 3, 4, -100.]], 1:[[2, 3, 4, -100.]]} lrcq = {} for per in range(periods_len): wel = [] wel = ml.wel.stress_period_data.__getitem__(kper=per) per_data = np_data[:,per] #存储每个应力期的数据 array2d = [] #well文件中 前562个网格不是井 # Layer= wel[i][0] # Row= wel[i][1] # Col = wel[i][2] # Q = wel[i][3] for i in range(0,562): array = [wel[i][0],wel[i][1], wel[i][2],wel[i][3]] array2d.append(array) for i in range(562,len(wel)): indexid = i-562 update_data=per_data[indexid] array = [wel[i][0],wel[i][1], wel[i][2],update_data] array2d.append(array) lrcq[per] = array2d flopy.modflow.ModflowWel(ml,ipakcb= ml.wel.ipakcb,dtype=ml.wel.dtype, options=ml.wel.options, stress_period_data=lrcq) ml.write_input() #解析边界的数据 def ibound_excel(model_name): paths = base.model_dir + model_name +"\\边界.xlsx" data=[] if not os.path.exists(paths): return data wb = load_workbook(filename = paths) ws = wb[wb.sheetnames[0]] for row in ws.iter_rows(): tmp =[] for cell in row: tmp.append(cell.value) data.append(tmp) wb.close() result =[] for i in range(1,len(data)): result.append(data[i]) np_array = ibound_analysis(result,model_name) #返回页面展示的内容数组 view_data = ibound_view_data(np_array) periods =CalHead.get_model_period(model_name) periods_len= len(periods) ibound_update(np_array,model_name,periods_len) return view_data #更新边界 def ibound_update(np_array,model_name,periods_len): if model_name == base.not_allowed_model or model_name in base.archive_models: return "本模型为验证模型,不允许修改！" model_ws = base.model_dir + model_name ml = flopy.modflow.Modflow.load("modflow.nam", model_ws = model_ws, exe_name="mf2005", verbose=True, version="mf2005", check=False) #循环设置wel文件，更新数据 # lrcq = {0:[[2, 3, 4, -100.]], 1:[[2, 3, 4, -100.]]} lrcq = {} for per in range(periods_len): wel = [] wel = ml.wel.stress_period_data.__getitem__(kper=per) per_data = np_array[:,per] #存储每个应力期的数据 array2d = [] count = 1 for Layer, Row, Column, Q in wel: array = [] if count>= ibound_seg[0][0] and count <= ibound_seg[0][1]: array = [Layer, Row, Column, per_data[0]] elif count>= ibound_seg[1][0] and count <= ibound_seg[1][1]: array = [Layer, Row, Column, per_data[1]] elif count>= ibound_seg[2][0] and count <= ibound_seg[2][1]: array = [Layer, Row, Column, per_data[2]] elif count>= ibound_seg[3][0] and count <= ibound_seg[3][1]: array = [Layer, Row, Column, per_data[3]] elif count>= ibound_seg[4][0] and count <= ibound_seg[4][1]: array = [Layer, Row, Column, per_data[4]] elif count>= ibound_seg[5][0] and count <= ibound_seg[5][1]: array = [Layer, Row, Column, per_data[5]] elif count>= ibound_seg[6][0] and count <= ibound_seg[6][1]: array = [Layer, Row, Column, per_data[6]] elif count>= ibound_seg[7][0] and count <= ibound_seg[7][1]: array = [Layer, Row, Column, per_data[7]] elif count>= ibound_seg[8][0] and count <= ibound_seg[8][1]: array = [Layer, Row, Column, per_data[8]] elif count>= ibound_seg[9][0] and count <= ibound_seg[9][1]: array = [Layer, Row, Column, per_data[9]] else: array = [Layer, Row, Column, Q] array2d.append(array) count +=1 lrcq[per] = array2d flopy.modflow.ModflowWel(ml,ipakcb= ml.wel.ipakcb, dtype=ml.wel.dtype, options=ml.wel.options, stress_period_data=lrcq) ml.write_input() #边界的分析 def ibound_analysis(data_array,model_name): array = np.asarray(data_array) data = array[1:10,2:14] #大台 南辛房 上苇甸 三家店 麻裕 温泉 海淀 阳坊 沙河 float_data= np.asarray(data,dtype=float) # 流入1-1 (1-6)大台 南辛房 上苇甸 三家店 麻裕 温泉 result =[] data1= float_data[[0,1,2,3,4,5]] sum1 = data1.sum(axis = 0) result.append(sum1) # 流入1-2()上苇甸 data2= float_data[2] result.append(data2) #流入1-3() 大台 data3 = float_data[0] result.append(data3) # 流出1-1 大台 南辛房 data4 = float_data[[0,1]] sum4 = data4.sum(axis = 0) result.append((-sum4)) # 流出1-2 大台 data5 = float_data[0] result.append((-data5)) #流入3-1 大台 data6 = float_data[0] result.append(data6) # 流入3-2 大台 南辛房 上苇甸 data7 = float_data[[0,1,2]] sum7 = data7.sum(axis = 0) result.append(sum7) # 流入3-3 大台 南辛房 上苇甸 # data8= float_data[[0,1,2]] result.append(sum7) # 流入3-4 大台 南辛房 上苇甸 # data9= float_data[[0,1,2]] result.append(sum7) # 流入3-5 大台 南辛房 上苇甸 # data10= float_data[[0,1,2]] result.append(sum7) np_data = np.asarray(result,dtype=float) np_data = np.around(np_data, decimals=2) return np_data def ibound_view_data(np_data): names=['流入1-1','流入1-2','流入1-3','流出1-1','流出1-2', '流入3-1','流入3-2','流入3-3','流入3-4','流入3-5'] row_sums = np.sum(np_data, axis=1) row_sums= np.around(row_sums,2) params=[1,1,1,2.5,1,0.5,0.5,0.5,0.5,0.5] cells=[86,47,31,28,50,42,72,39,43,17] per=12 x=30 water=[] for i in range(0,len(names)): tmp = round( row_sums[i] * params[i] * cells[i] * per *x , 2) water.append(tmp) arr=[] arr.append(names) arr.append(row_sums) arr.append(params) arr.append(water) str_np = np.asarray(arr,dtype=str) zz= str_np.transpose() title =['边界','降雨量','系数','输入值'] result =[] result.append(title) for item in zz: result.append(item) result = np.asarray(result,dtype=str).tolist() return result