天气学诊断实习六计算水汽流函数和水汽势函数

最新推荐文章于 2023-02-15 16:22:38 发布
mgl666
最新推荐文章于 2023-02-15 16:22:38 发布
阅读量4.4k
点赞数 5
分类专栏：天气学诊断实习文章标签： shell 算法 linux
本文链接：https://blog.csdn.net/m0_51153866/article/details/121619482
版权
天气学诊断实习专栏收录该内容
6 篇文章 62 订阅
订阅专栏
这段代码实现了计算水汽流函数和水汽势函数，并用Python的matplotlib库进行图像绘制。算法基于利布曼法和四边差分法，涉及气象学中的水汽通量、散度和涡度概念，用于分析特定时次的500hPa水平分布情况。
摘要由CSDN通过智能技术生成
一、实习目的:

熟悉水汽流函数和水汽势函数程序的实际编程计算。
二、实习内容:

编制计算水汽流函数和水汽势函数的程序，并且绘制出两个时次25日20时，26日20时的水汽流函数和水汽势函数分布（850hPa,500hPa)
三、算法原理

势函数初值： 0
流函数初值：
在这里插入图片描述
利布曼法：
在这里插入图片描述
四、代码实现

# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
import math
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
from pylab import *                                 #支持中文
mpl.rcParams['font.sans-serif'] = ['SimHei']

def micaps(a):
    data=np.zeros((29,53))
    for i in range(0,29):
        data[i,0:10]=a[i*6,0:10]
        data[i,10:20]=a[i*6+1,0:10]
        data[i,20:30]=a[i*6+2,0:10]
        data[i,30:40]=a[i*6+3,0:10]
        data[i,40:50]=a[i*6+4,0:10]
        data[i,50:53]=a[i*6+5,0:3]
    return data

##角度转弧度
def hd(x):
    a=math.pi/180*x
    return a

def sd(u,v,leftlon, rightlon, lowerlat, upperlat,f):
    a=6371000
    b=hd(f)
    n1=len(u[:,0])-1
    n2=len(u[0,:])-1
    N=int((upperlat-lowerlat)/f+1)
    lat=np.linspace(lowerlat, upperlat,N)
    lat= lat[::-1]
    data=np.zeros((29,53))
    ##四边差分
    for i in range(1,n1):
        for j in range (1,n2):
            data[0,j]=(1/2/a)*((u[1,j+1]-u[0,j-1])/(math.cos(hd(lat[0]))*b)+(v[1,j]-v[0,j])/b-2*v[0,j]*math.tan(hd(lat[0])))
            data[n1,j]=(1/2/a)*((u[n1,j+1]-u[n1,j-1])/(math.cos(hd(lat[n1]))*b)+(v[n1-1,j]-v[n1,j])/b-2*v[n1,j]*math.tan(hd(lat[n1])))
            data[i,0]=(1/2/a)*((u[i,1]-u[i,0])/(math.cos(hd(lat[i]))*b)+(v[i+1,1]-v[i-1,0])/b-2*v[i,0]*math.tan(hd(lat[i])))
            data[i,n2]=(1/2/a)*((u[i,n2-1]-u[i,n2])/(math.cos(hd(lat[i]))*b)+(v[i+1,n2]-v[i-1,n2])/b-2*v[i,n2]*math.tan(hd(lat[i])))
    ##中间部分差分
    for i in range(1,n1):
        for j in range (1,n2):
            data[i,j]=(1/2/a)*((u[i,j+1]-u[i,j-1])/(math.cos(hd(lat[i]))*b)+(v[i+1,j]-v[i-1,j])/b-2*v[i,j]*math.tan(hd(lat[i])))
    ##四角差分
    data[0,0]=(1/2/a)*((u[0,1]-u[0,0])/(math.cos(hd(lat[0]))*b)+(v[1,0]-v[0,0])/b-2*v[0,0]*math.tan(hd(lat[0])))
    data[0,n2]=(1/2/a)*((u[0,n2-1]-u[0,n2])/(math.cos(hd(lat[0]))*b)+(v[0,n2-1]-v[0,n2])/b-2*v[0,n2]*math.tan(hd(lat[0])))
    data[n1,0]=(1/2/a)*((u[n1,1]-u[n1,0])/(math.cos(hd(lat[n1]))*b)+(v[n1-1,0]-v[n1,0])/b-2*v[n1,0]*math.tan(hd(lat[n1])))
    data[n1,n2]=(1/2/a)*((u[n1-1,n2]-u[n1,n2])/(math.cos(hd(lat[n1]))*b)+(v[n1-1,n2]-v[n1,n2])/b-2*v[n1,n1]*math.tan(hd(lat[n1])))
    return data

