这是我采用的不太优雅的方法,使用pyproj
首先计算点之间的距离和方位角(使用pyproj.Geod.inv
),然后通过必要的距离插值/外推该角度(使用pyproj.Geod.fwd
)到psi位置。
代码:
def calc_psi_coords(lons, lats):
''' Calcuate psi points from centered grid points'''
import numpy as np
import pyproj
g = pyproj.Geod(ellps='WGS84')
ydim, xdim = lons.shape
lons_psi = np.zeros((ydim+1, xdim+1))
lats_psi = np.zeros((ydim+1, xdim+1))
for j in range(ydim-1):
for i in range(xdim-1):
lon1 = lons[j,i]
lat1 = lats[j,i]
lon2 = lons[j+1,i+1]
lat2 = lats[j+1,i+1]
fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2)
lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*0.5)
lons_psi[j+1,i+1] = lon_psi
lats_psi[j+1,i+1] = lat_psi
for j in range(ydim):
lon1 = lons_psi[j+1,2]
lat1 = lats_psi[j+1,2]
lon2 = lons_psi[j+1,1]
lat2 = lats_psi[j+1,1]
fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2)
lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.)
lons_psi[j+1,0] = lon_psi
lats_psi[j+1,0] = lat_psi
lon1 = lons_psi[j+1,-3]
lat1 = lats_psi[j+1,-3]
lon2 = lons_psi[j+1,-2]
lat2 = lats_psi[j+1,-2]
fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2)
lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.)
lons_psi[j+1,-1] = lon_psi
lats_psi[j+1,-1] = lat_psi
for i in range(xdim):
lon1 = lons_psi[2,i+1]
lat1 = lats_psi[2,i+1]
lon2 = lons_psi[1,i+1]
lat2 = lats_psi[1,i+1]
fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2)
lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.)
lons_psi[0,i+1] = lon_psi
lats_psi[0,i+1] = lat_psi
lon1 = lons_psi[-3,i+1]
lat1 = lats_psi[-3,i+1]
lon2 = lons_psi[-2,i+1]
lat2 = lats_psi[-2,i+1]
fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2)
lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.)
lons_psi[-1,i+1] = lon_psi
lats_psi[-1,i+1] = lat_psi
lon1 = lons_psi[2,2]
lat1 = lats_psi[2,2]
lon2 = lons_psi[1,1]
lat2 = lats_psi[1,1]
fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2)
lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.)
lons_psi[0,0] = lon_psi
lats_psi[0,0] = lat_psi
lon1 = lons_psi[2,-3]
lat1 = lats_psi[2,-3]
lon2 = lons_psi[1,-2]
lat2 = lats_psi[1,-2]
fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2)
lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.)
lons_psi[0,-1] = lon_psi
lats_psi[0,-1] = lat_psi
lon1 = lons_psi[-3,2]
lat1 = lats_psi[-3,2]
lon2 = lons_psi[-2,1]
lat2 = lats_psi[-2,1]
fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2)
lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.)
lons_psi[-1,0] = lon_psi
lats_psi[-1,0] = lat_psi
lon1 = lons_psi[-3,-3]
lat1 = lats_psi[-3,-3]
lon2 = lons_psi[-2,-2]
lat2 = lats_psi[-2,-2]
fwd_az, bck_az, dist = g.inv(lon1,lat1,lon2,lat2)
lon_psi, lat_psi, bck_az = g.fwd(lon1,lat1,fwd_az,dist*2.)
lons_psi[-1,-1] = lon_psi
lats_psi[-1,-1] = lat_psi
return lons_psi, lats_psi
示例图片(位于丹麦/瑞典南部附近):
![输入图像描述](https://istack.dev59.com/MtZ9R.webp)
F2PY
。这个工具在其他方面也可能很有用。我正在绘制NORWECOM.e2e
模型的模块输出。我会等待其他输入,看看是否还有其他可用的东西。谢谢! - ryanjdillon