识别最近的网格点

Posted 2023-02-16

【中文标题】识别最近的网格点

【英文标题】：Identifying the nearest grid point

【发布时间】：2015-06-16 17:03:40

【问题描述】：

我有三个数组

lat=[15,15.25,15.75,16,....30]
long=[91,91.25,91.75,92....102]

data=
array([[  0. ,   0. ,   0. , ...,   0. ,   0. ,   0. ],
       [  0. ,   0. ,   0. , ...,   0. ,   0. ,   0. ],
       [  0. ,   0. ,   0. , ...,   0. ,   0. ,   0. ],
       ..., 
       [-99.9, -99.9, -99.9, ...,   0. ,   0. ,   0. ],
       [-99.9, -99.9, -99.9, ...,   0. ,   0. ,   0. ],
       [-99.9, -99.9, -99.9, ...,   0. ,   0. ,   0. ]])

它是 [44 cols and 60 rows] 与 long x lat 相同

如果我输入任何点 (16.3,101.6)，我需要找出最近的网格并从第三个数组中从该网格中提取数据。我怎样才能在 python 中使用 numpy 呢？这里我举一个例子，但在实际问题中，我有几点。

这个功能我试过了，

def getclosest_ij(lats,lons,latpt,lonpt):
    dis_sq1=(lats-latpt)
    dis_sq2=(lons-lonpt)
    minidex_lat=dis_sq1.argmin()
    minidex_lon=dis_sq2.argmin()
    return minidex_lon,minidex_lat

【问题讨论】：

请向我们展示您解决此问题的尝试

我试图计算每个 lat 和 long 值的绝对距离，并返回该索引以提取该数据，但我得到了错误的值。我在我的问题中添加了我的功能。

只需为差值np.abs(lats-latpt) 等添加绝对值，您的实现应该可以正常工作。

【参考方案1】：

您正在寻找的算法是规则网格上的最近邻插值。例如，您可以使用，

from scipy.interpolate import RegularGridInterpolator

itp = RegularGridInterpolator( (lat, lon), data, method='nearest') 
res = itp(some_new_point)

另外，如果您设置method='linear'，此函数还可以执行更精确的线性插值。

【讨论】：

用 numpy 自己能找出来吗？

【参考方案2】：

这应该可行：

import numpy as np

lat = np.array(lat)
long = np.array(long)

def get_data(lat_input, long_input):

    lat_index  = np.nanargmin((lat-lat_input)**2)
    long_index = np.nanargmin((long-long_input)**2)
    return data[long_index][lat_index]

您需要 numpy 数组格式的 lat 和 long 数据才能与 nanargmin 函数一起使用。如果您确定数据数组中不会有任何 nan 值，您也可以使用 argmin 而不是 nanargmin。

我将差平方而不是取绝对值，因为前者稍微快一些，并且无论如何都会找到相同的索引。

编辑：正如 rth 在 cmets 中对 tom 的回答所指出的那样，argmin 明显快于nanargmin。如果您的数据可能有nan 值，您可以简单地事先更正它们，然后安全地使用argmin。此外，正如汤姆所说，如果您的数据已排序，searchsorted 确实是更好的选择。

【参考方案3】：

您可以为此使用numpy.searchsorted：

import numpy as np

lat=np.linspace(15,30,61)
long=np.linspace(91,102,45)

def find_index(x,y):
    xi=np.searchsorted(lat,x)
    yi=np.searchsorted(long,y)
    return xi,yi

thisLat, thisLong = find_index(16.3,101.6)
print thisLat, thisLong
>>> 6, 43

# You can then access the `data` array like so:
print data[thisLat,thisLong]

注意：这将找到低于兴趣点的lat 和long 数组的索引。如果你需要最近的点，你可以比较lat[thisLat]和lat[thisLat+1]的值

【讨论】：

+1 非常好。使用上面定义的变量，np.searchsorted(lat,x) 比我计算机上的等效调用 np.nanargmin((lat-x)**2) 快 16 倍。

好吧，实际上这不是 NaN 搜索开销，在这个例子中它只比 np.argmin(..) 快 x2。数组大小 (n) 的缩放仍然要好得多：x5 表示 n=1000，x20 表示 n=10000，等等。