def wd(u,v,leftlon, rightlon, lowerlat, upperlat,f):
    a=6371000
    b=hd(f)
    n1=len(u[:,0])-1
    n2=len(u[0,:])-1
    N=int((upperlat-lowerlat)/f+1)
    lat=np.linspace(lowerlat, upperlat,N)
    lat= lat[::-1]
    data=np.zeros((29,53))
    ##四边差分
    for i in range(1,n1):
        for j in range (1,n2):
            data[0,j]=(1/2/a)*((v[1,j]-v[0,j])/(math.cos(hd(lat[0]))*b)-(u[0,j+1]-u[0,j-1])/b+2*u[0,j]*math.tan(hd(lat[0])))
            data[n1,j]=(1/2/a)*((v[n1-1,j]-v[n1,j])/(math.cos(hd(lat[n1]))*b)-(u[n1,j+1]-u[n1,j-1])/b+2*u[n1,j]*math.tan(hd(lat[n1])))
            data[i,0]=(1/2/a)*((v[i+1,0]-v[i-1,0])/(math.cos(hd(lat[i]))*b)-(u[i,1]-u[i,0])/b+2*u[i,0]*math.tan(hd(lat[i])))
            data[i,n2]=(1/2/a)*((v[i+1,n2]-v[i-1,n2])/(math.cos(hd(lat[i]))*b)-(u[i,n2]-u[i,n2-1])/b+2*u[i,n2]*math.tan(hd(lat[i])))
    ##中间部分差分
    for i in range(1,n1):
        for j in range (1,n2):
            data[i,j]=(1/2/a)*((v[i+1,j]-v[i-1,j])/(math.cos(hd(lat[i]))*b)-(u[i,j+1]-u[i,j-1])/b+2*u[i,j]*math.tan(hd(lat[i])))
    ##四角差分
    data[0,0]=(1/2/a)*((v[1,0]-v[0,0])/(math.cos(hd(lat[0]))*b)-(u[0,1]-u[0,0])/b+2*u[0,0]*math.tan(hd(lat[0])))
    data[0,n2]=(1/2/a)*((v[1,n2]-v[0,n2])/(math.cos(hd(lat[0]))*b)-(u[0,n2-1]-u[0,n2])/b+2*u[0,n2]*math.tan(hd(lat[0])))
    data[n1,0]=(1/2/a)*((v[n1-1,0]-v[n1,0])/(math.cos(hd(lat[n1]))*b)-(u[n1,0]-u[n1,1])/b+2*u[n1,0]*math.tan(hd(lat[n1])))
    data[n1,n2]=(1/2/a)*((v[n1-1,n2]-v[n1,n2])/(math.cos(hd(lat[n1]))*b)-(u[n1,n2-1]-u[n1,n2])/b+2*u[n1,n1]*math.tan(hd(lat[n1])))
    return data

