[Python知识库]python 如何通过海表面高度数据计算海表地转流速、并绘制流线图

最近，学习海气相互作用时，老师布置了一个小作业。通过卫星高度计测得的海表面高度异常数据，计算其表层地转流速，并研究其与海表面高度异常的关系。
刚刚看到作业内容时，还是有点一头雾水的，不晓得怎么通过海表面高度数据推算出地转流速。挠头思考时，无意中看到物理海洋学，想起里面曾经有关系地转流部分的内容，索性拿出来看看找一下。惊喜大于意外，书里果然有推导的公式，直接贴在下面了：

其中，$\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的海洋菜鸡
												                 水平有限，欢迎指正！！！
												                    欢迎点赞、评论、收藏。