python气象信息系统工程 第九章

第9章 基本绘图类型与气象绘图

9.1 折线图

9.1.1 基本折线图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
#读取数据
ao = pd.read_csv("/home/mw/input/pythonbook9857/AO.txt",sep='\s+',header=None, names=['year','month','AO'])  
ao_jan = ao[ao.month==1]
#创建Figure
fig = plt.figure(figsize=(10, 8))
#创建Axes
ax1 = fig.add_subplot(1,1,1)
#绘制折线图
ax1.plot(np.arange(1950,2020,1),ao_jan.AO, 'ko-')
#添加图题
ax1.set_title('1950-2019 January AO Index')
#添加y=0水平参考线
ax1.axhline(0,ls=':',c='r')
#添加x=1990垂直参考线
ax1.axvline(1990,ls='--',c='g')
plt.show()

9.1.2 多折线图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
ao = pd.read_csv("/home/mw/input/pythonbook9857/AO.txt",sep='\s+',header=None, names=['year','month','AO'])  
ao_jan = ao[ao.month==1]
ao_feb = ao[ao.month==2]
fig = plt.figure(figsize=(10, 8))
ax1 = fig.add_subplot(2,1,1)
#绘制第一个序列
ax1.plot(np.arange(1950,2020,1),ao_jan.AO,'ko-',label='Jan')
#绘制第二个序列
ax1.plot(np.arange(1950,2020,1),ao_feb.AO,'rs-',label='Feb')
ax1.set_title('1950-2019 AO Index')
#添加图例
ax1.legend()
plt.show()

9.1.3 多y轴折线图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import matplotlib.pyplot as plt
import numpy as np
#读取数据
ao = pd.read_csv("/home/mw/input/pythonbook9857/AO.txt",sep='\s+',header=None, names=['year','month','AO'])  
ao_jan = ao[ao.month==1]
ao_feb = ao[ao.month==2]
#创建Figure
fig = plt.figure()
#绘制单y轴折线图
ax1 = fig.add_subplot(1,1,1)
ax1.plot(np.arange(1950,2020,1),ao_jan.AO,'ko-',label='Jan')
ax1.set_ylabel('January',c='r')
ax1.set_title('1950-2019 AO Index')
#创建第二个y轴
ax2 = ax1.twinx()  
ax2.plot(np.arange(1950,2020,1),ao_feb.AO,'rs-',label='Feb')
ax2.set_ylim(-4,4)
ax2.set_ylabel('February',c='k')

Text(0, 0.5, 'February')

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
from mpl_toolkits.axisartist.parasite_axes import HostAxes, ParasiteAxes
import matplotlib.pyplot as plt
import numpy as np
#读取数据
ao = pd.read_csv("/home/mw/input/pythonbook9857/AO.txt",sep='\s+',header=None, names=['year','month','AO'])  
ao_dec = ao[ao.month==12]
ao_jan = ao[ao.month==1]
ao_feb = ao[ao.month==2]
ao_djf = (ao_jan.AO.values+ao_feb.AO.values+ao_dec.AO.values)/3
#创建Figure
fig = plt.figure() 
#创建主轴
ax = HostAxes(fig, [0, 0, 0.9, 0.9])  #[left, bottom, weight, height]
#创建共享x轴的其他y轴
ax1 = ParasiteAxes(ax, sharex=ax)
ax2 = ParasiteAxes(ax, sharex=ax)
ax3 = ParasiteAxes(ax, sharex=ax)
ax.parasites.append(ax1)
ax.parasites.append(ax2)
ax.parasites.append(ax3)
#将主轴的右轴隐藏，同时开始设置第二个轴的可见性
ax.axis['right'].set_visible(False)
ax1.axis['right'].set_visible(True)
ax1.axis['right'].major_ticklabels.set_visible(True)
ax1.axis['right'].label.set_visible(True)
#设置各轴的标签
ax.set_ylabel('DJF')
ax.set_xlabel('year')
ax1.set_ylabel('Dec')
ax2.set_ylabel('Jan')
ax3.set_ylabel('Feb')
#设置第三个和第四个y轴的位置
axisline2 = ax2.get_grid_helper().new_fixed_axis
axisline3 = ax3.get_grid_helper().new_fixed_axis
ax2.axis['right2'] = axisline2(loc='right', axes=ax2, offset=(40,0))
ax3.axis['right3'] = axisline3(loc='right', axes=ax3, offset=(80,0))
#将设置好的主轴的Axes放在Figure上
fig.add_axes(ax)
#绘制折线
ax.plot(np.arange(1950,2020,1), ao_djf, label="DJF", color='black')
ax1.plot(np.arange(1950,2020,1), ao_dec.AO, label="Dec", color='red')
ax2.plot(np.arange(1950,2020,1), ao_jan.AO, label="Jan", color='green')
ax3.plot(np.arange(1950,2020,1), ao_feb.AO, label="Feb", color='orange')
ax1.set_ylim(-4,4)
ax2.set_ylim(-5,5)
ax3.set_ylim(-6,6)
ax.legend()
#设置各个轴及其刻度的颜色
ax1.axis['right'].major_ticks.set_color('red')
ax2.axis['right2'].major_ticks.set_color('green')
ax3.axis['right3'].major_ticks.set_color('blue')
ax1.axis['right'].major_ticklabels.set_color('red')
ax2.axis['right2'].major_ticklabels.set_color('green')
ax3.axis['right3'].major_ticklabels.set_color('blue')
ax1.axis['right'].line.set_color('red')
ax2.axis['right2'].line.set_color('green')
ax3.axis['right3'].line.set_color('blue')
plt.show()

