来自粉丝“青羽”的投稿
文章很长,可以直接看末尾的用法总结。
参考:
Metpy官方文档
Metpy Skew-T Tutorials
Metpy Skew-T Complex Layout
Siphon Wyoming Upper Air Data Request
用于天气绘图的Metpy包更新(0.8版本)了,他们要逐渐抛弃Python2.X,转到Python>=3.6的版本上。所以,之前(越2018年6月以前,0.7版本)的一些脚本就无法使用了。在大气科学专业,我们主要使用 Metpy 绘制以怀俄明大学高空探测数据为基础的斜-T图(Skew-T)。受更新影响,原本的 upperair_sounding.py脚本也有改动。
老版本的第17~41行为:
from datetime import datetime
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import numpy as np
import metpy.calc as mpcalc
from metpy.io import get_upper_air_data
from metpy.plots import Hodograph, SkewT
#########################################################################
# Getting Data
# ------------
#
# We will download data from the
# `University of Wyoming sounding data page <http://weather.uwyo.edu/upperair/sounding.html>`_
# , which has an extensive archive of data available, as well as current data.
#
# In this case, we will download the sounding data from the Veterans Day
# tornado outbreak in 2002 by passing a ``datetime`` object and station name to the
# ``get_upper_air_data`` function.
dataset = get_upper_air_data(datetime(2018, 5, 7, 0), 'ZSQD')
##########################################################################
而新版本的第18~47行为:
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import numpy as np
import pandas as pd
import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.plots import Hodograph, SkewT
from metpy.units import units
#########################################################################
# Getting Data
# ------------
#
# Upper air data can be obtained using the siphon package, but for this tutorial we will use
# some of MetPy's sample data. This event is the Veterans Day tornado outbreak in 2002.
col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
df = pd.read_fwf(get_test_data('nov11_sounding.txt', as_file_obj=False),
skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
df['u_wind'], df['v_wind'] = mpcalc.wind_components(df['speed'],
np.deg2rad(df['direction']))
# Drop any rows with all NaN values for T, Td, winds
df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed',
'u_wind', 'v_wind'), how='all').reset_index(drop=True)
##########################################################################
对比可以发现,用于绘图的模块基本没变,而是用于读取数据的模块变了,无论是库还是代码调用方法。
from metpy.io import get_upper_air_data
---
from metpy.cbook import get_test_data
dataset = get_upper_air_data(datetime(2018, 5, 7, 0), 'ZSQD')
---
df = pd.read_fwf(get_test_data('nov11_sounding.txt', as_file_obj=False),
skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
而从官方说明文档中得知,获取数据这个功能被独立出来,用Siphon这个库实现。让我们来看看Siphon时何方神圣?搜索Siphon文档得知,这货就是个专门爬气象数据的包。他提供
之前我喜欢用怀俄明大学的数据,所以这次就来看看要怎么用怀俄明大学的数据画Skew-T图。
分析新版本的数据读取的语句
df = pd.read_fwf(get_test_data('nov11_sounding.txt', as_file_obj=False),
skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
得知,nov11_sounding.txt 是官方内置的一个示例数据文件,这行代码主要是用pandos的功能读取一个列表数据。
如果能用Siphon下载同格式的数据的话,那么只要改改读取的文件名就可以了。我们看一下这个神奇的txt文件在哪里 运行一下脚本
python3 upperair_soundings.py
---
/home/bugatti/anaconda3/lib/python3.6/site-packages/pooch/core.py:260: UserWarning: Downloading data file 'nov11_sounding.txt' from remote data store 'https://github.com/Unidata/MetPy/raw/v0.10.0/staticdata/' to '/home/bugatti/.cache/metpy/v0.10.0'.
action, fname, self.base_url, str(self.path)
923.0963220202037 hectopascal 15.595422727825456 degC
/home/bugatti/anaconda3/lib/python3.6/site-packages/pint/quantity.py:1377: UnitStrippedWarning: The unit of the quantity is stripped.
warnings.warn("The unit of the quantity is stripped.", UnitStrippedWarning)
能够成功出图,然后看到命令行窗口的第二行显示这个文件在'/home/bugatti/.cache/metpy/v0.10.0' 这个目录里,bugatti是我的用户名,换在你的电脑中应该也是类似的路径。打开这个文件格式如下:
-----------------------------------------------------------------------------
PRES HGHT TEMP DWPT RELH MIXR DRCT SKNT THTA THTE THTV
hPa m C C % g/kg deg knot K K K
-----------------------------------------------------------------------------
1000.0 -12
978.0 180 20.4 16.5 78 12.22 180 16 295.4 330.7 297.6
964.1 305 22.2 17.1 73 12.92 185 29 298.5 336.3 300.8
954.0 397 23.6 17.6 69 13.45 188 35 300.8 340.5 303.2
931.0 610 22.5 16.5 69 12.84 195 49 301.7 339.8 304.1
925.0 667 22.2 16.2 69 12.68 200 48 302.0 339.7 304.3
898.9 914 20.2 14.5 70 11.68 205 49 302.4 337.2 304.5
867.6 1219 17.7 12.4 71 10.54 215 52 302.9 334.4 304.8
850.0 1396 16.2 11.2 72 9.92 220 55 303.1 332.9 304.9
这个和在怀俄明无线电探空数据网站 的数据格式一模一样(毕竟就是从这爬取的数据嘛)
from datetime import datetime
from metpy.units import units
from siphon.simplewebservice.wyoming import WyomingUpperAir
date = datetime(2017, 9, 10, 6)
station = 'MFL'
df = WyomingUpperAir.request_data(date, station)
就是这么简单粗暴!
在Spyder中运行上述代码,获取数据,在变量窗口查看数据值
与怀俄明大学网站的数据列表作对比
发现 (2017, 9, 10, 6)表示的就是2017年9月10日6时。6所在位置的参数只有(0, 6, 12)三个数字可选。这里使用的是世界时,0和12也就是中国地区早八晚八时放气球,美国迈阿密(站号MFL,时区西五区)大概是晚七早七,6属于补测的时间,一般是由于天气过程复杂,或者早上的数据没测好才会有补测的6时。
回到数据读取问题上来,示例文件中的语句是把一个写好的txt文件读取为类似nc文件那样带变量描述的值的一个东西。
df = pd.read_fwf(get_test_data('nov11_sounding.txt', as_file_obj=False),
skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
而siphon爬去网络数据就直接是这样的格式,所以如果要改为自己想要的某天日期也十分简单。就是把示例文件中的语句注释掉,然后添加如下代码:
#df = pd.read_fwf(get_test_data('nov11_sounding.txt', as_file_obj=False),
# skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
from datetime import datetime
from siphon.simplewebservice.wyoming import WyomingUpperAir
date = datetime(2019, 3, 25, 0)
station = 'MFL'
df = WyomingUpperAir.request_data(date, station)
运行得到本文写作当天早八点的探空图。嗯,看样子多云天还得持续几天。
print('欢迎使用 Metpy 和 Siphon 获取怀俄明大学无线电探空数据,使用方法如下:')
print('1.输入’年月日时‘和站号,如‘2019032500ZSQD’,中间无间隔。')
print('2.输入图片保存路径,回车默认脚本所在文件夹。')
print('提示:探空时段采用世界时,一般为00或12,06为补测数据,请自行换算时区。')
print(' 中国区常用站号:北京-ZBAA,青岛-ZSQD,昆明-ZPPP,南京-ZSNJ')
dateNum=input('请输入日期和站号')
figpath=input('请输入保存路径')
print('正在获取探空图',dateNum[:4], dateNum[4:6], dateNum[6:8], dateNum[8:10],dateNum[10:])
print('保存路径为',figpath)
点击下载代码
(PS:Siphon的意思是‘虹吸管’,很适合它,吸数据的虹吸管)
转载请注:Metpy新版功能下载TLnP图设置 - Bugatii100Peagle
更多详见:https://zhuanlan.zhihu.com/p/102032513