準備¶

必要なパッケージをインストールする。

$ conda install -c conda-forge cartopy metpy

標準大気¶

高層観測データを扱う前に，国際（米国）標準大気の気温分布をプロットして，matploblibの簡単な使い方について学ぶ。

In [1]:

%matplotlib inline
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(5, 5))
T = [15.0, -56.5, -56.5, -44.5, -2.5, -2.5, -58.5, -86.2]
h = [0., 11, 20., 32., 47, 51., 71, 84.852]
ax.plot(T, h)

Out[1]:

[<matplotlib.lines.Line2D at 0x1049628d0>]

高層気象観測データ¶

ここでは，Metpy Air Sounding Tutorial に基づいて高層気象観測の解析と描画を行う。

ラジオゾンデによる高層気象観測

データの取得¶

気象庁からも，日本の高層気象観測データは取得できるが，扱いにくいので，この実習ではWyoming大学にアーカイブされている観測データを取得する。手順は次の通り。

urllibでサーバからデータを取得する。
Beautifulsoup4でHTMLを解析する。
ioで文字列をテキストストリームに変換する。
pandasでデータを読み取る。

In [1]:

from urllib.request import urlopen
from bs4 import BeautifulSoup
import io
import pandas as pd
url = "http://weather.uwyo.edu/cgi-bin/sounding?region=seasia&TYPE=TEXT%3ALIST&YEAR=2018&MONTH=09&FROM=0400&TO=0400&STNM=47778"
response = urlopen(url)
html = response.read().decode("utf8")
soup = BeautifulSoup(html, "html.parser")
s = soup.pre.string
f = io.StringIO(s)
col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
df = pd.read_fwf(f, skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)

In [2]:

df

Out[2]:

	pressure	height	temperature	dewpoint	direction	speed
0	1000.0	8	NaN	NaN	NaN	NaN
1	990.0	75	27.4	24.7	120.0	24.0
2	972.0	236	26.1	23.8	125.0	47.0
3	925.0	671	22.6	21.4	145.0	54.0
4	850.0	1405	17.8	17.3	155.0	66.0
5	775.0	2192	13.2	12.4	163.0	72.0
6	754.0	2422	12.4	10.9	165.0	74.0
7	700.0	3046	10.4	6.9	165.0	74.0
8	696.0	3094	10.4	5.5	165.0	74.0
9	654.0	3608	6.2	2.1	165.0	78.0
10	628.0	3939	4.3	2.8	165.0	80.0
11	619.0	4057	3.6	3.1	165.0	80.0
12	572.0	4696	1.4	1.3	163.0	80.0
13	514.0	5550	-2.6	-2.8	160.0	80.0
14	500.0	5770	-3.7	-3.8	160.0	82.0
15	474.0	6190	-5.9	-6.5	160.0	79.0
16	473.0	6207	-6.1	-12.1	160.0	79.0
17	470.0	6257	-6.9	-14.9	160.0	79.0
18	429.0	6962	-11.9	-14.9	160.0	74.0
19	425.0	7034	-11.0	-13.8	160.0	74.0
20	420.0	7124	-9.9	-12.5	164.0	74.0
21	400.0	7500	-10.9	-13.4	180.0	76.0
22	387.0	7751	-12.2	-14.8	180.0	74.0
23	369.0	8112	-14.1	-16.8	185.0	68.0
24	325.0	9077	-19.1	-22.2	194.0	63.0
25	300.0	9670	-22.7	-26.4	200.0	60.0
26	277.0	10250	-26.7	-32.7	200.0	53.0
27	250.0	10980	-32.7	-37.4	200.0	43.0
28	236.0	11384	-36.3	-40.0	202.0	40.0
29	222.0	11806	-40.1	-43.1	203.0	36.0
30	210.0	12181	-43.4	NaN	205.0	33.0
31	200.0	12510	-46.3	NaN	210.0	31.0
32	180.0	13197	-52.7	NaN	207.0	35.0
33	165.0	13742	-57.9	NaN	205.0	39.0
34	150.0	14340	-63.5	NaN	205.0	43.0
35	133.0	15063	-69.0	NaN	210.0	43.0
36	125.0	15436	-71.9	NaN	220.0	31.0
37	119.0	15725	-72.4	NaN	225.0	19.0
38	113.0	16028	-72.8	NaN	255.0	14.0
39	107.0	16348	-73.4	NaN	265.0	10.0
40	101.0	16686	-73.9	NaN	NaN	NaN

データの解析¶

In [3]:

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.cbook import get_test_data
from metpy.plots import Hodograph, SkewT
from metpy.units import units

In [4]:

df['u_wind'], df['v_wind'] = mpcalc.wind_components(df['speed'], np.deg2rad(df['direction']))
df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed',
                       'u_wind', 'v_wind'), how='all').reset_index(drop=True)