9.1.4 非等比坐标轴图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import pandas as pd
import matplotlib.ticker as ticker
#读取数据
f = xr.open_dataset('/home/mw/input/pythonbook9857/data.nc')
t = f['air'].loc['2005-07-01',:,45,120]
lev = t.level
#创建Figure
fig = plt.figure(figsize=(15,5))
#设置为log
ax1 = fig.add_subplot(1,3,1)
ax1.plot(t,lev)
ax1.set_ylim(1000,100)
ax1.set_yscale('log')
ax1.set_title('log')
ax1.set_yticks([1000,850,700,500,300,200,100])
ax1.set_yticklabels([1000,850,700,500,300,200,100])
#隐藏次坐标刻度
ax1.yaxis.set_minor_formatter(ticker.NullFormatter())
#设置为symlog
ax2 = fig.add_subplot(1,3,2)
ax2.plot(t,lev)
ax2.set_ylim(1000,100)
ax2.set_yscale('symlog')
ax2.set_title('symlog')
ax2.set_yticks([1000,850,700,500,300,200,100])
ax2.set_yticklabels([1000,850,700,500,300,200,100])
#设置为linear
ax3 = fig.add_subplot(1,3,3)
ax3.plot(t,lev)
ax3.invert_yaxis()
ax3.set_yscale('linear')
ax3.set_title('linear')
ax3.set_yticks([1000,850,700,500,300,200,100])
ax3.set_yticklabels([1000,850,700,500,300,200,100])
plt.show()

Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.14.1

9.2.1 基础散点图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
import matplotlib.pyplot as plt
import numpy as np
rng = np.random.RandomState(0).rand(5)
markers = ['o','+','*','s','<']
labels = ['CNRM-CM3','GFDL-CM2.0','GISS-AOM','MIROC3.2','CCSM3',]
#绘制散点图
fig = plt.figure()
ax1 = fig.add_subplot(1,1,1)
#通过循环，每次绘制一种标记的散点，循环5次
for i in range(5):
    ax1.scatter(rng,rng,marker=markers[i],label=labels[i])
ax1.legend()
plt.show()

9.2 散点图

