最近，学习海气相互作用时，老师布置了一个小作业。通过卫星高度计测得的海表面高度异常数据，计算其表层地转流速，并研究其与海表面高度异常的关系。
刚刚看到作业内容时，还是有点一头雾水的，不晓得怎么通过海表面高度数据推算出地转流速。挠头思考时，无意中看到物理海洋学，想起里面曾经有关系地转流部分的内容，索性拿出来看看找一下。惊喜大于意外，书里果然有推导的公式，直接贴在下面了：
其中，δ \deltaδ为海表面高度，f为科氏力（f=2Ωsinψ）ψ纬度 ;，g为重力加速度,y x 就是经纬度。
根据公式很容易知道计算的过程，细节在于对于数据的处理，包括科氏力、经纬度转换等等。
# # 首先导入相关的库 import cartopy.feature as cfeature#陆地、河流等 from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter#经纬度 import matplotlib.pyplot as plt#绘图 import numpy as np #数组 import cartopy.crs as ccrs #投影 import xarray as xr#读取nc文件 # 读取nc文件数据 path='D:\\sla.nc' data=xr.open_dataset(path) sla=data.sla[0] .data#海表面高度异常数据 x=data.longitude.data #lon y=data.latitude.data #lat # 计算将弧度与角度转换 、科氏力、 a=6400000 omega=7.292e-5 g = 9.8 coslat=np.cos(y*np.pi/180).reshape((-1,1)) coslon=np.cos(x*np.pi/180).reshape((-1,1)) sinlat = (np.sin(y*np.pi/180)).reshape((y.shape[0],1)) f=2*omega*sinlat # np.gradient()用来求导 dx=(np.gradient(x)*np.pi/180).reshape((-1,1))*coslat*a dy=(np.gradient(y)*np.pi/180).reshape((-1,1))*a dslady=np.gradient(sla,axis=0)/dy dsladx=np.gradient(sla,axis=1)/dx #带入公式 us=-1*(g/f)*dslady vs=(g/f)*dsladx #==============================绘图================================= fig = plt.figure(figsize=(10, 9)) ax=fig.add_subplot(111,projection=ccrs.PlateCarree()) plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号 ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', \ edgecolor='white', facecolor='white',zorder=2)) # # # 设置经纬度范围、显示区域 ax.set_xticks(np.arange(143, 146, 0.5),crs=ccrs.PlateCarree()) ax.set_yticks(np.arange(33, 36, 0.5),crs=ccrs.PlateCarree()) ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label =False))#经度0不加标识 ax.yaxis.set_major_formatter(LatitudeFormatter()) ax.set_extent([143,145,33,35]) # # #绘制流线 、海表面高度填色图 ax.streamplot(x, y, us,vs,density=1)#python 中绘制流线的函数 cb=ax.contourf(x,y,sla) #绘制填色图 # # # 设置轴、标签相关大小、名称 ax.tick_params(which='major', direction='out', length=10, width=0.99, pad=0.2, bottom=True, left=True, right=False, top=False) ax.set_title('streamline plot',pad=10,fontsize=20)#标题 ax.set_xlabel('Longtitude',fontsize=20)#x轴标签 ax.set_ylabel('Latitude',fontsize=20)#y轴标签 # # #保存图片 fig.savefig('D:\\'+'oceanas.png',format='png',dpi=150)
看一下成图：
通过选取任意一个涡旋，进行绘图，可以很明显的发现，流速的流线图与海表面高度异常数据吻合很好，可以暂定认为，流线的分布一定程度上反应了海表面高度的分布。
可能计算过程稍显粗糙，简单记录一下，感兴趣的小伙伴可以去测试一下。有问题的话欢迎评论留言交流补充，peace~~
一个努力学习python的海洋菜鸡 水平有限，欢迎指正！！！ 欢迎点赞、评论、收藏。
