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
|
|
|
|
|
|
|
|