def shi(D,leftlon, rightlon, lowerlat, upperlat,f):
    kapa=np.zeros((29,53))
    R=np.zeros((29,53))
    a=6371000
    b=hd(f)
    N=int((upperlat-lowerlat)/f+1)
    lat=np.linspace(lowerlat, upperlat,N)
    lat= lat[::-1]
    for k in range(0,500):                  ##迭代次数
        t=kapa[1,1]
        #中间差分迭代部分
        for i in range(1,28):
            for j in range(1,52):
                R[i,j]=(kapa[i+1,j]+kapa[i-1,j])/((a*math.cos(hd(lat[i]))*b)**2)+(kapa[i,j+1]+kapa[i,j-1])/((a*b)**2)+D[i,j]-(2/((a*math.cos(hd(lat[i]))*b)**2)+2/(a*b)**2)*kapa[i,j]
                kapa[i,j]=kapa[i,j]+R[i,j]/2/(1/((a*math.cos(hd(lat[i]))*b)**2)+1/(a*b)**2)             
     
        if (k>1) & (R[1,1]<pow(10,-10)):
            print(k)
            return kapa
    return kapa
##流函数边界修正
def liu_bj(data,h,leftlon, rightlon, lowerlat, upperlat,f):
    n1=len(data[:,0])
    n2=len(data[0,:])
    N=int((upperlat-lowerlat)/f+1)
    lat=np.linspace(lowerlat, upperlat,N)
    lat= lat[::-1]
    for i in range(0,n1):
        for j in range (0,n2):
            if i==0 or i==n1-1 or j==0 or j==n2-1 :
                data[i,j]=h[i,j]/(2*7.19*10e-5*math.sin(hd(lat[i])))
    return data
         
##流函数
def liu(D,leftlon, rightlon, lowerlat, upperlat,f):
    kapa=np.zeros((29,53))
    kapa=liu_bj(kapa,h,leftlon, rightlon, lowerlat, upperlat,f)
    R=np.zeros((29,53))
    a=6371000
    b=hd(f)
    N=int((upperlat-lowerlat)/f+1)
    lat=np.linspace(lowerlat, upperlat,N)
    lat= lat[::-1]
    for k in range(0,500):                  ##迭代次数
        t=kapa[1,1]
        #中间差分迭代部分
        for i in range(1,28):
            for j in range(1,52):
                R[i,j]=(kapa[i+1,j]+kapa[i-1,j])/((a*math.cos(hd(lat[i]))*b)**2)+(kapa[i,j+1]+kapa[i,j-1])/((a*b)**2)+D[i,j]-(2/((a*math.cos(hd(lat[i]))*b)**2)+2/(a*b)**2)*kapa[i,j]
                kapa[i,j]=kapa[i,j]+R[i,j]/2/(1/((a*math.cos(hd(lat[i]))*b)**2)+1/(a*b)**2)             
        if (k>1) & (R[1,1]<pow(10,-10)):
            print(k)
            return kapa

    return kapa

##画图函数
def draw(data,leftlon, rightlon, lowerlat, upperlat,a,name):
    lon=np.arange(leftlon, rightlon+0.01,a)
    lat=np.arange(lowerlat, upperlat+0.01,a)
    data=data[::-1, :]##坐标反转
    #建立画布
    proj = ccrs.PlateCarree() # 设置投影
    fig, f2_ax1 = plt.subplots(figsize=(15,15), subplot_kw=dict(projection=proj))
    leftlon, rightlon, lowerlat, upperlat = (leftlon, rightlon, lowerlat, upperlat)
    #绘制
    data1=f2_ax1.contourf(lon,lat,data)
    # data2=f2_ax1.contour(lon,lat,data,colors='k', linewidths=1, linestyles='solid',levels=np.linspace(-32,60,23))
    # plt.clabel(data2,fontsize=10,colors='r',fmt='%.2f')
    #在画布的绝对坐标建立子图
    f2_ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
    #海岸线，50m精度
    f2_ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
    #湖泊数据
    f2_ax1.add_feature(cfeature.LAKES, alpha=0.5)
    #以下6条语句是定义地理坐标标签格式
    f2_ax1.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())
    f2_ax1.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())
    lon_formatter = cticker.LongitudeFormatter()
    lat_formatter = cticker.LatitudeFormatter()
    f2_ax1.xaxis.set_major_formatter(lon_formatter)
    f2_ax1.yaxis.set_major_formatter(lat_formatter)
    f2_ax1.set_title(name,loc='center',fontsize =20)
    # shrink 控制 colorbar 长度，pad 控制colorbar和图的距离
    plt.rcParams['axes.unicode_minus'] = False##负号显示问题
    plt.colorbar(data1, shrink=0.3, pad=0.02)#orientation='horizontal'位置参数
