天气学诊断实习三 计算涡度、散度、涡度平流和温度平流

天气学诊断实习三 计算涡度、散度、涡度平流和温度平流 一、实习目的:

熟悉使用差分方法解决方程中的导数、偏导数，计算涡度、散度、涡度平流和温度平流的程序

二、实习内容:

编制计算涡度、散度、涡度平流和温度平流的程序，并且绘制出两个时次25日20时，26日20时的涡度、散度、涡度平流和温度平流分布（850hPa,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 hd(x):
    a=math.pi/180*x
    return a
##micaps读数据函数
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 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 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 pl(T,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]=u[0,j]*(T[0,j+1]-T[0,j-1])/(math.cos(hd(lat[0]))*b*a)+v[0,j]*(T[1,j]-T[0,j])/(a*b)
            data[n1,j]=u[n1,j]*(T[n1,j+1]-T[n1,j-1])/(math.cos(hd(lat[n1]))*b*a)+v[n1,j]*(T[n1-1,j]-T[n1,j])/(a*b)
            data[i,0]=u[i,0]*(T[i,1]-T[i,0])/(math.cos(hd(lat[i]))*b*a)+v[i,0]*(T[i+1,j]-T[i-1,j])/(a*b)
            data[i,n2]=u[i,n2]*(T[i,n2-1]-T[0,n2])/(math.cos(hd(lat[i]))*b*a)+v[i,n2]*(T[i+1,n2]-T[i-1,n2])/(a*b)
    ##中间部分差分
    for i in range(1,n1):
        for j in range (1,n2):
            data[i,j]=u[i,j]*(T[i,j+1]-T[i,j-1])/(math.cos(hd(lat[i]))*b*a)+v[i,j]*(T[i+1,j]-T[i-1,j])/(a*b)
    ##四角差分    
    data[0,0]=u[0,0]*(T[0,1]-T[0,0])/(math.cos(hd(lat[0]))*b*a)+v[0,0]*(T[1,0]-T[0,0])/(a*b)
    data[0,n2]=u[0,n2]*(T[0,n2-1]-T[0,n2])/(math.cos(hd(lat[0]))*b*a)+v[0,n2]*(T[1,n2]-T[0,n2])/(a*b)
    data[n1,0]=u[n1,0]*(T[n1,1]-T[n1,0])/(math.cos(hd(lat[n1]))*b*a)+v[n1,0]*(T[n1-1,0]-T[n1,0])/(a*b)
    data[n1,n2]=u[n1,n2]*(T[n1,n2-1]-T[n1,n2])/(math.cos(hd(lat[n1]))*b*a)+v[n1,n2]*(T[n1-1,n2]-T[n1,n2])/(a*b)
    return data

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'位置参数
##读数据
uvdata500_25=pd.read_table(r'C:Users马冠龙Desktopmicaps 11层micaps 11层uv50013052520.000',header=None,skiprows=3,sep='s+')
uvdata500_25=np.array(uvdata500_25)
udata500_25=micaps(uvdata500_25[0:174,:])
vdata500_25=micaps(uvdata500_25[174:,:])
uvdata500_26=pd.read_table(r'C:Users马冠龙Desktopmicaps 11层micaps 11层uv50013052620.000',header=None,skiprows=3,sep='s+')
uvdata500_26=np.array(uvdata500_26)
udata500_26=micaps(uvdata500_26[0:174,:])
vdata500_26=micaps(uvdata500_26[174:,:])
uvdata850_25=pd.read_table(r'C:Users马冠龙Desktopmicaps 11层micaps 11层uv85013052520.000',header=None,skiprows=3,sep='s+')
uvdata850_25=np.array(uvdata850_25)
udata850_25=micaps(uvdata850_25[0:174,:])
vdata850_25=micaps(uvdata850_25[174:,:])
uvdata850_26=pd.read_table(r'C:Users马冠龙Desktopmicaps 11层micaps 11层uv85013052620.000',header=None,skiprows=3,sep='s+')
uvdata850_26=np.array(uvdata850_26)
udata850_26=micaps(uvdata850_26[0:174,:])
vdata850_26=micaps(uvdata850_26[174:,:])
tdata500_25=pd.read_table(r'C:Users马冠龙Desktopmicaps 11层micaps 11层temper50013052520.000',header=None,skiprows=4,sep='s+')
tdata500_25=np.array(tdata500_25)
tdata500_25=micaps(tdata500_25)
tdata500_26=pd.read_table(r'C:Users马冠龙Desktopmicaps 11层micaps 11层temper50013052620.000',header=None,skiprows=4,sep='s+')
tdata500_26=np.array(tdata500_26)
tdata500_26=micaps(tdata500_26)
tdata850_25=pd.read_table(r'C:Users马冠龙Desktopmicaps 11层micaps 11层temper85013052520.000',header=None,skiprows=4,sep='s+')
tdata850_25=np.array(tdata850_25)
tdata850_25=micaps(tdata850_25)
tdata850_26=pd.read_table(r'C:Users马冠龙Desktopmicaps 11层micaps 11层temper85013052620.000',header=None,skiprows=4,sep='s+')
tdata850_26=np.array(tdata850_26)
tdata850_26=micaps(tdata850_26)
##涡度，散度，涡度平流，温度平流的计算
wd500_25=wd(udata500_25,vdata500_25,30,160,10,80,2.5)
sd500_25=sd(udata500_25,vdata500_25,30,160,10,80,2.5)
wdpl500_25=pl(wd500_25,udata500_25,vdata500_25,30,160,10,80,2.5)
tpl500_25=pl(tdata500_25,udata500_25,vdata500_25,30,160,10,80,2.5)
wd500_26=wd(udata500_26,vdata500_26,30,160,10,80,2.5)
sd500_26=sd(udata500_26,vdata500_26,30,160,10,80,2.5)
wdpl500_26=pl(wd500_26,udata500_26,vdata500_26,30,160,10,80,2.5)
tpl500_26=pl(tdata500_26,udata500_26,vdata500_26,30,160,10,80,2.5)
wd850_25=wd(udata850_25,vdata850_25,30,160,10,80,2.5)
sd850_25=sd(udata850_25,vdata850_25,30,160,10,80,2.5)
wdpl850_25=pl(wd850_25,udata850_25,vdata850_25,30,160,10,80,2.5)
tpl850_25=pl(tdata850_25,udata850_25,vdata850_25,30,160,10,80,2.5)
wd850_26=wd(udata850_26,vdata850_26,30,160,10,80,2.5)
sd850_26=sd(udata850_26,vdata850_26,30,160,10,80,2.5)
wdpl850_26=pl(wd850_26,udata850_26,vdata850_26,30,160,10,80,2.5)
tpl850_26=pl(tdata850_26,udata850_26,vdata850_26,30,160,10,80,2.5)
##画图
draw(wd500_25,30,160,10,80,2.5,'500hPa涡度图  2013年5月25日20时')
draw(wd500_26,30,160,10,80,2.5,'500hPa涡度图  2013年5月26日20时')
draw(wd850_25,30,160,10,80,2.5,'850hPa涡度图  2013年5月25日20时')
draw(wd850_26,30,160,10,80,2.5,'850hPa涡度图  2013年5月26日20时')
draw(sd500_25,30,160,10,80,2.5,'500hPa散度图  2013年5月25日20时')
draw(sd500_26,30,160,10,80,2.5,'500hPa散度图  2013年5月26日20时')
draw(sd850_25,30,160,10,80,2.5,'850hPa散度图  2013年5月25日20时')
draw(sd850_26,30,160,10,80,2.5,'850hPa散度图  2013年5月26日20时')
draw(wdpl500_25,30,160,10,80,2.5,'500hPa涡度平流图  2013年5月25日20时')
draw(wdpl500_26,30,160,10,80,2.5,'500hPa涡度平流图  2013年5月26日20时')
draw(wdpl850_25,30,160,10,80,2.5,'850hPa涡度平流图  2013年5月25日20时')
draw(tpl850_26,30,160,10,80,2.5,'850hPa涡度平流图  2013年5月26日20时')
draw(tpl500_25,30,160,10,80,2.5,'500hPa温度平流图  2013年5月25日20时')
draw(tpl500_26,30,160,10,80,2.5,'500hPa温度平流图  2013年5月26日20时')
draw(tpl850_25,30,160,10,80,2.5,'850hPa温度平流图  2013年5月25日20时')
draw(tpl850_26,30,160,10,80,2.5,'850hPa温度平流图  2013年5月26日20时')

五、实习结果

（图片顺序有点乱）