9.2.2 带有地图投影的散点图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
import cartopy.io.shapereader as shpreader
data = pd.read_csv("/home/mw/input/pythonbook9857/825station.txt",sep='\s+',header=None, names=['stat','lat','lon','high'])
fig = plt.figure(figsize=(10, 8))
ax1 = fig.add_subplot(1,1,1,projection = ccrs.PlateCarree(central_longitude=105) )
#设置图形范围及刻度
ax1.set_extent([70,110,30,60], crs=ccrs.PlateCarree())
ax1.set_xticks(np.arange(70,120,10), crs=ccrs.PlateCarree())
ax1.set_yticks(np.arange(30,70,10), crs=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#绘制站点的三维分布图
s = ax1.scatter(data.lon,data.lat,s = 5,c = data.high,cmap='jet',transform=ccrs.PlateCarree())
#添加色标，fraction参数用于设置色标缩放比例
fig.colorbar(s,ax=ax1,fraction=0.034)
#添加海岸线等特征
# ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
# china = shpreader.Reader('/data/home/zenggang/yxy/shp/bou2_4l.dbf').geometries()
# ax1.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
plt.show()

9.3 柱状图

9.3.1 单变量柱状图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import pandas as pd
import matplotlib.pyplot as plt
#读取数据
ao = pd.read_csv("/home/mw/input/pythonbook9857/AO.txt",sep='\s+',header=None, names=['year','month','AO'])  
ao_jan = ao[ao.month==1]
#创建Figure及Axes
fig = plt.figure(figsize=(10, 8))
ax1 = fig.add_subplot(1,1,1)
ax1.set_title('1950-2019 January AO Index')
#绘制单变量柱状图
ax1.bar(np.arange(1950,2020,1),ao_jan.AO)
#添加0值参考线
ax1.axhline(0,c='k')
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import pandas as pd
import matplotlib.pyplot as plt
#读取数据
ao = pd.read_csv("/home/mw/input/pythonbook9857/AO.txt",sep='\s+',header=None, names=['year','month','AO'])  
ao_jan = ao[ao.month==1]
#创建颜色数组
colors = np.zeros(ao_jan.AO.shape,dtype=np.str)
colors[ao_jan.AO>=0] = 'red'
colors[ao_jan.AO<0] = 'blue'
#创建Figure及Axes
fig = plt.figure(figsize=(10, 8))
ax1 = fig.add_subplot(1,1,1)
ax1.set_title('1950-2019 January AO Index')
#绘制单变量柱状图
ax1.bar(np.arange(1950,2020,1),ao_jan.AO,color=colors)
#添加0值参考线
ax1.axhline(0,c='k')
plt.show()

/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:7: DeprecationWarning: `np.str` is a deprecated alias for the builtin `str`. To silence this warning, use `str` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.str_` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  import sys

9.3.2 多变量柱状图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import pandas as pd
import matplotlib.pyplot as plt
#数据读取
ao = pd.read_csv("/home/mw/input/pythonbook9857/AO.txt",sep='\s+',header=None, names=['year','month','AO'])  
ao_dec = ao[ao.month==12]
ao_jan = ao[ao.month==1]
ao_feb = ao[ao.month==2]
#创建Figure及Axes
fig = plt.figure(figsize=(10, 8))
ax1 = fig.add_subplot(1,1,1)
ax1.set_title('2000—2019 DJF AO Index')
#绘制柱状图，ao_dec.AO[50:]表明从2000年开始绘制
ax1.bar(np.arange(2000,2020,1)-0.25, ao_dec.AO[50:], width=0.25, color='r', label='Dec')
ax1.bar(np.arange(2000,2020,1),ao_jan.AO[50:],width=0.25,color='b',label='Jan')
ax1.bar(np.arange(2000,2020,1)+0.25,ao_feb.AO[50:],width=0.25,color='g',label='Feb')
#添加0值参考线
ax1.axhline(0,c='k')
ax1.legend()
#由于多变量柱状图对应的x轴的刻度不一定为整数，因此我们指定x轴刻度以避免小数出现
ax1.set_xticks([2000,2005,2010,2015,2020])
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
import numpy as np
import matplotlib.pyplot as plt
#生成数据
labels = ['A', 'B', 'C', 'D']
v1 = np.array([20, 35, 30, 35])
v2 = np.array([25, 32, 34, 20])
v3 = np.array([11, 14, 20, 18])
v1_std = np.array([2, 3, 4, 1])
v2_std = np.array([3, 5, 2, 3])
v3_std = np.array([2, 3, 4, 1])
width = 0.35
fig = plt.figure(figsize=(10, 8))
ax1 = fig.add_subplot(1,1,1)
#绘制垂直堆叠的多变量柱状图
ax1.bar(labels, v1, width, yerr=v1_std, label='v1')
ax1.bar(labels, v2, width, yerr=v2_std, bottom=v1,label='v2')
ax1.bar(labels, v3, width, yerr=v3_std, bottom=v2+v1,label='v3')
ax1.legend()
plt.show()

9.4 箱线图

1
2
3
4
5
6
7
8
9
10
import matplotlib.pyplot as plt
import numpy as np
#生成数据
data=[np.random.normal(0,std,100) for std in range(1,4)]
#绘制箱线图
fig = plt.figure(figsize=(10, 8))
ax1 = fig.add_subplot(1,1,1)
ax1.boxplot(data,vert=True,whis=1.2,patch_artist=True,boxprops={'color':'black','facecolor':'red'})
plt.show()

9.5 等值线图

9.5.1 基本等值线图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import numpy as np
import matplotlib.pyplot as plt
#生成数据
x = np.arange(-10,10,0.01)
y = np.arange(-10,10,0.01)
x_grid,y_grid = np.meshgrid(x,y)
z = x_grid**2+y_grid**2
#创建Figure及Axes
fig = plt.figure(figsize=(10, 8))
ax1 = fig.add_subplot(2,2,1)
ax1.contour(x,y,z)
ax1.set_title('a')
ax2 = fig.add_subplot(2,2,2)
ax2.contour(x,y,z,linestyles=':')
ax2.set_title('b')
ax3 = fig.add_subplot(2,2,3)
ax3.contour(x,y,z,levels=np.arange(60,150,30),linewidths=5)
ax3.set_title('c')
ax4 = fig.add_subplot(2,2,4)
contour = ax4.contour(x,y,z,levels=np.arange(60,150,30))
ax4.clabel(contour,fontsize=10,colors='k')
ax4.set_title('d')
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
import pandas as pd
#读取数据
f = xr.open_dataset('/home/mw/input/pythonbook9857/data.nc')
z = f['hgt'].loc['2005-01-01':'2005-12-01',500,:,120]
time = z.time
lat = z.lat
#创建Figure
fig = plt.figure(figsize=(16, 13))
#绘制500hPa位势高度场
ax1 = fig.add_subplot(1,1,1)
#绘制等值线
q1 = ax1.contour(range(time.shape[0]), lat, z.T)
#设置y轴为对数坐标轴，并设置相关标签
ax1.set_xticks(np.arange(0,12,1))
ax1.set_yticks(np.arange(-90,120,30))
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#设置x轴标签
ax1.set_xlim(0,11)
ax1.set_xticks(np.arange(0,12,1))
ax1.set_xticklabels(pd.date_range(start='2005-01',periods=12,freq='M').date )
#添加图题
ax1.set_title('2005 500hPa Z',loc='center',fontsize=18)
plt.show()

9.5.2 带有地图投影的等值线图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
#读取数据
f = xr.open_dataset('/home/mw/input/pythonbook9857/data.nc')
z = f['hgt'].loc['2005-07-01',500,:,:]
t = f['air'].loc['2005-07-01',700,:,:]
lat = f['lat']
lon = f['lon']
#创建Figure
fig = plt.figure(figsize=(15, 12))
#绘制500hPa位势高度场
ax1 = fig.add_subplot(1,2,1, projection=ccrs.PlateCarree())
#设置ax1的范围
ax1.set_extent([60,180,0,90])
#为ax1添加海岸线
#ax1.coastlines()
#为ax1添加地理经纬度标签及刻度
ax1.set_xticks(np.arange(60,210,30), crs=ccrs.PlateCarree())
ax1.set_yticks(np.arange(0,90,30), crs=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#从5000到5800每间隔80绘制一条黑色的等值线
c1 = ax1.contour(lon, lat, z,levels=np.arange(5000,5880,80), colors='k',transform=ccrs.PlateCarree())
#为c1添加数值标签，标签格式为整型
ax1.clabel(c1, fmt='%d', fontsize=14.5,inline=True)
#单独绘制加粗并且用虚线表示的5880副高特征线
c2 = ax1.contour(lon, lat, z,levels = [5880],colors='k',transform=ccrs.PlateCarree(),linewidths=3,linestyles='--')
#为c2添加数值标签，标签格式为浮点型（保留一位小数）
ax1.clabel(c2, fmt='%d', fontsize=14.5,inline=True)
#添加图题
ax1.set_title('July 2005 500hPa Z',loc='center',fontsize=18)
ax1.set_title('unit: gpm',loc='right',fontsize=18)
#绘制700hPa温度场
ax2 = fig.add_subplot(1,2,2, projection=ccrs.PlateCarree())
ax2.set_extent([60,180,0,90])
#ax2.coastlines()
ax2.set_xticks(np.arange(60,210,30), crs=ccrs.PlateCarree())
ax2.set_yticks(np.arange(0,90,30), crs=ccrs.PlateCarree())
ax2.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax2.yaxis.set_major_formatter(cticker.LatitudeFormatter())
c3 = ax2.contour(lon, lat, t, levels=np.arange(-30,33,3),colors='k',transform=ccrs.PlateCarree())
ax2.clabel(c3, fmt='%.1f', fontsize=14.5,inline=True)
ax2.set_title('July 2005 700hPa T',loc='center',fontsize=18)
ax2.set_title('unit: $^\circ$C',loc='right',fontsize=18)
plt.show()

9.5.3 垂直剖面等值线图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
import pandas as pd
import numpy as np
#读取数据
f = xr.open_dataset('/home/mw/input/pythonbook9857/data.nc')
t = f['air'].loc['2005-07-01',:,:,80:150].mean('lon')
lev = t.level
lat = t.lat
#创建Figure
fig = plt.figure(figsize=(10, 8))
#绘制500hPa位势高度场
ax1 = fig.add_subplot(1,1,1)
#绘制温度场垂直剖面
c1 = ax1.contour(lat, lev, t,colors='k')
ax1.clabel(c1, fmt='%d', fontsize=14.5,inline=True)
#设置y轴为对数坐标轴，并设置相关标签
ax1.set_yscale('symlog')
ax1.set_xticks(np.arange(0,12,1))
ax1.set_yticks([1000,850,700,500,300,200,100])
ax1.set_yticklabels([1000,850,700,500,300,200,100])
ax1.set_ylim(1000,100)
#为ax1的x轴添加地理经纬度标签及刻度
ax1.set_xticks(np.arange(-90,120,30))
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
#添加图题
ax1.set_title('July 2005 T',loc='center',fontsize=18)

Text(0.5, 1.0, 'July 2005 T')

9.6 填色图

9.6.1 contourf()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
import numpy as np
import matplotlib.pyplot as plt
#生成数据
x = np.arange(-10,10,0.01)
y = np.arange(-10,10,0.01)
x_grid,y_grid = np.meshgrid(x,y)
z = x_grid**2+y_grid**2
#创建Figure及Axes
fig1 = plt.figure(figsize=(10, 8))
#基础填色图
ax1 = fig1.add_subplot(2,2,1)
ca = ax1.contourf(x,y,z)
ax1.set_title('a')
fig1.colorbar(ca,ax=ax1)
#填充等值线图
ax2 = fig1.add_subplot(2,2,2)
cb1 = ax2.contour(x,y,z,colors="k")
cb2 = ax2.contourf(x,y,z,hatches=['.',',,','-', '/', '\\', '//'],colors="none")
ax2.set_title('b')
fig1.colorbar(ca,ax=ax2)
#特定level值的填色图1
ax3 = fig1.add_subplot(2,2,3)
cc = ax3.contourf(x,y,z,levels=np.arange(60,150,30))
ax3.set_title('c')
fig1.colorbar(cc,ax=ax3)
#特定level值的填色图2
ax4 = fig1.add_subplot(2,2,4)
cd = contour = ax4.contourf(x,y,z,levels=np.arange(60,150,30),extend='both')
ax4.set_title('d')
fig1.colorbar(cd,ax=ax4)
plt.show()

/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:18: MatplotlibDeprecationWarning: hatch must consist of a string of "*+-./OX\ox|" or None, but found the following invalid values ",". Passing invalid values is deprecated since 3.4 and will become an error in 3.7.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
#读取数据
f = xr.open_dataset('/home/mw/input/pythonbook9857/data.nc')
z = f['hgt'].loc['2005-07-01',500,:,:]
t = f['air'].loc['2005-07-01',700,:,:]
lat = f['lat']
lon = f['lon']
#创建Figure
fig = plt.figure(figsize=(15, 12))
#绘制500hPa位势高度场
ax1 = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
#设置ax1的范围
ax1.set_extent([60,180,0,90])
#为ax1添加海岸线
#ax1.coastlines()
#为ax1添加地理经纬度标签及刻度
ax1.set_xticks(np.arange(60,210,30), crs=ccrs.PlateCarree())
ax1.set_yticks(np.arange(0,90,30), crs=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#将700hPa温度场绘制为填色等值线图
c1 = ax1.contourf(lon, lat, t,levels=np.arange(-21,22,1),extend='both',cmap='bwr',zorder=1,transform=ccrs.PlateCarree())
#把温度大于(小于)10(-10)℃的区域覆盖上阴影，100(-100)可以换为更大或更小的数，只是用于表明区间
c2 = ax1.contourf(lon, lat, t,levels=[-100,-10,10,100],hatches=['//', None,'//'],zorder=2,colors='none',transform=ccrs.PlateCarree())
#从5000到5880每间隔80绘制一条黑色的等值线
c3 = ax1.contour(lon, lat, z,levels=np.arange(5000,5880,80), colors='k',zorder=3,transform=ccrs.PlateCarree())
#为c1添加数值标签，标签格式为整型
ax1.clabel(c3, fmt='%d', fontsize=14.5,inline=True)
#单独绘制加粗并且用虚线表示的5880副高特征线
c4 = ax1.contour(lon, lat, z,levels = [5880],colors='k',zorder=4,transform=ccrs.PlateCarree(),linewidths=3,linestyles='--')
#为c2添加数值标签，标签格式为浮点型（保留一位小数）
ax1.clabel(c4, fmt='%d', fontsize=14.5,inline=True)
#添加图题
ax1.set_title('July 2005 500hPa Z&T',loc='center',fontsize=18)
#添加色标
fig.colorbar(c1,ax=ax1,fraction=0.032)
plt.show()

9.6.2 pcolor()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
#读取数据
f = xr.open_dataset('/home/mw/input/pythonbook9857/data.nc')
z = f['hgt'].loc['2005-07-01',500,:,:]
t = f['air'].loc['2005-07-01',700,:,:]
lat = f['lat']
lon = f['lon']
#创建Figure
fig = plt.figure(figsize=(15, 12))
#绘制500hPa位势高度场
ax1 = fig.add_subplot(2,2,1, projection=ccrs.PlateCarree())
#设置ax1的范围
ax1.set_extent([60,180,0,90])
#为ax1添加海岸线
#ax1.coastlines()
#为ax1添加地理经纬度标签及刻度
ax1.set_xticks(np.arange(60,210,30), crs=ccrs.PlateCarree())
ax1.set_yticks(np.arange(0,90,30), crs=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#将700hPa温度场绘制为填色等值线图
c1 = ax1.pcolor(lon, lat, t,vmin=-20, vmax=20,cmap='bwr',zorder=1,transform=ccrs.PlateCarree())
#添加图题
ax1.set_title('(a) July 2005 500hPa T',loc='center',fontsize=18)
#添加色标
fig.colorbar(c1,ax=ax1,fraction=0.032)
#ax2
ax2 = fig.add_subplot(2,2,2, projection=ccrs.PlateCarree())
#设置ax2的范围
ax2.set_extent([60,180,0,90])
#为ax2添加海岸线
#ax2.coastlines()
#为ax2添加地理经纬度标签及刻度
ax2.set_xticks(np.arange(60,210,30), crs=ccrs.PlateCarree())
ax2.set_yticks(np.arange(0,90,30), crs=ccrs.PlateCarree())
ax2.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax2.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#将700hPa温度场绘制为填色等值线图
c2 = ax2.pcolor(lon, lat, t, vmin=-20, vmax=20, cmap='bwr', zorder=1, edgecolors='k', transform=ccrs.PlateCarree())
#添加图题
ax2.set_title('(b) July 2005 500hPa T',loc='center',fontsize=18)
#添加色标
fig.colorbar(c2,ax=ax2,fraction=0.032)
#创建一个t_copy数组并将大于15的数据设置为缺测
t_copy = t.values.copy()
t_copy[t_copy>15] = np.nan
#ax3
ax3 = fig.add_subplot(2,2,3, projection=ccrs.PlateCarree())
#设置ax3的范围
ax3.set_extent([60,180,0,90])
#为ax3添加海岸线
#ax3.coastlines()
#为ax3添加地理经纬度标签及刻度
ax3.set_xticks(np.arange(60,210,30), crs=ccrs.PlateCarree())
ax3.set_yticks(np.arange(0,90,30), crs=ccrs.PlateCarree())
ax3.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax3.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#将700hPa温度场绘制为填色等值线图
c3 = ax3.contourf(lon, lat, t_copy,levels=np.arange(-21,22,1), extend='both', cmap='bwr', zorder=1, transform=ccrs.PlateCarree())
#添加图题
ax3.set_title('(c) July 2005 500hPa T',loc='center',fontsize=18)
#添加色标
fig.colorbar(c3,ax=ax3,fraction=0.032)
#ax4
ax4 = fig.add_subplot(2,2,4, projection=ccrs.PlateCarree())
#设置ax4的范围
ax4.set_extent([60,180,0,90])
#为ax4添加海岸线
#ax4.coastlines()
#为ax4添加地理经纬度标签及刻度
ax4.set_xticks(np.arange(60,210,30), crs=ccrs.PlateCarree())
ax4.set_yticks(np.arange(0,90,30), crs=ccrs.PlateCarree())
ax4.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax4.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#将700hPa温度场绘制为填色等值线图
c4 = ax4.pcolor(lon, lat, t_copy,vmin=-20, vmax=20,cmap='bwr',zorder=1,edgecolors='k',transform=ccrs.PlateCarree())
#添加图题
ax4.set_title('(d) July 2005 500hPa T',loc='center',fontsize=18)
#添加色标
fig.colorbar(c4,ax=ax4,fraction=0.032)
plt.show()

9.7 轨迹绘制（以台风路径的绘制为例）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
import os
import pandas as pd
import numpy as np
from pathlib import Path
from typing import List
from typing import Union
from typing import Tuple
from matplotlib.collections import LineCollection
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
#读取CMA热带气旋最佳路径数据集
def reader(
   typhoon_txt: os.PathLike, code: Union[str, int]
) -> Tuple[List[str], pd.DataFrame]:
   typhoon_txt = Path(typhoon_txt)
   if isinstance(code, int):
      code = "{:04}".format(code)
   with open(typhoon_txt, "r") as txt_handle:
      while True:
         header = txt_handle.readline().split()
         if not header:
            raise ValueError(f"没有在文件里找到编号为{code}的台风的数据")
         if header[4].strip() == code:
            break
         [txt_handle.readline() for _ in range(int(header[2]))]
      data_path = pd.read_table(
         txt_handle,
         sep=r"\s+",
         header=None,
         names=["TIME", "I", "LAT", "LONG", "PRES", "WND", "OWD"],
         nrows=int(header[2]),
         dtype={
            "I": np.int,
            "LAT": np.float32,
            "LONG": np.float32,
            "PRES": np.float32,
            "WND": np.float32,
            "OWD": np.float32,
         },
         parse_dates=True,
         date_parser=lambda x: pd.to_datetime(x, format=f"%Y%m%d%H"),
         index_col="TIME",
      )
      data_path["LAT"] = data_path["LAT"] / 10
      data_path["LONG"] = data_path["LONG"] / 10
      return header, data_path
#读取0608号台风的数据（数据为超6小时的占风相关信息）
head, dat = reader(r"/home/mw/input/pythonbook9857/CH2006BST.txt",'0608')
lat = dat.LAT
lon = dat.LONG
level = dat.I
pressure = dat.PRES
#创建Figure
fig = plt.figure(figsize=(15, 12))
#绘制台风路径
ax1 = fig.add_subplot(1,2,1, projection=ccrs.PlateCarree())
#设置ax1的范围
ax1.set_extent([100,160,-10,40])
#为ax1添加海岸线
#ax1.coastlines()
#ax1.add_feature(cfeature.LAND) #添加大陆特征
#为ax1添加地理经纬度标签及刻度
ax1.set_xticks(np.arange(100,170,10), crs=ccrs.PlateCarree())
ax1.set_yticks(np.arange(-10,50,10), crs=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#绘制台风路径，并将标记逐6小时坐标点及其对应的台风强度
ax1.plot(lon,lat,linewidth=2)
s1 = ax1.scatter(lon,lat,c=pressure,s=(level+1)*13,cmap='Reds_r',vmax=1050,vmin=900,alpha=1)
fig.colorbar(s1,ax=ax1,fraction=0.04)
#绘制台风路径
ax2 = fig.add_subplot(1,2,2, projection=ccrs.PlateCarree())
#设置ax2的范围
ax2.set_extent([100,160,-10,40])
#为ax2添加海岸线
#ax2.coastlines()
#ax2.add_feature(cfeature.LAND) #添加大陆特征
#为ax2添加地理经纬度标签及刻度
ax2.set_xticks(np.arange(100,170,10), crs=ccrs.PlateCarree())
ax2.set_yticks(np.arange(-10,50,10), crs=ccrs.PlateCarree())
ax2.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax2.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#将经纬度数据存入同一数组
points = np.array([lon, lat]).T.reshape(-1, 1, 2)
segments = np.concatenate([points[:-1], points[1:]], axis=1)
#设置色标的标准化范围(即将Z维度的数据对应为颜色数组，Z维度指dat.WND[-1])
norm = plt.Normalize(0, 80)
#设置颜色线条
lc = LineCollection(segments, cmap='jet', norm=norm,transform=ccrs.PlateCarree())
lc.set_array(dat.WND[:-1])
#绘制线条
line = ax2.add_collection(lc)   
fig.colorbar(lc,ax=ax2,fraction=0.04)
plt.show()

/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:34: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations

9.8 流线图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
#读取数据
f = xr.open_dataset('/home/mw/input/pythonbook9857/data.nc')
u = f['uwnd'].loc['2005-07-01',500,:,:].values
v = f['vwnd'].loc['2005-07-01',500,:,:].values
lat = f['lat']
lon = f['lon']
#创建Figure
fig = plt.figure(figsize=(16, 13))
#绘制500hPa位势高度场
ax1 = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
#设置ax1的范围
ax1.set_extent([80,140,0,40])
#为ax1添加海岸线
#ax1.coastlines()
#为ax1添加地理经纬度标签及刻度
ax1.set_xticks(np.arange(80,160,20), crs=ccrs.PlateCarree())
ax1.set_yticks(np.arange(0,60,20), crs=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#绘制500hPa风矢量流线图
q1 = ax1.streamplot(lon, lat, u, v,density=[1, 5],transform=ccrs.PlateCarree())
plt.show()

9.9 矢量箭头图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(15, 3))
ax1 = fig.add_subplot(1,2,1)
# 默认
q1 = ax1.quiver(0, 0, 1, 1)
# scale：缩放比例，值越大，箭头越短
q2 = ax1.quiver(0.1, 0, 1, 1, scale=15)
# width：主轴宽度
q3 = ax1.quiver(0.2, 0, 1, 1, width=0.05)
# headwidth：矢量符号箭头头部的宽度
q4 = ax1.quiver(0.3, 0, 1, 1, headwidth=8)
# headlength：矢量符号箭头头部的长度，从头部到尾部
q5 = ax1.quiver(0.4, 0, 1, 1, headlength=8)
# headaxislength 矢量符号箭头的头部到矢量符号轴和箭头的交接处的长度
q6 = ax1.quiver(0.5, 0, 1, 1, headwidth=8, headlength=8, headaxislength=2)
# minlength 最小箭头长度，当矢量长度小于此值时，矢量将被替换为正六边形
q7 = ax1.quiver(0.65, 0, 1, 1, headwidth=4, headlength=4, minlength=10)
ax1.set_axis_off()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
#读取数据
f = xr.open_dataset('/home/mw/input/pythonbook9857/data.nc')
u = f['uwnd'].loc['2005-07-01',500,:,:].values
v = f['vwnd'].loc['2005-07-01',500,:,:].values
lat = f['lat']
lon = f['lon']
#创建Figure
fig = plt.figure(figsize=(16, 13))
#绘制500hPa位势高度场
ax1 = fig.add_subplot(1,2,1, projection=ccrs.Mercator())
#设置ax1的范围
ax1.set_extent([80,140,0,40])
#为ax1添加海岸线
ax1.coastlines()
#为ax1添加地理经纬度标签及刻度
ax1.set_xticks(np.arange(80,160,20), crs=ccrs.PlateCarree())
ax1.set_yticks(np.arange(0,60,20), crs=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#绘制500hPa水平风场
q1 = ax1.quiver(lon, lat, u,v,scale=200,transform=ccrs.PlateCarree())
#添加参考矢量
ax1.quiverkey(q1,0.9,1.01,U=30,label='30m/s')
#添加图题
ax1.set_title('(a) July 2005 500hPa UV',loc='center',fontsize=18)
#绘制ax2
ax2 = fig.add_subplot(1,2,2, projection=ccrs.Mercator())
#设置ax2的范围
ax2.set_extent([80,140,0,40])
#为ax2添加海岸线
ax2.coastlines()
#为ax2添加地理经纬度标签及刻度
ax2.set_xticks(np.arange(80,160,20), crs=ccrs.PlateCarree())
ax2.set_yticks(np.arange(0,60,20), crs=ccrs.PlateCarree())
ax2.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax2.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#绘制500hPa水平风场，[::2,::2]表示每隔两个栅格绘制一次，避免箭头过于密集
q1 = ax2.quiver(lon[::2], lat[::2], u[::2,::2], v[::2,::2], scale=200, transform=ccrs.PlateCarree())
#添加参考矢量
ax1.quiverkey(q1,0.9,1.01,U=30,label='30m/s')
#添加图题
ax2.set_title('(b) July 2005 500hPa UV',loc='center',fontsize=18)
plt.show()

/opt/conda/lib/python3.7/site-packages/cartopy/io/__init__.py:260: DownloadWarning: Downloading: https://naciscdn.org/naturalearth/50m/physical/ne_50m_coastline.zip
  warnings.warn('Downloading: {}'.format(url), DownloadWarning)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import pandas as pd
#读取数据
f = xr.open_dataset('/home/mw/input/pythonbook9857/data.nc')
u = f['uwnd'].loc['2005-01-01':'2005-12-01',:,35,120]
w = f['omega'].loc['2005-01-01':'2005-12-01',:,35,120]
time = u.time
lev = u.level
#创建Figure
fig = plt.figure(figsize=(16, 13))
#绘制500hPa位势高度场
ax1 = fig.add_subplot(1,1,1)
#绘制500hPa水平风场
q1 = ax1.quiver(range(time.shape[0]), lev, u,w*200,pivot='mid')
#添加参考矢量
ax1.quiverkey(q1,0.9,1.01,U=30,label='30m/s')
#设置y轴为对数坐标轴，并设置相关标签
ax1.set_yscale('symlog')
ax1.set_xticks(np.arange(0,12,1))
ax1.set_yticks([1000,850,700,500,300,200,100])
ax1.set_yticklabels([1000,850,700,500,300,200,100])
ax1.set_ylim(1050,90)
#设置x轴标签
ax1.set_xlim(-1,12)
ax1.set_xticks(np.arange(0,12,1))
ax1.set_xticklabels(pd.date_range(start='2005-01',periods=12,freq='M').date )
#添加图题
ax1.set_title('2005 U&OMEGA',loc='center',fontsize=18)
plt.show()

9.10 风向杆图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
#读取数据
f = xr.open_dataset('/home/mw/input/pythonbook9857/data.nc')
u = f['uwnd'].loc['2005-07-01',500,:,:].values
v = f['vwnd'].loc['2005-07-01',500,:,:].values
lat = f['lat']
lon = f['lon']
#创建Figure
fig = plt.figure(figsize=(16, 13))
#绘制500hPa位势高度场
ax1 = fig.add_subplot(1,1,1, projection=ccrs.Mercator())
#设置ax1的范围
ax1.set_extent([80,140,0,40])
#为ax1添加海岸线
#ax1.coastlines()
#为ax1添加地理经纬度标签及刻度
ax1.set_xticks(np.arange(80,160,20), crs=ccrs.PlateCarree())
ax1.set_yticks(np.arange(0,60,20), crs=ccrs.PlateCarree())
ax1.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#绘制500hPa水平风场
#Matplotlib原始的风向杆等级划分不适用于国内标准，所以需要额外指定符合标准的风向杆等级
barb_increments = {'half': 2, 'full': 4, 'flag': 20}
q1 = ax1.barbs(lon, lat, u,v,barb_increments=barb_increments,transform=ccrs.PlateCarree())
#添加图题
ax1.set_title('(a) July 2005 500hPa Wind',loc='center',fontsize=18)
plt.show()

9.11 探空图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
import numpy as np
# MetPy中用于绘制探空图的类
from metpy.plots import SkewT
from metpy.units import units
from metpy.calc import lcl
from metpy.calc import parcel_profile
# 探空数据，原始的np.ndarray数组需要带上单位，转换为带单位的数据
# 气压层
p = np.array([1000, 925, 850, 700, 600, 500, 450, 400, 300, 250]) * units.hPa
# 气压层对应温度
t = np.array([4, 8, 3, -11, -21, -26, -33, -38, -55, -60]) * units.degC
# 气压层对应露点
td = np.array([-8, -9, -14, -18, -25, -34, -38, -43, -61, -67]) * units.degC
# 气压层对应风的U分量
u = np.array([-0.39, 0.11, 3.1, 10.7, 16.61, 24.0, 20.31, 33.43, 49.32, 59.21]) * units('m/s')
# 气压层对应风的V分量
v = np.array([-0.57, -0.75, -1.09, -0.79, -0.48, -0.04, -0.26, 0.96, 2.87, 4.14]) * units('m/s')
# 计算气块绝热抬升参数
prof = parcel_profile(p, t[0], td[0]).to('degC')
# 用最底层气压(这里是1000hPa)、温度和露点，计算抬升凝结高度对应的气压和温度
lcl_pressure, lcl_temperature = lcl(p[0], t[0], td[0])
fig = plt.figure(figsize=(8, 7))
# 用于绘制探空图的绘图实例
skew = SkewT(fig)
# 绘制气块绝热抬升路径
skew.plot(p, prof, 'k')
# 在图上标出抬升凝结高度和对应温度所在的点
skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')
# 绘制CIN阴影部分
skew.shade_cin(p, t, prof)
# 绘制CAPE阴影部分
skew.shade_cape(p, t, prof)
# 绘制0℃等温线
skew.ax.axvline(0, color='c', linestyle='--', linewidth=1)
# 环境温度垂直廓线
skew.plot(p, t, 'r')
# 环境露点垂直廓线
skew.plot(p, td, 'g')
# 高度层对应水平风场(原始风向杆等级划分不适用于国内标准)
barb_increments = {'half': 2, 'full': 4, 'flag': 20}
skew.plot_barbs(p, u, v, barb_increments=barb_increments)
# 绘制干绝热线
skew.plot_dry_adiabats(t0=np.arange(233, 533, 10) * units.K, alpha=0.5, color='orangered', linewidth=0.7)
# 绘制湿绝热线
skew.plot_moist_adiabats(t0=np.arange(233, 400, 5) * units.K, alpha=0.5, color='tab:green', linewidth=0.7)
# 绘制混合比线
skew.plot_mixing_lines(pressure=np.arange(1000, 99, -20) * units.hPa, linestyle='dotted', color='tab:blue', linewidth=0.7)
#设置Y轴(高度层)范围
skew.ax.set_ylim(1000, 250)
#设置X轴(温度)范围
skew.ax.set_xlim(-40, 50)
plt.show()

9.12 泰勒图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
from matplotlib.projections import PolarAxes
from mpl_toolkits.axisartist import floating_axes
from mpl_toolkits.axisartist import grid_finder
import numpy as np
import matplotlib.pyplot as plt
#绘制泰勒图坐标系
def set_tayloraxes(fig, location):
#新建极坐标系
    trans = PolarAxes.PolarTransform()
#相关系数轴
    r1_locs = np.hstack((np.arange(1,10)/10.0,[0.95,0.99]))
    t1_locs = np.arccos(r1_locs)        
    gl1 = grid_finder.FixedLocator(t1_locs)    
    tf1 = grid_finder.DictFormatter(dict(zip(t1_locs, map(str,r1_locs))))
#标准差轴
    r2_locs = np.arange(0,2,0.25)
    r2_labels = ['0 ', '0.25 ', '0.50 ', '0.75 ', 'REF ', '1.25 ', '1.50 ', '1.75 ']
    gl2 = grid_finder.FixedLocator(r2_locs)
    tf2 = grid_finder.DictFormatter(dict(zip(r2_locs, map(str,r2_labels))))
    ghelper = floating_axes.GridHelperCurveLinear(trans,extremes=(0,np.pi/2,0,1.75),
                                                  grid_locator1=gl1,tick_formatter1=tf1,
                                                  grid_locator2=gl2,tick_formatter2=tf2)
    ax = floating_axes.FloatingSubplot(fig, location, grid_helper=ghelper)
    fig.add_subplot(ax)
#设置各个轴的格式
    ax.axis["top"].set_axis_direction("bottom")  
    ax.axis["top"].toggle(ticklabels=True, label=True)
    ax.axis["top"].major_ticklabels.set_axis_direction("top")
    ax.axis["top"].label.set_axis_direction("top")
    ax.axis["top"].label.set_text("Correlation")
    ax.axis["top"].label.set_fontsize(14)
    ax.axis["left"].set_axis_direction("bottom") 
    ax.axis["left"].label.set_text("Standard deviation")
    ax.axis["left"].label.set_fontsize(14)
    ax.axis["right"].set_axis_direction("top")   
    ax.axis["right"].toggle(ticklabels=True)
    ax.axis["right"].major_ticklabels.set_axis_direction("left")
    ax.axis["bottom"].set_visible(False)         
    ax.grid(True)
    polar_ax = ax.get_aux_axes(trans)   

    rs,ts = np.meshgrid(np.linspace(0,1.75,100),
                            np.linspace(0,np.pi/2,100))
    rms = np.sqrt(1 + rs**2 - 2*rs*np.cos(ts))
    CS = polar_ax.contour(ts, rs,rms,colors='gray',linestyles='--')
    plt.clabel(CS, inline=1, fontsize=10)
    t = np.linspace(0,np.pi/2)
    r = np.zeros_like(t) + 1
    polar_ax.plot(t,r,'k--')
#将垂直轴的REF改为1.00
    polar_ax.text(np.pi/2+0.032,1.02, " 1.00", size=10.3,ha="right", va="top",
                  bbox=dict(boxstyle="square",ec='w',fc='w'))
    return polar_ax
#在泰勒图上绘制数据点
def plot_taylor(axes, refsample, sample, *args, **kwargs):
    std = np.std(refsample)/np.std(sample)
    corr = np.corrcoef(refsample, sample) 
    theta = np.arccos(corr[0,1])
    t,r = theta,std
    d = axes.plot(t,r, *args, **kwargs) 
    return d
1
2
3
4
5
6
7
8
9
10
11
12
13
14
x = np.linspace(0,100*np.pi,100)
#生成观测数据
data = np.sin(x)    
#生成3组与原始数据形状相同的对比数据                 
m1 = data + 0.4*np.random.randn(len(x))   
m2 = 0.3*data + 0.6*np.random.randn(len(x)) 
m3 = np.sin(x-np.pi/10)           
#绘图     
fig = plt.figure(figsize=(10,4))
ax1 = set_tayloraxes(fig, 121)
d1 = plot_taylor(ax1,data,m1, 'bo')
d2 = plot_taylor(ax1,data,m2, 'ro')
d3 = plot_taylor(ax1,data,m3, 'go')
plt.show()
/小书匠/python