uv= pd.read_table(r'E:\实习\天气学诊断\micaps 11层\micaps 11层\uv\500\13052520.000',header=None,skiprows=3,sep='\s+')
uv= np.array(uv)
u=micaps(uv[0:174,:])
v=micaps(uv[174:,:])
h= pd.read_table(r'E:\实习\天气学诊断\micaps 11层\micaps 11层\height\500\13052520.000',header=None,skiprows=4,sep='\s+')
h= np.array(h)
h=micaps(h)
a= pd.read_table(r'E:\实习\天气学诊断\micaps 11层\micaps 11层\temper\500\13052520.000',header=None,skiprows=4,sep='\s+')
a= np.array(a)
a=micaps(a)
b= pd.read_table(r'E:\实习\天气学诊断\micaps 11层\micaps 11层\t-td\500\13052520.000',header=None,skiprows=4,sep='\s+')
b= np.array(b)
b=micaps(b)
td=a-b

e=np.zeros((29,53))
for i in range (0,29):
        for j in range(0,53):
            e[i,j]=6.1078*math.exp(17.26*td[i,j]/(273.16+td[i,j]-35.86))
q=622*e/(500-0.378*e)
##水汽通量的计算
qtl_500_u=1/9.8*u*q
qtl_500_v=1/9.8*v*q
    
##水汽通量散度涡度的计算
qtlsd_500=sd(qtl_500_u,qtl_500_v,30,160,10,80,2.5)
qtlwd_500=wd(qtl_500_u,qtl_500_v,30,160,10,80,2.5)

kapa=shi(qtlsd_500,30,160,10,80,2.5)
phi=liu(qtlwd_500,30,160,10,80,2.5)
draw(kapa,30,160,10,80,2.5,'25日20时 500hPa势函数')  
draw(phi,30,160,10,80,2.5,'25日20时 500hPa流函数')    
# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
import math
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
from pylab import *                                 #支持中文
mpl.rcParams['font.sans-serif'] = ['SimHei']

def micaps(a):
    data=np.zeros((29,53))
    for i in range(0,29):
        data[i,0:10]=a[i*6,0:10]
        data[i,10:20]=a[i*6+1,0:10]
        data[i,20:30]=a[i*6+2,0:10]
        data[i,30:40]=a[i*6+3,0:10]
        data[i,40:50]=a[i*6+4,0:10]
        data[i,50:53]=a[i*6+5,0:3]
    return data

##角度转弧度
def hd(x):
    a=math.pi/180*x
    return a

