读取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文件 

#打开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']#涡度 

注：这里读取的数据是所有的格点数据，但我们不需要那么多的绘图，所以我们分割数据。只需要一部分

#设置经纬度范围，不需要画多余的数据 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

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

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

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

#画填色图
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()

#处理数据的包
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()