# 导入Flask类 import numpy as np import flopy import flopy.utils.binaryfile as bf import os import json import Base as base import CalHead import ModelPeriod import OpenExcel base_init_year=["2020","2021","2022"] river_start_index = 454 river_end_index =562 #预测周期数 predict_per = 12 #降水量 base_water = base.prefix + 'base_water.ini' def predict_water_chart(base_year,start_time ,end_time,value): water_array = np.loadtxt(base_water, dtype=str,encoding='utf-8') y_data=[] x_data= ModelPeriod.get_months_in_year() water= water_array[int(base_year)] for data in water: float_data = float(data)/9 float_data= round(float_data*value,2) y_data.append(float_data) result = {"y_data": y_data, "x_data": x_data} return result #河流的折线图 base_river = base.prefix + 'base_river.ini' def predict_river_chart(base_year,start_time ,end_time,value): river_array = np.loadtxt(base_river, dtype=str,encoding='utf-8') y_data=[] x_data= ModelPeriod.get_months_in_year() for data in river_array: #万立方米 float_data = float(data)/4/10000 float_data= round(float_data*value,2) y_data.append(float_data) result = {"y_data": y_data, "x_data": x_data} return result #开采量折线图 base_mining = base.prefix + 'base_mining.ini' def predict_well_chart(base_year,start_time ,end_time,data): # 2022/海淀/昌平/门头沟/石景山 area_names=["全部区域","海淀区","昌平区","门头沟区","石景山区"] mining_array = np.loadtxt(base_mining, dtype=str,encoding='utf-8') x_data= ModelPeriod.get_months_in_year() result = {"x_data": x_data} y_data=[] for item in data: name = item["area_name"] value = float(item["value"]) if name =='朝阳区': continue index = area_names.index(name) row_data = mining_array[index] float_data=[] for item in row_data: x = round(float(item)/10000*value,2) float_data.append(x) dicts={"name":name,"data":float_data} y_data.append(dicts) result["y_data"] = y_data return result def run_model(model_name): if model_name == base.not_allowed_model or model_name in base.archive_models: return "本模型为验证模型,不允许修改！" dicts= run_model_predict(model_name) if dicts["code"] == 400: return dicts #导出csv文件 CalHead.exportCsV(model_name) # #更新模型三维网格配置 base.updateModelConfig(model_name) # #创建模型的三维网格 filedir = base.model3d_path + model_name if not os.path.exists(filedir): os.makedirs(filedir, exist_ok=True) base.callModelexe() #计算水资源量和水均衡 CalHead.run_zonebudget_bal(model_name) CalHead.run_zonebudget_res(model_name) dicts= {"code":200,"msg":"预测模型运行成功！" } return dicts def run_model_predict(model_name): predictiondata="" prediction_path = base.model_dir + model_name +"\\prediction.json" if os.path.exists(prediction_path): with open(prediction_path,encoding='utf-8') as f: predictiondata = json.load(f) periods =CalHead.get_model_period(model_name) periods_len= len(periods) if predictiondata: updateDisFile(model_name,periods_len) try: flag = updateBase6File(model_name,predictiondata) if bool(flag)==False: dicts= {"code":400,"msg":"请检查模型的初始水是否设置！" } return dicts except: dicts= {"code":400,"msg":"请检查模型的初始水是否设置！" } return dicts try: updateRchFile(model_name,predictiondata) except: print("RchFile无预测参数，无需修改！") try: updateRiverFile(model_name,predictiondata) except: print("RiverFile无预测参数，无需修改！") else: print("prediction.json 预测场景文件为空，无需更改相应文件") 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) ml.run_model(report = True) dicts= {"code":200,"msg":"预测模型运行成功！" } return dicts #获取 area的 name--> ratio 的结构 def get_area_dict(area): result ={} for i in range(len(area)): name = area[i]["name"] rt = area[i]["ratio"] result[name]= rt return result #获取区县的 row+column --> name结构 def get_distric_dict(): data = base.district result = {} for row ,column ,id ,name in data: key = str(row)+","+str(column) result[key]= name return result #根据 row clomn 获取 ratio def get_row_column_ratio(row, column ,district_dict, area_dict ): key = str(row) +"," + str(column) if area_dict.__contains__("全部区域"): return area_dict["全部区域"] if district_dict.__contains__(key): name = district_dict[key] ratio = area_dict[name] return float(ratio) return float(1.0) def updateRiverFile(model_name,predictiondata): flag = check_rain_param(predictiondata) if flag == "true": rain_ratio = float(predictiondata["rain"]["ratio"]) rain_base_year = predictiondata["rain"]["base_year"] if rain_base_year=='4': rain_base_year="1" river_ratio= float(predictiondata["river"]["ratio"]) area= predictiondata["mine"]["area"] # ws = base.predictParamModel + rain_base_year ws="" start_time = predictiondata["start_time"] end_time = predictiondata["end_time"] count = ModelPeriod.get_months_in_range_count(start_time, end_time) if count==12: ws= base.predictParamModel + rain_base_year else: ws= base.predictModel60 baseMdoel = flopy.modflow.Modflow.load("modflow.nam", model_ws= ws, exe_name="mf2005", verbose=True, version="mf2005", check=False) update_model_ws = base.model_dir + model_name updateMdoel = flopy.modflow.Modflow.load("modflow.nam", model_ws= update_model_ws, exe_name="mf2005", verbose=True, version="mf2005", check=False) district_dict = get_distric_dict() area_dict = get_area_dict(area) lrcq = {} periods =CalHead.get_model_period(model_name) periods_len= len(periods) for per in range(periods_len): wel = [] array2d = [] wel = baseMdoel.wel.stress_period_data.__getitem__(kper = per) wel_len = len(wel) #侧向边界 for i in range (0,453): wel[i][3] = wel[i][3] * rain_ratio #河流 for i in range(453, 562): wel[i][3] = wel[i][3] * river_ratio #抽水井 for i in range(562,wel_len): r = (float) (get_row_column_ratio(wel[i][1], wel[i][2], district_dict, area_dict)) wel[i][3] = wel[i][3] * r #重置数组 for Layer, Row, Column, Q in wel: array = [Layer, Row, Column, Q] array2d.append(array) #追加井 flex_data= getFlexdata(model_name) print("==============") print(flex_data) for i in range(len(flex_data)): array2d.append(flex_data[i]) lrcq[per] = array2d flopy.modflow.ModflowWel(updateMdoel, ipakcb= baseMdoel.wel.ipakcb, dtype=baseMdoel.wel.dtype, options=baseMdoel.wel.options, stress_period_data=lrcq) updateMdoel.write_input() else: print("Well--River文件无需修改！") #追加复杂原汇项信息 def getFlexdata(model_name): welldata="" well_path = base.model_dir + model_name +"\\pump_well.json" data=[] if os.path.exists(well_path): with open(well_path,encoding='utf-8') as f: welldata = json.load(f) wel= welldata["well"] for i in range (len(wel)): layer = int (wel[i]['layer'])-1 row= int(wel[i]['row'])-1 column = int(wel[i]['column'])-1 v = float(wel[i]['value']) arr = [layer,row, column, v] data.append(arr) return data def updateRchFile(model_name,predictiondata): flag = check_rain_param(predictiondata) if flag == "true": #丰水年 枯水年 base_year = predictiondata["rain"]["base_year"] if base_year =="1" or base_year =="2" or base_year =="3": updateRchBaseYear(model_name,predictiondata) elif base_year =="4": #接入传感器数据 updateRchRealData(model_name,predictiondata) else: print("Rch文件无需修改！") #更新平水年 丰水年 枯水年数据 def updateRchBaseYear(model_name,predictiondata): #丰水年 枯水年 base_year = predictiondata["rain"]["base_year"] ratio= float(predictiondata["rain"]["ratio"]) #数据来源的模型文件夹 base_ws="" start_time = predictiondata["start_time"] end_time = predictiondata["end_time"] count = ModelPeriod.get_months_in_range_count(start_time, end_time) if count==12: base_ws= base.predictParamModel + base_year else: base_ws= base.predictModel60 baseMdoel = flopy.modflow.Modflow.load("modflow.nam", model_ws= base_ws, exe_name="mf2005", verbose=True, version="mf2005", check=False) update_model_ws = base.model_dir + model_name updateMdoel = flopy.modflow.Modflow.load("modflow.nam", model_ws= update_model_ws, exe_name="mf2005", verbose=True, version="mf2005", check=False) periods =CalHead.get_model_period(model_name) periods_len= len(periods) for per in range(periods_len): item = baseMdoel.rch.rech.__getitem__(kper = per) array2d = item.get_value() array2d_len = len(array2d) for i in range(array2d_len): array_len = len(array2d[i]) for j in range(array_len): if str(base.area_array[i][j]) != '-9999': array2d[i][j] = array2d[i][j] * ratio updateMdoel.rch.rech.__setitem__(key = per, value=array2d) rch = flopy.modflow.ModflowRch(updateMdoel, nrchop=baseMdoel.rch.nrchop, ipakcb=baseMdoel.rch.ipakcb, rech=updateMdoel.rch.rech, irch =baseMdoel.rch.irch) rch.write_file(check=False) #接入传感器数据 def updateRchRealData(model_name,predictiondata): # 取平水年的rch文件 base_year = "1" ratio= float(predictiondata["rain"]["ratio"]) if not ratio: ratio = 1 excel_data = OpenExcel.read_jyl_excel(model_name) array_data =[] if not excel_data: print("降水模板数据未获取到！") return "降水模板数据未获取到！" for i in range(1,len(excel_data)): temp =[] for j in range(1,len(excel_data[i])): data = round(float(excel_data[i][j]),8) temp.append(data) array_data.append(temp) base_ws="" start_time = predictiondata["start_time"] end_time = predictiondata["end_time"] count = ModelPeriod.get_months_in_range_count(start_time, end_time) if count==12: base_ws= base.predictParamModel + base_year else: base_ws= base.predictModel60 baseMdoel = flopy.modflow.Modflow.load("modflow.nam", model_ws= base_ws, exe_name="mf2005", verbose=True, version="mf2005", check=False) update_model_ws = base.model_dir + model_name updateMdoel = flopy.modflow.Modflow.load("modflow.nam", model_ws= update_model_ws, exe_name="mf2005", verbose=True, version="mf2005", check=False) #降雨分区 area_dicts = base.getAreas() periods =CalHead.get_model_period(model_name) periods_len= len(periods) for per in range(periods_len): #16个分区某个周期的 数据 下标(0,15) 和area_dicts字典降雨分区的key对应 water_per_data = get_Array2_column(array_data,per) item = baseMdoel.rch.rech.__getitem__(kper = per) array2d = item.get_value() for key in area_dicts: tuples= area_dicts[key] values = water_per_data[int(key)] for i in range(len(tuples)): x = tuples[i][0] y = tuples[i][1] array2d[x][y]= values*ratio updateMdoel.rch.rech.__setitem__(key = per, value=array2d) rch = flopy.modflow.ModflowRch(updateMdoel, nrchop=baseMdoel.rch.nrchop, ipakcb=baseMdoel.rch.ipakcb, rech=updateMdoel.rch.rech, irch =baseMdoel.rch.irch) rch.write_file(check=False) print("降水分区数据更新完毕!") return "降水分区数据更新完毕！" #获取列数据 def get_Array2_column(array_data,column): arr = np.array(array_data) column_data = arr[:, column] return column_data def check_rain_param(predictiondata): rain = predictiondata["rain"] if not rain: print("Rch预测参数为空，无需要修改") return "false" base_year = predictiondata["rain"]["base_year"] if not base_year : print(" Rch : base_year预测参数为空，无需要修改") return "false" ratio= predictiondata["rain"]["ratio"] if not ratio : print(" Rch : ratio预测参数为空，无需要修改") return "false" return "true" #更新bas6文件 初始水头信息 def updateBase6File(model_name,predictdata): 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) #初始水头 if "initHeader" not in predictdata: print("=============has no initHeader=============") return False init_header = predictdata["initHeader"] dir = base.model_dir + init_header + "\\modflow.head" head = bf.HeadFile(dir) alldata = head.get_alldata() last_index = len(alldata)-1 last_array3= alldata[last_index] strt = ml.bas6.strt # strs = ml.bas6.strt.__getitem__(2) # print(strs.get_value()) strt.__setitem__(0,last_array3[0]) strt.__setitem__(1,last_array3[1]) strt.__setitem__(2,last_array3[2]) mfBase6 = flopy.modflow.ModflowBas( ml, strt= strt, ibound=ml.bas6.ibound, hnoflo=ml.bas6.hnoflo, extension="bas6",) mfBase6.write_file(check=False) return True #修改dis 文件 def updateDisFile(model_name, per): 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) mfDis = flopy.modflow.ModflowDis( ml, nlay=ml.dis.nlay, nrow=ml.dis.nrow, ncol=ml.dis.ncol, nper=per, delr=ml.dis.delr, delc=ml.dis.delc, top=ml.dis.top, botm=ml.dis.botm, perlen=ml.dis.perlen, nstp=ml.dis.nstp, tsmult=ml.dis.tsmult, steady=ml.dis.steady, itmuni=ml.dis.itmuni, lenuni=ml.dis.lenuni, extension="dis") mfDis.write_file(check=False)