从中心坐标重采样到外部（即角落）坐标
有没有简单的方法可以从网格中心的位置（红点）推算出网格角落的位置（蓝点）呢？
我现在用的网格不是矩形的，所以普通的双线性插值方法可能不太合适。不过，这只是为了让我能用 pyplot.pcolormesh() 来绘制我的数据，所以也许这并不是特别重要。
示例网格数据：
import numpy as np
lons = np.array([[ 109.93299681, 109.08091365, 108.18301276, 107.23602539],
[ 108.47911382, 107.60397996, 106.68325946, 105.71386119],
[ 107.06790187, 106.17259769, 105.23214707, 104.2436463 ],
[ 105.69908292, 104.78633156, 103.82905363, 102.82453812]])
lats = np.array([[ 83.6484245 , 83.81088466, 83.97177823, 84.13098916],
[ 83.55459198, 83.71460466, 83.87294803, 84.02950188],
[ 83.4569054 , 83.61444708, 83.77022192, 83.92410637],
[ 83.35554612, 83.51060313, 83.6638013 , 83.81501464]])
2 个回答
2
这是我用 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
# Create Geod object with WGS84 ellipsoid
g = pyproj.Geod(ellps='WGS84')
# Get grid field dimensions
ydim, xdim = lons.shape
# Create empty coord arrays
lons_psi = np.zeros((ydim+1, xdim+1))
lats_psi = np.zeros((ydim+1, xdim+1))
# Calculate internal points
#--------------------------
for j in range(ydim-1):
for i in range(xdim-1):
lon1 = lons[j,i] # top left point
lat1 = lats[j,i]
lon2 = lons[j+1,i+1] # bottom right point
lat2 = lats[j+1,i+1]
# Calc distance between points, find position at half of dist
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)
# Assign to psi interior positions
lons_psi[j+1,i+1] = lon_psi
lats_psi[j+1,i+1] = lat_psi
# Caclulate external points (not corners)
#----------------------------------------
for j in range(ydim):
# Left external points
#~~~~~~~~~~~~~~~~~~~~~
lon1 = lons_psi[j+1,2] # left inside point
lat1 = lats_psi[j+1,2]
lon2 = lons_psi[j+1,1] # left outside point
lat2 = lats_psi[j+1,1]
# Calc dist between points, find position at dist*2 from pos1
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
# Right External points
#~~~~~~~~~~~~~~~~~~~~~~
lon1 = lons_psi[j+1,-3] # right inside point
lat1 = lats_psi[j+1,-3]
lon2 = lons_psi[j+1,-2] # right outside point
lat2 = lats_psi[j+1,-2]
# Calc dist between points, find position at dist*2 from pos1
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):
# Top external points
#~~~~~~~~~~~~~~~~~~~~
lon1 = lons_psi[2,i+1] # top inside point
lat1 = lats_psi[2,i+1]
lon2 = lons_psi[1,i+1] # top outside point
lat2 = lats_psi[1,i+1]
# Calc dist between points, find position at dist*2 from pos1
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
# Bottom external points
#~~~~~~~~~~~~~~~~~~~~~~~
lon1 = lons_psi[-3,i+1] # bottom inside point
lat1 = lats_psi[-3,i+1]
lon2 = lons_psi[-2,i+1] # bottom outside point
lat2 = lats_psi[-2,i+1]
# Calc dist between points, find position at dist*2 from pos1
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
# Calculate corners:
#-------------------
# top left corner
#~~~~~~~~~~~~~~~~
lon1 = lons_psi[2,2] # bottom right point
lat1 = lats_psi[2,2]
lon2 = lons_psi[1,1] # top left point
lat2 = lats_psi[1,1]
# Calc dist between points, find position at dist*2 from pos1
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
# top right corner
#~~~~~~~~~~~~~~~~~
lon1 = lons_psi[2,-3] # bottom left point
lat1 = lats_psi[2,-3]
lon2 = lons_psi[1,-2] # top right point
lat2 = lats_psi[1,-2]
# Calc dist between points, find position at dist*2 from pos1
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
# bottom left corner
#~~~~~~~~~~~~~~~~~~~
lon1 = lons_psi[-3,2] # top right point
lat1 = lats_psi[-3,2]
lon2 = lons_psi[-2,1] # bottom left point
lat2 = lats_psi[-2,1]
# Calc dist between points, find position at dist*2 from pos1
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
# bottom right corner
#~~~~~~~~~~~~~~~~~~~~
lon1 = lons_psi[-3,-3] # top left point
lat1 = lats_psi[-3,-3]
lon2 = lons_psi[-2,-2] # bottom right point
lat2 = lats_psi[-2,-2]
# Calc dist between points, find position at dist*2 from pos1
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
示例图片（丹麦/南瑞典附近）:
4
我不知道有没有什么特别好的方法可以用matplotlib来解决你提到的问题，但我有一个不同的解决方案。通常我需要填补或推测一些网格区域的信息，这些地方可能缺少数据。为此，我使用一个Fortran程序，并通过F2PY（这个工具和numpy一起提供）把它编译成一个python模块。如果你有Intel Fortran编译器，可以用以下命令来编译它：f2py --verbose --fcompiler=intelem -c -m extrapolate fill.f90。然后你可以在python中调用这个程序（完整示例可以查看这里）：
import extrapolate as ex
undef=2.0e+35
tx=0.9*undef
critx=0.01
cor=1.6
mxs=100
field = Zg
field=np.where(abs(field) > 50 ,undef,field)
field=ex.extrapolate.fill(int(1),int(grdROMS.xi_rho),
int(1),int(grdROMS.eta_rho),
float(tx), float(critx), float(cor), float(mxs),
np.asarray(field, order='Fortran'),
int(grdROMS.xi_rho),
int(grdROMS.eta_rho))
这个程序通过迭代的方法解决了拉普拉斯方程，并使用了Neumann边界条件（dA/dn = 0），在矩形坐标系中填补那些值为“undef”的网格点，给出合理的数值。对我来说这个方法效果很好，也许你也会觉得有用。这个程序可以在我的github账号上找到，链接在这里。