Python：从纬度/经度坐标获取网格尺寸的有效方法

Question

我需要一种非常有效的方法来提取一组纬度/经度numpy数组的网格尺寸（即2d数组的x / y索引）。过去，我是通过计算所有网格单元和经纬度坐标之间的大圆距离（使用harversine公式），然后找到最小值指数（基本上是寻找最近的点）来实现的。这是指向包含该问题的numpy网格数组的zip文件的链接。这些数组最初来自具有曲线网格的netCDF文件。

Download numpy grid arrays

这很好，但是我需要为此做三亿点，所以这种方法太慢了（将花费一年的时间）。这是我目前的做法。

import numpy as np
import math
import time

# load the numpy lat/lon grids
lat_array = np.load('lat_array.npy')
lon_array = np.load('lon_array.npy')
# get the array shape dimensions for later conversion
grid_shape = lat_array.shape

# test for lat/lon coordinate 
lat = -32
lon = 154
N = 3e8 # how many lat/lon paris I need to do

# Harversine method
def harversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance between two points 
    on the earth (specified in decimal degrees)
    """
    # convert decimal degrees to radians
    lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
    # harversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = math.sin(dlat/2.)**2. + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2.)**2.
    c = 2. * math.asin(math.sqrt(a))
    km = 6371. * c # radius of earth
    return km

# test harversine method
tic = time.perf_counter()
# create a list of distance from the point for all grid cells
dist_array = np.asarray([harversine(lon, lat, grid_lon, grid_lat) for grid_lon, grid_lat in zip(lon_array.flatten(), lat_array.flatten())])
# get the index of the minium value
min_idx = np.argmin(dist_array)
# transform the index back into grid cell dimensions
grid_dims = np.unravel_index(min_idx, grid_shape)
toc = time.perf_counter()

# report results
print('Single iteration time in seconds:', round(toc - tic, 2))
print('N iterations time in days:', round(((toc - tic)*N)/60/60/24, 2))
print('Grid coordinate:', grid_dims)

输出

Single iteration time in seconds: 0.13
N iterations time in days: 448.04
Grid coordinate: (179, 136)

我还尝试了一种隧道距离方法，该方法速度更快，但要完成所有操作仍需要15天左右。

# tunnel distance method
def tunnel_fast(latvar,lonvar,lat0,lon0):
    rad_factor = math.pi/180.0 # for trignometry, need angles in radians
    # Read latitude and longitude from file into numpy arrays
    latvals = latvar[:] * rad_factor
    lonvals = lonvar[:] * rad_factor
    ny,nx = latvals.shape
    lat0_rad = lat0 * rad_factor
    lon0_rad = lon0 * rad_factor
    # Compute numpy arrays for all values, no loops
    clat,clon = np.cos(latvals),np.cos(lonvals)
    slat,slon = np.sin(latvals),np.sin(lonvals)
    delX = np.cos(lat0_rad)*np.cos(lon0_rad) - clat*clon
    delY = np.cos(lat0_rad)*np.sin(lon0_rad) - clat*slon
    delZ = np.sin(lat0_rad) - slat;
    dist_sq = delX**2 + delY**2 + delZ**2
    minindex_1d = dist_sq.argmin()  # 1D index of minimum element
    iy_min,ix_min = np.unravel_index(minindex_1d, latvals.shape)
    return (iy_min,ix_min)

# test tunnel distance method
tic = time.perf_counter()
# create a list of distance from the point for all grid cells
grid_dims = tunnel_fast(lat_array, lon_array, lat, lon)
toc = time.perf_counter()

# report results
print('Single iteration time in seconds:', round(toc - tic, 5))
print('N iterations time in days:', round(((toc - tic)*N)/60/60/24, 2))
print('Grid coordinate:', grid_dims)

输出

Tunnel distance results
Single iteration time in seconds: 0.00445
N iterations time in days: 15.44
Grid coordinate: (179, 136)

我已经尝试了其他一些方法，例如kdtrees，但是它们也太慢了（每次迭代大约需要一秒钟）。

有人能以更快的方式获取网格单元的尺寸，或者也许我需要重新考虑整个方法吗？

Answer 1

如果您的纬度/经度网格等距分布，则无需进行任何计算，只需将当前纬度/经度除以网格大小，即可获得当前所在单元格的网格索引。如果某种方式对您不起作用，则需要解释（或提供示例）您正在谈论的网格。没有提供指向未知文件的链接，但实际上在问题中包括了网格的一部分，并解释了为什么不能通过简单除以网格大小来计算出它。