首先,码上希望实现的成果样式:
这里,我处理的数据是平流强度数据,具体的实现思路是:
1、将数据分为5个bin
2、找出每个数据点对应的经纬度坐标
3、对每组数据进行循环画图,每次画图采用不同的颜色
这里主要需要用到几个库,以及函数
1、读取nc文件的xarray库
2、实现绘图的 matplotlib库
3、实现地形投影的cartopy库
4、实现数组计算的numpy库
5、np.argwhere()函数:获取数组中数据对应的索引值
可以查看官网说明:np.argwhere()
下面贴上代码以及每一步的说明:
#导入相关的库 import cartopy.feature as cfeature from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter import matplotlib.pyplot as plt import cartopy.crs as ccrs import numpy as np import xarray as xr #读取数据 path='D:\\mse.nc'#数据的路径位置 dh=xr.open_dataset(path)#读取数据 lon=np.array(dh['lon'])#读取数据中的经度并转为array数组 lat=np.array(dh['lat'])#读取数据中的纬度并转为array数组 time=dh['time']#读取数据中的时间 time=time.loc['1982':'2012'][:]#选择数据的时间范围 lat_range = lat[(lat>-22.5) & (lat<22.5)]#选择数据的纬度范围 hadv_region =dh.sel(lon=lon, lat=lat_range,time=time).hadv#读取选取时间、经度范围内的数据 hadv =np.array(hadv_region.mean('time', skipna=True))#对数据进行平均处理 hadv_range=np.arange(-125,125+50,50)#随机生成一组数值用于筛选数据 hadv_bin=[]#生成一个空的list fig=plt.figure(figsize=(20,12))#产生画板 plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号 plt.rcParams['font.sans-serif']=['Fangsong']#显示中文 ax=fig.add_subplot(111,projection=ccrs.PlateCarree(central_longitude=180))#设置投影以及投影的中心经度,创建一个画纸,创建子图 colors=['r','b','y','g','k']#定义五个bin的颜色 lab=['-100 MSE','-50 MSE','0 MSE','50 MSE','100 MSE']#定义五个bin的标签 # 循环挑选不同的bin ,每50为一个bin,并绘制散点图 for i in range(len(hadv_range)-1): idx=np.argwhere((hadv>hadv_range[i])&(hadv<hadv_range[i+1])) lonr=lon[idx[:,1]] latr=lat_range[idx[:,0]] ax.scatter(lonr,latr,marker='o',c=colors[i],\ transform=ccrs.PlateCarree(central_longitude=180),label=lab[i])#绘制散点图 ax.legend(loc='upper right', bbox_to_anchor=(1, 1.7))#设置图例以及位置 ax.coastlines()#添加海岸线 ax.set_xticks(np.arange(0, 360+45, 45),crs=ccrs.PlateCarree(central_longitude=180))#设置x轴的经度范围 ax.set_yticks(np.arange(-20, 30, 10),crs=ccrs.PlateCarree())#设置y轴的纬度范围 #设置刻度格式为经纬度格式 ax.xaxis.set_major_formatter(LongitudeFormatter())#设置刻度格式 ax.yaxis.set_major_formatter(LatitudeFormatter()) ax.set_title('热带海域 MSE 空间水平分布图',fontsize=20)#添加标题 ax.set_xlabel('经度($°$)',fontsize=20)#添加x轴标签 ax.set_ylabel('纬度($°$)',fontsize=20)#添加y轴标签 ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', \ edgecolor='black', facecolor='grey'))#添加陆地 # fig.savefig('D:\\picture\\'+'热带海域 MSE 空间水平分布图.tiff',format='tiff',dpi=150)#保存数据
结果如下,非常的surprise~,感兴趣的小伙伴赶快尝试吧
一个努力学习python的海洋小白 水平有限,欢迎指正!!! 欢迎评论、收藏。