一、知识点

读取NC格式在分析资料里面需要的数据，画填色图和风羽图

小知识点：
1.数据的读取
2.数据的清洗
3.风羽图的设置

二、代码拆分

导入包

#处理数据的包
import xarray as xr
import numpy as np
#画图的包
import matplotlib.pyplot as plt
#地图的包
import cartopy.crs as ccrs
import my_class   #一个自己写的画地图的py文件
1
2
3
4
5
6
7
8

读取数据

#打开nc文件
ds = xr.open_dataset('download_201305.nc')
# print(ds)
#获取经度，纬度，时间，时间的作用下面讲
lat = ds.latitude#读取维度
lon = ds.longitude#读取经度
u = ds['u']#风场U分量
v = ds['v']#风场V分量
vo = ds['vo']#涡度
1
2
3
4
5
6
7
8
9

注意：这里读取的数据是全部的格点数据，但是我们画图用不了这么多，所以对数据做分割。只需要一部分

数据分割

#设置经纬度范围，我们并不需要画出多余的数据
lat_range = lat[(lat>=22) & (lat<=40)]
lon_range = lon[(lon>=102) & (lon<=126)]
#根据时间，经纬度得到数据，U,V,涡度
u_region = u.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-27T00:00:00.000000000')
v_region = v.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-27T00:00:00.000000000')
vo_region = vo.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-27T00:00:00.000000000')
#涡度乘以100000，因为是10的负5次方，方便呈现效果
vo_region = vo_region*100000
1
2
3
4
5
6
7
8
9

这里涡度乘以了100000，是为了方便观察，让数据绘制出来的效果更好。

数据处理

#下面是数据的简单清洗，本来没有这个，直接填色，发现有很多负值，画图太难看了
#设置条件，找到涡度大于等于2的位置
e = (vo_region >= 2)
#利用np包，把所有不满足条件的位置改写成0
vo_region = np.where(e,vo_region,0)
1
2
3
4
5

把涡度小于2的全部置零，这样在填色的时候小于2的地方就白了。

画图设置

#设置画布15*7
fig = plt.figure(figsize=(15,7))
#添加子图，如果是一个就是111，一行一列第一个，设置投影方式，因为要画地图了
ax = fig.add_subplot(111,projection = ccrs.PlateCarree())
#设置x,y轴标签
ax.set_xticks(np.arange(102,127,2),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(22,41,2),crs=ccrs.PlateCarree())
xticks_str = ['102', '', '106', '', '110', '', '114', '', '118','', '122', '', '126°E']
ax.set_xticklabels(xticks_str,fontsize = 18)
yticks_str = ['22   ','24    ','26    ','28    ','30    ','32    ','34    ','36    ','38    ','40°N']
ax.set_yticklabels(yticks_str,fontsize = 20)
#画地图
my_class.readshapefile('bou2_4l.shp',linewidth=1,ax=ax)
1
2
3
4
5
6
7
8
9
10
11
12
13

画图

#画填色图
cf_rh = ax.contourf(lon_range,#lon
                    lat_range ,#lat
                    vo_region,#vo
                    #挑了4款效果相对较好的色条，可以换着看看效果
                    cmap='hot_r',#加入色条   #hot_r     lGnBu    OrRd    YlGnBu
                    extend='both')#设置色条为两头尖
#色标的设置
cb=fig.colorbar(cf_rh,orientation='vertical',shrink=0.75)


#画风杆
ax.barbs(lon_range[::4],#纬度降尺度，0.25变成1
         lat_range[::4],#经度降尺度，0.25变成1
         u_region[::4,::4],#u降尺度，0.25变成1
         v_region[::4,::4]#v降尺度，0.25变成1
         ,barbcolor=['k'],#颜色
         linewidth=0.5, #宽度
         length=5, #长度
         barb_increments=dict(half=2, full=4, flag=20))#风羽图设置   半根2米，一根4米，三角20米



#写a
ax.text(0.05, 0.95, '(a)',transform=ax.transAxes, fontsize=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

各种画图参数都在代码中注释好了。

出图

plt.show()
1

三、完整代码

#处理数据的包
import xarray as xr
import numpy as np
#画图的包
import matplotlib.pyplot as plt
#地图的包
import cartopy.crs as ccrs
import my_class

#打开nc文件
ds = xr.open_dataset('download_201305.nc')
# print(ds)
#获取经度，纬度，时间，时间的作用下面讲
lat = ds.latitude
lon = ds.longitude
time = ds.time
u = ds['u']
v = ds['v']
vo = ds['vo']
'2013-06-27T00:00:00.000000000'
#设置经纬度范围，我们并不需要画出多余的数据
lat_range = lat[(lat>=22) & (lat<=40)]
lon_range = lon[(lon>=102) & (lon<=126)]


#根据时间，经纬度得到数据，U,V,涡度
u_region = u.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-27T00:00:00.000000000')
v_region = v.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-27T00:00:00.000000000')
vo_region = vo.sel(longitude=lon_range, latitude=lat_range,time = '2013-06-27T00:00:00.000000000')
#涡度乘以100000，因为是10的负5次方，方便呈现效果
vo_region = vo_region*100000
#设置画布15*7
fig = plt.figure(figsize=(15,7))
#添加子图，如果是一个就是111，一行一列第一个，设置投影方式，因为要画地图了
ax = fig.add_subplot(111,projection = ccrs.PlateCarree())



#设置x,y轴标签
ax.set_xticks(np.arange(102,127,2),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(22,41,2),crs=ccrs.PlateCarree())
xticks_str = ['102', '', '106', '', '110', '', '114', '', '118','', '122', '', '126°E']
ax.set_xticklabels(xticks_str,fontsize = 18)
yticks_str = ['22   ','24    ','26    ','28    ','30    ','32    ','34    ','36    ','38    ','40°N']
ax.set_yticklabels(yticks_str,fontsize = 20)
my_class.readshapefile('bou2_4l.shp',linewidth=1,ax=ax)





#下面是数据的简单清洗，本来没有这个，直接填色，发现有很多负值，画图太难看了
#设置条件，找到涡度大于等于2的位置
e = (vo_region >= 2)
#利用np包，把所有不满足条件的位置改写成0
vo_region = np.where(e,vo_region,0)


#画填色图
cf_rh = ax.contourf(lon_range,#lon
                    lat_range ,#lat
                    vo_region,#vo
                    #用了系统自带的色条，文章中的色条与你给的数据不匹配，效果不好
                    #挑了4款效果相对较好的色条，可以换着看看效果
                    cmap='hot_r',#加入色条   #hot_r     lGnBu    OrRd    YlGnBu
                    extend='both')#设置色条为两头尖
#色标的设置
cb=fig.colorbar(cf_rh,orientation='vertical',shrink=0.75)


#画风杆
ax.barbs(lon_range[::4],#纬度降尺度，0.25变成1
         lat_range[::4],#经度降尺度，0.25变成1
         u_region[::4,::4],#u降尺度，0.25变成1
         v_region[::4,::4]#v降尺度，0.25变成1
         ,barbcolor=['k'],#颜色
         linewidth=0.5, #宽度
         length=5, #长度
         barb_increments=dict(half=2, full=4, flag=20))#风羽图设置   半根2米，一根4米，三角20米



#写a
ax.text(0.05, 0.95, '(a)',transform=ax.transAxes, fontsize=11)

#出图
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
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