def sd(u,v,leftlon, rightlon, lowerlat, upperlat,f):
    a=6371000
    b=hd(f)
    n1=len(u[:,0])-1
    n2=len(u[0,:])-1
    N=int((upperlat-lowerlat)/f+1)
    lat=np.linspace(lowerlat, upperlat,N)
    lat= lat[::-1]
    data=np.zeros((29,53))
    ##四边差分
    for i in range(1,n1):
        for j in range (1,n2):
            data[0,j]=(1/2/a)*((u[1,j+1]-u[0,j-1])/(math.cos(hd(lat[0]))*b)+(v[1,j]-v[0,j])/b-2*v[0,j]*math.tan(hd(lat[0])))
            data[n1,j]=(1/2/a)*((u[n1,j+1]-u[n1,j-1])/(math.cos(hd(lat[n1]))*b)+(v[n1-1,j]-v[n1,j])/b-2*v[n1,j]*math.tan(hd(lat[n1])))
            data[i,0]=(1/2/a)*((u[i,1]-u[i,0])/(math.cos(hd(lat[i]))*b)+(v[i+1,1]-v[i-1,0])/b-2*v[i,0]*math.tan(hd(lat[i])))
            data[i,n2]=(1/2/a)*((u[i,n2-1]-u[i,n2])/(math.cos(hd(lat[i]))*b)+(v[i+1,n2]-v[i-1,n2])/b-2*v[i,n2]*math.tan(hd(lat[i])))
    ##中间部分差分
    for i in range(1,n1):
        for j in range (1,n2):
            data[i,j]=(1/2/a)*((u[i,j+1]-u[i,j-1])/(math.cos(hd(lat[i]))*b)+(v[i+1,j]-v[i-1,j])/b-2*v[i,j]*math.tan(hd(lat[i])))
    ##四角差分
    data[0,0]=(1/2/a)*((u[0,1]-u[0,0])/(math.cos(hd(lat[0]))*b)+(v[1,0]-v[0,0])/b-2*v[0,0]*math.tan(hd(lat[0])))
    data[0,n2]=(1/2/a)*((u[0,n2-1]-u[0,n2])/(math.cos(hd(lat[0]))*b)+(v[0,n2-1]-v[0,n2])/b-2*v[0,n2]*math.tan(hd(lat[0])))
    data[n1,0]=(1/2/a)*((u[n1,1]-u[n1,0])/(math.cos(hd(lat[n1]))*b)+(v[n1-1,0]-v[n1,0])/b-2*v[n1,0]*math.tan(hd(lat[n1])))
    data[n1,n2]=(1/2/a)*((u[n1-1,n2]-u[n1,n2])/(math.cos(hd(lat[n1]))*b)+(v[n1-1,n2]-v[n1,n2])/b-2*v[n1,n1]*math.tan(hd(lat[n1])))
    return data

def wd(u,v,leftlon, rightlon, lowerlat, upperlat,f):
    a=6371000
    b=hd(f)
    n1=len(u[:,0])-1
    n2=len(u[0,:])-1
    N=int((upperlat-lowerlat)/f+1)
    lat=np.linspace(lowerlat, upperlat,N)
    lat= lat[::-1]
    data=np.zeros((29,53))
    ##四边差分
    for i in range(1,n1):
        for j in range (1,n2):
            data[0,j]=(1/2/a)*((v[1,j]-v[0,j])/(math.cos(hd(lat[0]))*b)-(u[0,j+1]-u[0,j-1])/b+2*u[0,j]*math.tan(hd(lat[0])))
            data[n1,j]=(1/2/a)*((v[n1-1,j]-v[n1,j])/(math.cos(hd(lat[n1]))*b)-(u[n1,j+1]-u[n1,j-1])/b+2*u[n1,j]*math.tan(hd(lat[n1])))
            data[i,0]=(1/2/a)*((v[i+1,0]-v[i-1,0])/(math.cos(hd(lat[i]))*b)-(u[i,1]-u[i,0])/b+2*u[i,0]*math.tan(hd(lat[i])))
            data[i,n2]=(1/2/a)*((v[i+1,n2]-v[i-1,n2])/(math.cos(hd(lat[i]))*b)-(u[i,n2]-u[i,n2-1])/b+2*u[i,n2]*math.tan(hd(lat[i])))
    ##中间部分差分
    for i in range(1,n1):
        for j in range (1,n2):
            data[i,j]=(1/2/a)*((v[i+1,j]-v[i-1,j])/(math.cos(hd(lat[i]))*b)-(u[i,j+1]-u[i,j-1])/b+2*u[i,j]*math.tan(hd(lat[i])))
    ##四角差分
    data[0,0]=(1/2/a)*((v[1,0]-v[0,0])/(math.cos(hd(lat[0]))*b)-(u[0,1]-u[0,0])/b+2*u[0,0]*math.tan(hd(lat[0])))
    data[0,n2]=(1/2/a)*((v[1,n2]-v[0,n2])/(math.cos(hd(lat[0]))*b)-(u[0,n2-1]-u[0,n2])/b+2*u[0,n2]*math.tan(hd(lat[0])))
    data[n1,0]=(1/2/a)*((v[n1-1,0]-v[n1,0])/(math.cos(hd(lat[n1]))*b)-(u[n1,0]-u[n1,1])/b+2*u[n1,0]*math.tan(hd(lat[n1])))
    data[n1,n2]=(1/2/a)*((v[n1-1,n2]-v[n1,n2])/(math.cos(hd(lat[n1]))*b)-(u[n1,n2-1]-u[n1,n2])/b+2*u[n1,n1]*math.tan(hd(lat[n1])))
    return data

