import matplotlib.pyplot as plt from matplotlib import cm import numpy as np import math from docxtpl import DocxTemplate import flopy import flopy.utils.binaryfile as bf from matplotlib import cm, colors import os import time import Base as base import CalHead import WaterXBL #获取Z值 def get_flow_Z(model_name,per,layer): dir = base.model_dir + model_name + "\\modflow.head" mlist = CalHead.muiltyModelList() if model_name==base.not_allowed_model: dir = base.baseModel2 + "\\modflow.head" if model_name in mlist : dir = base.muiltyModel + model_name + "\\modflow.head" head = bf.HeadFile(dir) alldata = head.get_alldata() #获取周期数据 z = alldata[int(per+1)*3-1,int(layer),:,:] z[(z<=0)] = 0 return z #读取base6文件读取初始水位 def get_init_flowZ(model_name,layer): 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) strs = ml.bas6.strt.__getitem__(layer) z = strs.get_value() z[(z<=0)] = 0 #处理验证模型的无效边界 if model_name==base.not_allowed_model: arr = np.loadtxt(base.areapath, dtype=int) for i in range(len(arr)): for j in range(len(arr[i])): if arr[i][j] == -9999: z[i][j]=0 return z #获取流场图的打了等值线数组 def get_flow_levels(data): maxdata= np.max(data) levels=[] levels.append(0.1) lines = 10 while lines < maxdata: levels.append(lines) lines = lines + 10 levels.append(math.ceil(maxdata)) return levels #获取变幅比较小的等值线图 def get_flow_levels_small(min , max, line_count): step = round((max-min)/line_count,0) if step <= 0: step = 0.5 levels=[] levels.append(min) lines = min + step while lines < max: levels.append(lines) lines = lines + step levels.append(max) return levels #获取变幅的Z值 def get_bf_z(model_name,per1,per2): dir = base.model_dir + model_name + "\\modflow.head" mlist = CalHead.muiltyModelList() if model_name==base.not_allowed_model: dir = base.baseModel2 + "\\modflow.head" if model_name in mlist : dir = base.muiltyModel + model_name + "\\modflow.head" head = bf.HeadFile(dir) alldata = head.get_alldata() #变幅值 z1= alldata[int(per1+1)*3-1,0,:,:] z2 = alldata[int(per2+1)*3-1,0,:,:] res = np.subtract(z2,z1) return res #获取变幅的流场图 def getFlowFieldBF(model_name,per1,per2,outpath): res = get_bf_z(model_name,per1,per2) max = np.max(res) min = np.min(res) # 把无效区设置为9999 arr = np.loadtxt(base.areapath, dtype=int) for i in range(len(arr)): for j in range(len(arr[i])): if arr[i][j] == -9999: res[i][j]= 9999.0 levels= get_flow_levels_small(min,max,10) draw_flow(res,levels,outpath,cm.RdBu_r) str_levs=[] for item in levels: str_levs.append(str(round(item,2))) cols = get_gradient_color_list("RdBu_r",len(levels)) dicts ={"levels":str_levs,"colors":cols} xbl = WaterXBL.get_grd_storage(model_name,int(per1),int(per2)) dicts["xbl"]=xbl return dicts #获取变幅 埋深 def getFlowFieldDepth(model_name,per,layer,title,flow_types,out_path): Z= get_flow_Z(model_name,per,layer) dis = np.loadtxt(base.dis_top_path, dtype=str) areaMatrix = np.loadtxt(base.areapath, dtype=int) for i in range(len(areaMatrix)): for j in range(len(areaMatrix[i])): if areaMatrix[i][j] == -9999: Z[i][j] =0.0 else : Z[i][j] =float(dis[i][j])-Z[i][j] levels=[10,30,40,60,100,200,300,400,500,700,900] draw_flow(Z,levels,out_path,cm.RdBu_r) cols = get_gradient_color_list("RdBu_r",len(levels)) dicts ={"levels":levels,"colors":cols} return dicts def getWaterResFiled(model_name,per): pic = str(int(time.time())) +".png" outpath = base.flow_file + pic flow_field(model_name,per,0,"flow","online",outpath) return pic #流场图 def flow_field(model_name,per,layer,title,flow_types,out_path): Z= get_flow_Z(model_name,per,layer) levels = get_flow_levels(Z) draw_flow(Z,levels,out_path,cm.RdBu) cols = get_gradient_color_list("RdBu",len(levels)) dicts ={"levels":levels,"colors":cols} return dicts #初始流场图 def init_flow_field(model_name,layer,title,flow_types,out_path): #如果是验证模型直接读取base6的文件当做初始水头 if model_name == base.not_allowed_model: Z= get_init_flowZ(model_name,layer) levels = get_flow_levels(Z) draw_flow(Z,levels,out_path,cm.RdBu) else: pjson= CalHead.get_model_json(model_name) if "initHeader" in pjson: initHead = pjson["initHeader"] print(initHead) Z= get_flow_Z(initHead,0,0) levels = get_flow_levels(Z) draw_flow(Z,levels,out_path,cm.RdBu) else: return "#" print(out_path) return out_path #绘制流场 def draw_flow(Z,levels,out_path,colorbar): #开始绘制流场图 x = np.arange(0, 114, 1) y = np.arange(0, 104, 1) #网格 X, Y = np.meshgrid(x, y) #坐标轴翻转 左上角为起始 plt.gca().invert_yaxis() #等值线 C = plt.contour(X, Y, Z, levels=levels, linewidths=0.5, colors='white') #等值线样式 plt.clabel(C, inline=2,fmt='%.2f', fontsize=4,colors='black') #绘制等值线之间的填充 # cmap="RdBu_r" cm.RdBu_r plt.contourf(X, Y, Z,levels= levels,alpha = 0.75,cmap=colorbar) plt.axis("off") plt.colorbar().ax.set_visible(False) plt.xticks([]) plt.yticks([]) plt.savefig(out_path, dpi=300,transparent=True, bbox_inches='tight') plt.close('all') return out_path def get_gradient_color_list(m_color_name, m_num): m_color_list = [] m_color_map = plt.get_cmap(m_color_name, m_num) for m_i in range(m_num): m_color = tuple([int(_ * 256) for _ in list(m_color_map([m_i][0])[:-1])]) m_color_list.append(m_color) return m_color_list #导出模板 def exportReport(model_name,period): tpl = DocxTemplate(base.prefix + '地下水发布信息模板.docx') cont = archive_report_content(model_name,period) tpl.render(cont) save_path="" tpl.save(save_path) #报告信息 def archive_report_content(model_name,period): data1 = archive_grd_depth(model_name,period) data2= archive_grd_res(model_name,period) #初始流场 initpic = str(int(time.time())) +"1.png" outpath = base.flow_file + initpic flow_field(model_name,0,0,"初始流场信息", "online",outpath) currentpic = str(int(time.time())) +"2.png" outpath2 = base.flow_file + currentpic flow_field(model_name,int(period),0,"流场信息", "online",outpath2) content={ "py_avg_water":str(data1[0]), "py_m_water":str(data1[1]), "py_y_water":str(data1[2]), "sq_avg_water":str(data1[3]), "sq_m_water":str(data1[4]), "sq_y_water":str(data1[5]), "yq_avg_water":str(data1[6]), "yq_m_water":str(data1[7]), "yq_y_water":str(data1[8]), "w_m_res":str(data2[0]), "w_y_res":str(data2[1]), "flow1":"/xishan/xinshanFlow/"+initpic, "flow2":"/xishan/xinshanFlow/"+currentpic } if data1[1]>=0: content["py_m_water"]= "回升"+str(abs(data1[1])) else: content["py_m_water"]= "下降"+str(abs(data1[1])) if data1[2]>=0: content["py_y_water"]= "回升"+str(abs(data1[2])) else: content["py_y_water"]= "下降"+str(abs(data1[2])) if data1[4]>=0: content["sq_m_water"]= "回升"+str(abs(data1[4])) else: content["sq_m_water"]= "下降"+str(abs(data1[4])) if data1[5]>=0: content["sq_y_water"]= "回升"+str(abs(data1[5])) else: content["sq_y_water"]= "下降"+str(abs(data1[5])) if data1[7]>=0: content["yq_m_water"]= "回升"+str(abs(data1[7])) else: content["yq_m_water"]= "下降"+str(abs(data1[7])) if data1[8]>=0: content["yq_y_water"]= "回升"+str(abs(data1[8])) else: content["yq_y_water"]= "下降"+str(abs(data1[8])) if data2[1]>=0: content["w_y_res"]= "增加"+str(abs(data2[1])) else: content["w_y_res"]= "减少"+str(abs(data2[1])) return content #水资源量 def archive_grd_res(model_name,period): txt_path = base.model_dir + model_name + "\\water_res.txt" if not os.path.exists(txt_path): CalHead.run_zonebudget_res(model_name) monthdata = CalHead.water_res_month(model_name,txt_path,int(period)) monthdata2=[] if int(period) > 0: monthdata2 = CalHead.water_res_month(model_name,txt_path,int(period-1)) else: monthdata2 = monthdata water1 = monthdata[0] + monthdata[1]+ monthdata[2]- monthdata[3] water2 = monthdata2[0] + monthdata2[1] + monthdata2[2] - monthdata2[3] water1 = round(water1 ,4) wat= round(water1-water2 ,4) return [water1,wat] #生成月度报告的内容地下水深 def archive_grd_depth(model_name,period): model_dir = base.model_dir + model_name + "\\modflow.head" head = bf.HeadFile(model_dir) alldata = head.get_alldata() #获取月度数据 per = int(period) current_month_data=[] pre_month_data=[] pre_year_data=[] if per > 0: current_month_data = alldata[int(per+1)*3-1,0,:,:] pre_month_data= alldata[int(per)*3-1,0,:,:] else : current_month_data = alldata[int(per+1)*3-1,0,:,:] pre_month_data= current_month_data mpdict = getMPDict() YQSDict = getYQSDict() pingyuanArray = mpdict["1"] shanquArray = mpdict["2"] yqArray = YQSDict["1"] #高程坐标 gc_array= np.array(getTopDis()) depth_array = np.subtract(gc_array,current_month_data) py_data,py_data2 = 0,0 sq_data, sq_data2 = 0,0 yqs_data , yqs_data2 = 0,0 # 埋深 py_depth,sq_depth,yqs_depth =0,0,0 k,m,n =0,0,0 for item in pingyuanArray: i,j = item[0],item[1] if current_month_data[i,j]>0: py_depth += depth_array[i,j] py_data += current_month_data[i,j] k+=1 if pre_month_data[i,j]>0: py_data2 += pre_month_data[i,j] for item in shanquArray: i,j = item[0],item[1] if current_month_data[i,j]>0: sq_depth += depth_array[i,j] sq_data += current_month_data[i,j] m+=1 if pre_month_data[i,j]>0: sq_data2 += pre_month_data[i,j] for item in yqArray: i,j = item[0],item[1] if current_month_data[i,j]>0: yqs_depth += depth_array[i,j] yqs_data += current_month_data[i,j] n+=1 if pre_month_data[i,j]>0: yqs_data2 += pre_month_data[i,j] py_data = py_data/k sq_data = sq_data/m yqs_data= yqs_data/n py_data2 = py_data2/k sq_data2 = sq_data2/m yqs_data2= yqs_data2/n py_depth = py_depth/k sq_depth = sq_depth/m yqs_depth = yqs_depth/n result=[py_depth,py_data-py_data2,0, sq_depth,sq_data-sq_data2,0, yqs_depth,yqs_data-yqs_data2,0] for i in range(len(result)): result[i]= round(result[i],2) return result #获取研究区域 山区 平原区 玉泉山区的折线图 def getXs3LineChart(paths): head = bf.HeadFile(paths) alldata = head.get_alldata() months = int(len(alldata)/3) mpdict = getMPDict() YQSDict = getYQSDict() pingyuanArray = mpdict["1"] shanquArray = mpdict["2"] yqArray = YQSDict["1"] result1=[] result2=[] result3=[] for per in range(months): current_month = alldata[int(per+1)*3-1,0,:,:] yqs_data, py_data, sq_data = 0.0, 0.0, 0.0 M ,N,k = 0,0,0 for item in pingyuanArray: i,j = item[0],item[1] if current_month[i,j]>0: py_data += current_month[i,j] M+=1 py_data= round(py_data/M,2) result1.append(py_data) for item in shanquArray: i,j = item[0],item[1] if current_month[i,j]>0: sq_data += current_month[i,j] N+=1 sq_data= round(sq_data/N,2 ) result2.append(sq_data) for item in yqArray: i,j = item[0],item[1] if current_month[i,j]>0: yqs_data += current_month[i,j] k+=1 yqs_data= round(yqs_data/k,2 ) result3.append(yqs_data) dicts={"pyq":result1,"sq":result2,"yqs":result3} return dicts #山区/平原区的字典 #平原区=1 山区 =2 def getMPDict(): arr = np.loadtxt(base.xs_mp_path, dtype=int) dict ={} for i in range(len(arr)): for j in range(len(arr[i])): zb = str(arr[i][j]) if arr[i][j] == -9999: continue if zb not in dict: dict[zb] = [(i,j)] else: dict[zb].append((i,j)) return dict #玉泉山的字典 #玉泉山 = 1 def getYQSDict(): arr = np.loadtxt(base.yqs_path, dtype=int) dict ={} for i in range(len(arr)): for j in range(len(arr[i])): zb = str(arr[i][j]) if arr[i][j] != 1: continue if zb not in dict: dict[zb] = [(i,j)] else: dict[zb].append((i,j)) return dict #获取高程 def getTopDis(): arr = np.loadtxt(base.dis_top_path, dtype=str) float_array = np.array(arr).astype("float") return float_array