從中心座標到外部（即角點）座標的重新採樣網格

Question

是否有從網格中心位置（紅色點）外推網格角位置（藍色點）的現成方法？

我正在使用的網格不是矩形的，所以規則的雙線性插值似乎並不是最好的方式;不過，這只是我的情節我的數據使用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]]) 

Answer 1

我不知道這樣做你問任何強大的matplotlib技術的，但我可能有不同的解決方案。我經常需要填補/外推到我缺少信息的網格區域。爲此，我使用一個Fortran程序，我使用F2PY編譯（隨numpy一起提供），將其創建爲一個python模塊。假設您擁有英特爾Fortran編譯器，請使用以下命令編譯它：f2py --verbose --fcompiler=intelem -c -m extrapolate fill.f90。可以調用從蟒程序用（參見here爲完整的例子）：

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帳戶here上使用。

Answer 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 

實施例的圖像（大約丹麥/瑞典南部）：