|
##数据来源：NCEP再分析数据
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
#import scipy.stats as ss
import netCDF4 as nc
import math
file1=r'C:\Users\59799\Desktop\uwnd.mon.mean.nc'
file2=r'C:\Users\59799\Desktop\vwnd.mon.mean.nc'
file3=r'C:\Users\59799\Desktop\air.mon.mean.nc'
uwnd=nc.Dataset(file1);vwnd=nc.Dataset(file2)
air=nc.Dataset(file3)
###通用量读取
lat=np.array(uwnd.variables['lat'])
lon=np.array(uwnd.variables['lon'])
level=np.array(uwnd.variables['level'])
###u,v,airtemp获取
u=np.array(uwnd.variables['uwnd'])
v=np.array(vwnd.variables['vwnd'])
airtemp=np.array(air.variables['air'])
##水平温度平流计算公式为 -(udt/dx+vdt/dy)
##对于离散数据,无法做微分会偏微分运算,采取计算方法中
##提到的差分方法，这里利用前差
##对于level和time的选取,因为本code仅仅做示例
##我们选取level为850hPa,time为1948-01-01
##即level[2];time[0]
##还有一个重要问题:dx为某一纬度,通过2.5经度的距离
##即不同纬度dx不等(考虑球坐标);但是对于dy从赤道到极点是固定的
#--------------------------------------#
#首先计算dy
R=6371*1000#地球半径6371km
theta=math.radians(2.5)##角度转弧度,数据经纬度间隔为2.5°
dy=R*theta
#计算dx
dx=[]
###r为90N到0的小圆半径
for i in range(len(lat[0:73])):
print(lat[i])
arg=math.radians(lat[i])
r=R*math.cos(arg)
#print(r)
dx.append(r*math.radians(2.5))
dx=np.array(dx)##dx为90N到0每一纬圈经过2.5°经度的距离差
dx[0]=0;dx[72]=0
##先求dt/dx
###因为是圆形区域,所以前差并不会导致横向数据量减少
dtdx=np.empty((73,144))
for i in range(len(lat)):
for j in range(len(lon)):
if (i==0) or (i==72):
dtdx[i,j]=0
else:
if j==143:
dtdx[i,0]=(airtemp[1,2,i,0]-airtemp[1,2,i,j])/dx[i]
else:
dtdx[i,j+1]=(airtemp[1,2,i,j+1]-airtemp[1,2,i,j])/dx[i]
###再求dt/dy
dtdy=np.zeros((73,144))
for j in range(len(lon)):
for i in range(len(lat)):
if i==72:
dtdy[i,j]=0
else:
dtdy[i+1,j]=(airtemp[1,2,i+1,j]-airtemp[1,2,i,j])/dy
###得出温度平流 -(udt/dx+vdt/dy)
######!!!!!!!唯独87.5往上无数据
######!!!!!!!唯独87.5往上无数据
######!!!!!!!唯独87.5往上无数据
dTdt=np.zeros((73,144))
for i in range(len(lat)):
for j in range(len(lon)):
dTdt[i,j]=-(u[1,2,i,j]*dtdx[i,j]+v[1,2,i,j]*dtdy[i,j])
#绘图
[X,Y]=np.meshgrid(lon[:],lat[1:72])
LON_0=lon.mean()
LAT_0=lat.mean()
map=Basemap(lon_0=LON_0,lat_0=LAT_0,llcrnrlon=lon[0],llcrnrlat=lat[71],urcrnrlon=lon[143],urcrnrlat=lat[1])
map.drawcoastlines()
c=plt.contourf(X,Y,dTdt[1:72,:],cmap='coolwarm')
plt.colorbar(c)
plt.title(' 1948 JAN. Temperature Advection (87.5N-87.5S) ') | 
转播
