北大钱维宏钱教授的《瞬变涡扰动法在极端天气事件预报中的应用》，使用python对2017年台风“海棠”登陆后的天气形势、路径、水汽进行了一次详尽地分析，
知识链接：
瞬变扰动是指瞬变纬向波动，瞬变且纬向不对称部分，比如天气系统等。中高纬度、大气环流天气尺度变化十分明显，这种每日天气图上与时间有关的系统或流型变化主要 由大气中的瞬变波引起(丁一汇，2005)。台风降水是夏季常见的天气现象，对其进行 瞬变扰动分析符合瞬变扰动特点。
钱维宏（2012）认为，实际天气图上的系统由多尺度波和多尺度涡组成。Rossby对天气预报的重大贡献是提出了大气长波。全球数值天气预报谱模式把观测大气变量分解成了若干尺度的数学波，而没有天气尺度的涡。大气变量物理分解可得到逐日气候波、行星尺度的瞬变波和天气尺度的瞬变涡。实际的天气变化是由逐日气候波、行星尺度瞬变波和天气尺度瞬变涡的叠加或相互作用形成的。与极端天气事件直接有关的是那些瞬变涡，因此极端天气事件的预报可用瞬变涡扰动法。
数据获取：
欧洲中心(ECMWF, European Center for Medium-range Weather Forecasts)提供的全球大气再分析资料，时间间隔6小时，即每天00 UTC、06 UTC、12 UTC、18 UTC资料，水平格点为1°×1°，包括：降水、位势高度、气温、经纬向风速、相对湿度、位势涡度等资料
from mpl_toolkits.basemap import Basemap, cm
# requires netcdf4-python (netcdf4-python.googlecode.com)
from netCDF4 import Dataset as NetCDFFile
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm as CM
from matplotlib.colors import BoundaryNorm
from matplotlib.ticker import MaxNLocator
from mpl_toolkits.basemap import Basemap
import math
import matplotlib.tri as tri
import pyhdf.SD as hdf
import scipy.io as sio
from pyhdf import SD
import glob
import os
import shutil
import zipfile
from os.path import join, getsize
import sys, csv , operator
nc1 = NetCDFFile('F:\grads\mo.nc')
lon = nc1.variables['longitude'][:]
lat = nc1.variables['latitude'][:]
time = nc1.variables['time'][:]
level = nc1.variables['level'][:]
Geoooo = nc1.variables['z'][:,:,:,:].squeeze()
Geooo = Geoooo/10+120
Geoo = Geooo[66,21,:,:].squeeze()
Temmmm = nc1.variables['t'][:,:,:,:].squeeze()
Temmm = Temmmm-273
Temm = Temmm[66,30,:,:].squeeze()
Geo_original_field = Geoo/1
Tem_original_field = Temm/1
#我是华丽的分割线我是华丽的分割线我是华丽的分割线我是华丽的分割线我是华丽的分割线我是华丽的分割线我是华丽的分割线
Geo_time_ave = np.nanmean(Geooo[0:66,21,:,:],axis=0) #这里是9月1号到9月10号大约9天的时间平均数据
Geo_time_transient = Geoo - Geo_time_ave#这里是9月10号早上6点的数据减去时间平均得到的扰动数据
#这就是时间扰动啦！
#下面开始做时间扰动的纬向平均
Geo_time_lat_ave = np.nanmean(Geo_time_transient [100:150,:],axis=0)#纬向度我先选个全球的
#最后得出瞬变扰动，用时间扰动-时间扰动的纬向平均
Geo_transient = Geo_time_transient - Geo_time_lat_ave
#大功告成！ 这就是位势高度场的瞬变！
#同理温度场如下
Tem_time_ave = np.nanmean(Temmm[0:66,30,:,:],axis=0)
Tem_time_transient = Temm - Tem_time_ave
Tem_time_lat_ave = np.nanmean(Tem_time_transient [100:150,:],axis=0)#纬向度我先选个全球的
Tem_transient = Tem_time_transient - Tem_time_lat_ave
#大功告成！ 这就是温度场的瞬变！
###数据的顺序【时间（0-247 共248个时间点），层数（11层），纬度（101刀），经度（360刀）】
###中国的区域先看下面 不会编了再说
"""
pv是Potential vorticity 位势涡度
z是Geopotential 位势高度
t是Temperature 温度
w是Vertical velocity 垂直速度
r是Relative humidity 相对湿度
u是U 纬向风
v是V 经向风
"""
######ices = np.nanmean(ice[0:40,:,:],axis=0) 这是友军 用来对照写程序的
#你踏马能不能好好编程序别看动漫了？？？？？！！！！！
m = Basemap(projection='stere',lon_0=120,lat_0=40,lat_ts=0,\
llcrnrlat=0,urcrnrlat=50,\
llcrnrlon=80,urcrnrlon=160,\
rsphere=6371200.,resolution='l',area_thresh=10000)
"""
m = Basemap(projection='hammer',lon_0=180)
#m = Basemap(projection='ortho',lat_0=90,lon_0=0)
# llcrnrx = -4100000,llcrnry = -4100000,urcrnrx = 4100000,urcrnry = 4100000,
# resolution='l')
#北半球的哪条纬线没有看过？你敢不敢作经度的剖面平均？
"""
plt.figure(figsize=(10,12))#竖着8到10 躺着10到8
m.drawcoastlines() #海岸线
m.drawstates() #省市线
m.drawcountries() #郭姐线
parallels = np.arange(-100.,100,10.)
m.drawparallels(parallels,labels=[0,0,0,0],fontsize=12)
meridians = np.arange(-180.,180.,20.)
m.drawmeridians(meridians,labels=[1,1,0,1],fontsize=12)
lats, lons = np.meshgrid(lat, lon) # lat和lon排序略显混乱 可以试试改变先后顺序
x, y = m(lons, lats) # compute map proj coordinates.
csq = m.contourf(x,y,Tem_original_field.transpose(),np.arange(0,30,0.01),cmap = 'jet') #blue可以用于鄙视场
cbar = m.colorbar(csq,location='right',size="6%",pad="10%")
csr= m.contour(x, y,Geo_original_field.transpose(), np.arange(5000,6000,20), colors = 'black', linewidth = 0.08)
plt.clabel(csr, inline = True, fontsize = 10)# 显示各等高线的数据标签
plt.title('9_17_UTC12_00_500hPaGEO_850hPaTEM_original_field')
#plt.savefig("9_17_UTC12_00_500hPaGEO_850hPaTEM_original_field.png",dpi=1000)
#我是华丽的分割线我是华丽的分割线我是华丽的分割线我是华丽的分割线我是华丽的分割线我是华丽的分割线我是华丽的分割线
plt.figure(figsize=(10,12))
m.drawcoastlines()
m.drawstates()
m.drawcountries()
parallels = np.arange(-100.,100,10.)
m.drawparallels(parallels,labels=[0,0,0,0],fontsize=12)
meridians = np.arange(-180.,180.,20.)
m.drawmeridians(meridians,labels=[1,1,0,1],fontsize=12)
lats, lons = np.meshgrid(lat, lon) # get lat/lons of ny by nx evenly space grid.
x, y = m(lons, lats) # compute map proj coordinates.
csq = m.contourf(x,y,Tem_transient.transpose(),np.arange(-9,9,0.25),cmap = 'jet')
cbar = m.colorbar(csq,location='right',size="6%",pad="10%")
csr= m.contour(x, y,Geo_transient.transpose(), np.arange(-160,80,15), colors = 'black', linewidth = 0.1)
plt.clabel(csr, inline = True, fontsize = 15)# 显示各等高线的数据标签
plt.title('9_17_UTC12_00_500hPaGEO_850hPaTEM_transient_field')
#plt.savefig("9_17_UTC12_00_500hPaGEO_850hPaTEM_transient_field.png",dpi=1000)
网友评论