def shi(D,leftlon, rightlon, lowerlat, upperlat,f):
    kapa=np.zeros((29,53))
    R=np.zeros((29,53))
    a=6371000
    b=hd(f)
    N=int((upperlat-lowerlat)/f+1)
    lat=np.linspace(lowerlat, upperlat,N)
    lat= lat[::-1]
    for k in range(0,500):                  ##迭代次数
        t=kapa[1,1]
        #中间差分迭代部分
        for i in range(1,28):
            for j in range(1,52):
                R[i,j]=(kapa[i+1,j]+kapa[i-1,j])/((a*math.cos(hd(lat[i]))*b)**2)+(kapa[i,j+1]+kapa[i,j-1])/((a*b)**2)+D[i,j]-(2/((a*math.cos(hd(lat[i]))*b)**2)+2/(a*b)**2)*kapa[i,j]
                kapa[i,j]=kapa[i,j]+R[i,j]/2/(1/((a*math.cos(hd(lat[i]))*b)**2)+1/(a*b)**2)             
(kapa[1,j]+kapa[0,j])/((2*a*math.cos(hd(lat[i]))*b)**2)+(kapa[0,j+1]+kapa[0,j-1])/(2*a*b)**2+D[0,j]-(2/((2*a*math.cos(hd(lat[i]))*b)**2)+2/(2*a*b)**2)*kapa[0,j]

        if (k>1) & (R[1,1]<pow(10,-10)):
            print(k)
            return kapa
    return kapa
##流函数边界修正
def liu_bj(data,h,leftlon, rightlon, lowerlat, upperlat,f):
    n1=len(data[:,0])
    n2=len(data[0,:])
    N=int((upperlat-lowerlat)/f+1)
    lat=np.linspace(lowerlat, upperlat,N)
    lat= lat[::-1]
    for i in range(0,n1):
        for j in range (0,n2):
            if i==0 or i==n1-1 or j==0 or j==n2-1 :
                data[i,j]=h[i,j]/(2*7.19*10e-5*math.sin(hd(lat[i])))
    return data
         
##流函数
def liu(D,leftlon, rightlon, lowerlat, upperlat,f):
    kapa=np.zeros((29,53))
    kapa=liu_bj(kapa,h,leftlon, rightlon, lowerlat, upperlat,f)
    R=np.zeros((29,53))
    a=6371000
    b=hd(f)
    N=int((upperlat-lowerlat)/f+1)
    lat=np.linspace(lowerlat, upperlat,N)
    lat= lat[::-1]
    for k in range(0,500):                  ##迭代次数
        t=kapa[1,1]
        #中间差分迭代部分
        for i in range(1,28):
            for j in range(1,52):
                R[i,j]=(kapa[i+1,j]+kapa[i-1,j])/((a*math.cos(hd(lat[i]))*b)**2)+(kapa[i,j+1]+kapa[i,j-1])/((a*b)**2)+D[i,j]-(2/((a*math.cos(hd(lat[i]))*b)**2)+2/(a*b)**2)*kapa[i,j]
                kapa[i,j]=kapa[i,j]+R[i,j]/2/(1/((a*math.cos(hd(lat[i]))*b)**2)+1/(a*b)**2)             
        if (k>1) & (R[1,1]<pow(10,-10)):
            print(k)
            return kapa

    return kapa

