将网格从中心坐标重新采样到外部（即角）坐标

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 技术可以满足您的要求，但我可能有不同的解决方案。我经常需要填充/外推到我缺少信息的网格区域。为此，我使用了一个使用 F2PY（numpy 附带）编译的 Fortran 程序，该程序将其创建为 python 模块。假设您有 Intel Fortran 编译器，您可以使用以下命令对其进行编译：f2py --verbose --fcompiler=intelem -c -m extrapolate fill.f90。您可以从 python 调用程序（完整示例参见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))

程序通过迭代方法在 RECTANGULAR 坐标中求解具有 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

示例图片（丹麦周边/瑞典南部）：