In [5]:

p = df['pressure'].values * units.hPa
T = df['temperature'].values * units.degC
Td = df['dewpoint'].values * units.degC
wind_speed = df['speed'].values * units.knots
wind_dir = df['direction'].values * units.degrees
u, v = mpcalc.wind_components(wind_speed, wind_dir)

In [6]:

T[0]

Out[6]:

27.4 degC

In [7]:

p[0]

Out[7]:

990.0 hectopascal

持ち上げ凝結高度の計算

In [7]:

lcl_pressure, lcl_temperature = mpcalc.lcl(p[0], T[0], Td[0])
print(lcl_pressure, lcl_temperature)

951.8403007544475 hectopascal 24.04446396516281 degC

上昇気塊の鉛直分布の計算

In [8]:

parcel_prof = mpcalc.parcel_profile(p, T[0], Td[0]).to('degC')

matplotlibでプロット¶

In [9]:

fig, ax = plt.subplots(figsize=(5, 5))
ax.plot(T, p)
ax.plot(Td, p)

Out[9]:

[<matplotlib.lines.Line2D at 0x10bd46a90>]

y軸を対数にし，反転させて，目盛をつける。

In [10]:

from matplotlib.ticker import ScalarFormatter
fig, ax = plt.subplots(figsize=(5, 5))
ax.semilogy(T, p)
ax.semilogy(Td, p)
ax.set_yticks([1000, 850, 700, 500, 300, 200, 100])
ax.get_yaxis().set_major_formatter(ScalarFormatter())
ax.invert_yaxis()

簡単なSkew-T図¶

In [11]:

fig = plt.figure(figsize=(9, 9))
skew = SkewT(fig)
skew.plot(p, T, 'r', linewidth=2)
skew.plot(p, Td, 'g', linewidth=2)
skew.plot_barbs(p, u, v)

Out[11]:

<matplotlib.quiver.Barbs at 0x10beb5080>

Skew-T図を完成させる¶

In [12]:

fig = plt.figure(figsize=(9, 9))
skew = SkewT(fig)
skew.plot(p, T, 'r', linewidth=2)
skew.plot(p, Td, 'b', linewidth=2)
skew.plot_barbs(p, u, v)
skew.ax.set_ylim(1000, 100)
skew.ax.set_xlim(-40, 60)
skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')
skew.plot(p, parcel_prof, 'k', linewidth=2)
skew.shade_cin(p, T, parcel_prof)
skew.shade_cape(p, T, parcel_prof)
skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()

Out[12]:

<matplotlib.collections.LineCollection at 0x10c2160f0>

ホドグラフの追加¶

In [13]:

fig = plt.figure(figsize=(9, 9))
skew = SkewT(fig)
skew.plot(p, T, 'r', linewidth=2)
skew.plot(p, Td, 'b', linewidth=2)
skew.plot_barbs(p, u, v)
skew.ax.set_ylim(1000, 100)
skew.ax.set_xlim(-40, 60)
skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')
skew.plot(p, parcel_prof, 'k', linewidth=2)
skew.shade_cin(p, T, parcel_prof)
skew.shade_cape(p, T, parcel_prof)
skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()

ax_hod = inset_axes(skew.ax, '40%', '40%', loc=1)
h = Hodograph(ax_hod, component_range=80)
h.add_grid(increment=20)
h.plot_colormapped(u, v, wind_speed)

Out[13]:

<matplotlib.collections.LineCollection at 0x10c4c5da0>

関数にまとめる¶

In [14]:

from urllib.request import urlopen
from bs4 import BeautifulSoup
import io
import pandas as pd
def get_sounding(yyyy, mm, dd, hh, stnm):
    url = "http://weather.uwyo.edu/cgi-bin/sounding?TYPE=TEXT%3ALIST&YEAR={0}&MONTH={1:0=2}&FROM={2:0=2}{3:0=2}&TO={2:0=2}{3:0=2}&STNM={4}"
    response = urlopen(url.format(yyyy,mm,dd,hh,stnm))
    html = response.read().decode("utf8")
    soup = BeautifulSoup(html, "html.parser")
    s = soup.pre.string
    f = io.StringIO(s)
    col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
    df = pd.read_fwf(f, skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
    return df

In [15]:

df = get_sounding(2018, 8, 3, 0, 47778)

In [16]:

df

Out[16]:

	pressure	height	temperature	dewpoint	direction	speed
0	1001.0	75	28.2	25.8	260	12
1	1000.0	76	27.6	25.3	265	12
2	997.0	103	27.0	25.0	269	12
3	949.0	540	27.6	19.6	332	12
4	932.0	700	27.0	19.0	355	12
5	925.0	767	26.8	18.8	5	12
6	886.0	1147	24.1	16.1	5	4
7	881.0	1197	23.8	15.8	358	5
8	854.0	1469	21.7	15.4	320	10
9	850.0	1510	21.4	15.4	325	10
10	782.0	2227	17.4	7.5	330	12
11	765.0	2415	16.4	5.4	330	10
12	721.0	2915	13.1	3.4	330	6
13	700.0	3164	11.4	2.4	305	8
14	695.0	3224	11.1	1.7	305	8
15	675.0	3467	9.8	-1.2	328	7
16	641.0	3893	8.0	-6.9	10	6
17	631.0	4023	7.4	-8.6	4	7
18	628.0	4062	7.2	-2.8	2	7
19	617.0	4207	6.2	-3.0	355	8
20	588.0	4600	3.4	-3.6	9	9
21	586.0	4627	3.6	-13.4	10	9
22	576.0	4767	3.6	-20.4	15	9
23	569.0	4865	2.4	-8.6	19	9
24	564.0	4937	2.4	-13.6	22	9
25	554.0	5081	1.6	-3.1	27	10
26	539.0	5301	0.4	-6.8	35	10
27	517.0	5634	-1.5	-12.5	18	10
28	514.0	5680	-1.5	-8.5	16	10
29	510.0	5743	-1.3	-11.3	13	10
...	...	...	...	...	...	...
38	431.0	7064	-8.9	-16.3	45	10
39	415.0	7358	-10.9	-19.9	48	11
40	413.0	7395	-10.9	-25.9	48	11
41	400.0	7640	-12.5	-19.5	50	12
42	399.0	7659	-12.7	-19.7	50	12
43	368.0	8271	-16.3	-27.3	56	16
44	343.0	8795	-20.3	-27.3	61	20
45	324.0	9213	-23.3	-41.3	65	23
46	319.0	9327	-24.1	-45.1	68	24
47	300.0	9770	-27.9	-45.9	80	29
48	278.0	10305	-32.5	-48.8	75	31
49	263.0	10694	-35.8	-51.0	75	25
50	250.0	11050	-38.9	-52.9	65	27
51	245.0	11188	-40.1	-53.1	61	27
52	200.0	12540	-50.7	NaN	25	31
53	194.0	12734	-52.0	NaN	30	31
54	170.0	13574	-57.7	NaN	55	35
55	150.0	14370	-63.1	NaN	45	27
56	149.0	14411	-63.3	NaN	45	29
57	145.0	14577	-64.3	NaN	47	33
58	130.0	15228	-67.3	NaN	55	49
59	121.0	15656	-69.2	NaN	55	47
60	114.0	16011	-70.8	NaN	80	47
61	106.0	16444	-72.8	NaN	90	31
62	101.0	16732	-74.1	NaN	70	25
63	100.0	16790	-73.7	NaN	65	27
64	93.0	17216	-72.3	NaN	60	41
65	89.8	17421	-71.7	NaN	73	37
66	86.0	17676	-72.3	NaN	90	31
67	83.2	17870	-72.7	NaN	71	25

68 rows × 6 columns

In [17]:

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.cbook import get_test_data
from metpy.plots import Hodograph, SkewT
from metpy.units import units

def plot_sounding(df):
    df['u_wind'], df['v_wind'] = mpcalc.wind_components(df['speed'], np.deg2rad(df['direction']))
    df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed',
                       'u_wind', 'v_wind'), how='all').reset_index(drop=True)
    p = df['pressure'].values * units.hPa
    T = df['temperature'].values * units.degC
    Td = df['dewpoint'].values * units.degC
    wind_speed = df['speed'].values * units.knots
    wind_dir = df['direction'].values * units.degrees
    u, v = mpcalc.wind_components(wind_speed, wind_dir)
    lcl_pressure, lcl_temperature = mpcalc.lcl(p[0], T[0], Td[0])
    
    parcel_prof = mpcalc.parcel_profile(p, T[0], Td[0]).to('degC')
    fig = plt.figure(figsize=(9, 9))
    skew = SkewT(fig)
    skew.plot(p, T, 'r', linewidth=2)
    skew.plot(p, Td, 'b', linewidth=2)
    skew.plot_barbs(p, u, v)
    skew.ax.set_ylim(1000, 100)
    skew.ax.set_xlim(-40, 60)
    skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')
    skew.plot(p, parcel_prof, 'k', linewidth=2)
    skew.shade_cin(p, T, parcel_prof)
    skew.shade_cape(p, T, parcel_prof)
    skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)
    skew.plot_dry_adiabats()
    skew.plot_moist_adiabats()
    skew.plot_mixing_lines()

    ax_hod = inset_axes(skew.ax, '40%', '40%', loc=1)
    h = Hodograph(ax_hod, component_range=80)
    h.add_grid(increment=20)
    h.plot_colormapped(u, v, wind_speed)

In [18]:

plot_sounding(df)

In [ ]: