一、实习内容

  编制计算位温和相当位温的程序，并且绘制出两个时次25日20时，26日20时的位温和相当位温分布（850hPa,500hPa)

二、代码实现

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
import math
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 ww(t,P):
    theta=(t+273.16)*(1000/P)**0.286-273.16
    return theta
##相当位温计算函数,代入，输出数据摄氏度
def xdww(t,td,P):
    e=6.1078*math.exp(17.26*td/(273.16+td-35.86))
    q=622*e/(P-0.378*e)
    cp=1004*(1+0.85*q)
    Rd=287
    L=2260
    t=t+273.16
    td=td+273.16
    theta_e=t*math.exp((Rd/cp)*math.log(1000/(P-e))+L*q/cp/td+q/0.622*math.log(t/td))-273.16
    return theta_e

##画图函数
##data画图数据矩阵leftlon, rightlon, lowerlat经纬度范围,a数据经度,name图片名称
def draw1(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,levels=np.linspace(-32,60,24),extend='both')
    # 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'位置参数
    
def draw2(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,levels=np.linspace(-30,80,20),extend='both')
    # data2=f2_ax1.contour(lon,lat,data,colors='k', linewidths=1, linestyles='solid')
    # 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'位置参数   

##读数据
tdata500_25=pd.read_table(r'E:\实习\天气学诊断\micaps\micaps\t\500\13052520.000',header=None,skiprows=4,sep='\s+')
tdata500_25=np.array(tdata500_25)
tdata500_25=micaps(tdata500_25)
tdata500_26=pd.read_table(r'E:\实习\天气学诊断\micaps\micaps\t\500\13052620.000',header=None,skiprows=4,sep='\s+')
tdata500_26=np.array(tdata500_26)
tdata500_26=micaps(tdata500_26)
tdata850_25=pd.read_table(r'E:\实习\天气学诊断\micaps\micaps\t\850\13052520.000',header=None,skiprows=4,sep='\s+')
tdata850_25=np.array(tdata850_25)
tdata850_25=micaps(tdata850_25)
tdata850_26=pd.read_table(r'E:\实习\天气学诊断\micaps\micaps\t\850\13052620.000',header=None,skiprows=4,sep='\s+')
tdata850_26=np.array(tdata850_26)
tdata850_26=micaps(tdata850_26)
tddata500_25=pd.read_table(r'E:\实习\天气学诊断\micaps\micaps\t-td\500\13052520.000',header=None,skiprows=4,sep='\s+')
tddata500_25=np.array(tddata500_25)
tddata500_25=micaps(tddata500_25)
tddata500_25=tdata500_25-tddata500_25
tddata500_26=pd.read_table(r'E:\实习\天气学诊断\micaps\micaps\t-td\500\13052620.000',header=None,skiprows=4,sep='\s+')
tddata500_26=np.array(tddata500_26)
tddata500_26=micaps(tddata500_26)
tddata500_26=tdata500_26-tddata500_26
tddata850_25=pd.read_table(r'E:\实习\天气学诊断\micaps\micaps\t-td\850\13052520.000',header=None,skiprows=4,sep='\s+')
tddata850_25=np.array(tddata850_25)
tddata850_25=micaps(tddata850_25)
tddata850_25=tdata850_25-tddata850_25
tddata850_26=pd.read_table(r'E:\实习\天气学诊断\micaps\micaps\t-td\850\13052620.000',header=None,skiprows=4,sep='\s+')
tddata850_26=np.array(tddata850_26)
tddata850_26=micaps(tddata850_26)
tddata850_25=tdata850_25-tddata850_25
##位温计算
thetadata500_25=np.zeros((29,53))
thetadata500_26=np.zeros((29,53))
thetadata850_25=np.zeros((29,53))
thetadata850_26=np.zeros((29,53))
for i in range(0,29):
    for j in range(0,53): 
        thetadata500_25[i,j]=ww(tdata500_25[i,j],500)
        thetadata500_26[i,j]=ww(tdata500_26[i,j],500)
        thetadata850_25[i,j]=ww(tdata500_25[i,j],850)
        thetadata850_26[i,j]=ww(tdata500_26[i,j],850)
##位温画图
draw1(thetadata500_25,30,160,10,80,2.5,'500hPa位温图  2013年5月25日20时  单位：℃')
draw1(thetadata500_26,30,160,10,80,2.5,'500hPa位温图  2013年5月26日20时  单位：℃')
draw1(thetadata850_25,30,160,10,80,2.5,'850hPa位温图  2013年5月25日20时  单位：℃')
draw1(thetadata850_25,30,160,10,80,2.5,'850hPa位温图  2013年5月26日20时  单位：℃')
##相当位温计算
theta_edata500_25=np.zeros((29,53))
theta_edata500_26=np.zeros((29,53))
theta_edata850_25=np.zeros((29,53))
theta_edata850_26=np.zeros((29,53))
for i in range(0,29):
    for j in range(0,53): 
        theta_edata500_25[i,j]=xdww(tdata500_25[i,j],tddata500_25[i,j],500)
        theta_edata500_26[i,j]=xdww(tdata500_26[i,j],tddata500_26[i,j],500)
        theta_edata850_25[i,j]=xdww(tdata500_25[i,j],tddata500_25[i,j],850)
        theta_edata850_26[i,j]=xdww(tdata500_26[i,j],tddata500_26[i,j],850)        
##相当位温画图
draw2(theta_edata500_25,30,160,10,80,2.5,'500hPa相当位温图  2013年5月25日20时  单位：℃')
draw2(theta_edata500_26,30,160,10,80,2.5,'500hPa相当位温图  2013年5月26日20时  单位：℃')
draw2(theta_edata850_25,30,160,10,80,2.5,'850hPa相当位温图  2013年5月25日20时  单位：℃')
draw2(theta_edata850_25,30,160,10,80,2.5,'850hPa相当位温图  2013年5月26日20时  单位：℃')

三、计算结果

例图：