##画图函数
def draw(data,leftlon, rightlon, lowerlat, upperlat,a,name):
    lon=np.arange(leftlon, rightlon+0.01,a)
    lat=np.arange(lowerlat, upperlat+0.01,a)
    data=data[::-1, :]##坐标反转
    #建立画布
    proj = ccrs.PlateCarree() # 设置投影
    fig, f2_ax1 = plt.subplots(figsize=(15,15), subplot_kw=dict(projection=proj))
    leftlon, rightlon, lowerlat, upperlat = (leftlon, rightlon, lowerlat, upperlat)
    #绘制
    data1=f2_ax1.contourf(lon,lat,data)
    # data2=f2_ax1.contour(lon,lat,data,colors='k', linewidths=1, linestyles='solid',levels=np.linspace(-32,60,23))
    # plt.clabel(data2,fontsize=10,colors='r',fmt='%.2f')
    #在画布的绝对坐标建立子图
    f2_ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
    #海岸线，50m精度
    f2_ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
    #湖泊数据
    f2_ax1.add_feature(cfeature.LAKES, alpha=0.5)
    #以下6条语句是定义地理坐标标签格式
    f2_ax1.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())
    f2_ax1.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())
    lon_formatter = cticker.LongitudeFormatter()
    lat_formatter = cticker.LatitudeFormatter()
    f2_ax1.xaxis.set_major_formatter(lon_formatter)
    f2_ax1.yaxis.set_major_formatter(lat_formatter)
    f2_ax1.set_title(name,loc='center',fontsize =20)
    # shrink 控制 colorbar 长度，pad 控制colorbar和图的距离
    plt.rcParams['axes.unicode_minus'] = False##负号显示问题
    plt.colorbar(data1, shrink=0.3, pad=0.02)#orientation='horizontal'位置参数
uv= pd.read_table(r'E:\实习\天气学诊断\micaps 11层\micaps 11层\uv\500\13052520.000',header=None,skiprows=3,sep='\s+')
uv= np.array(uv)
u=micaps(uv[0:174,:])
v=micaps(uv[174:,:])
h= pd.read_table(r'E:\实习\天气学诊断\micaps 11层\micaps 11层\height\500\13052520.000',header=None,skiprows=4,sep='\s+')
h= np.array(h)
h=micaps(h)
a= pd.read_table(r'E:\实习\天气学诊断\micaps 11层\micaps 11层\temper\500\13052520.000',header=None,skiprows=4,sep='\s+')
a= np.array(a)
a=micaps(a)
b= pd.read_table(r'E:\实习\天气学诊断\micaps 11层\micaps 11层\t-td\500\13052520.000',header=None,skiprows=4,sep='\s+')
b= np.array(b)
b=micaps(b)
td=a-b

e=np.zeros((29,53))
for i in range (0,29):
        for j in range(0,53):
            e[i,j]=6.1078*math.exp(17.26*td[i,j]/(273.16+td[i,j]-35.86))
q=622*e/(500-0.378*e)
##水汽通量的计算
qtl_500_u=1/9.8*u*q
qtl_500_v=1/9.8*v*q
    
##水汽通量散度涡度的计算
qtlsd_500=sd(qtl_500_u,qtl_500_v,30,160,10,80,2.5)
qtlwd_500=wd(qtl_500_u,qtl_500_v,30,160,10,80,2.5)

kapa=shi(qtlsd_500,30,160,10,80,2.5)
phi=liu(qtlwd_500,30,160,10,80,2.5)
draw(kapa,30,160,10,80,2.5,'25日20时 500hPa势函数')  
draw(phi,30,160,10,80,2.5,'25日20时 500hPa流函数')