我正在与 Pygrib 合作，尝试使用 NBM grib 数据获取特定纬度/经度坐标的表面温度（如果有帮助，可以在 here 获取）。
我一直试图获取一个索引值来与特定纬度和经度的代表性数据一起使用。我能够导出索引，但问题是纬度和经度似乎各有 2 个坐标。我将使用佛罗里达州迈阿密（北纬 25.7617°，西经 80.1918°）作为示例来说明这一点。如果提供了 grib 文件，则格式化为最小可重复性。
def get_grib_data(self, gribfile, shortName):
grbs = pygrib.open(gribfile)
# Temp needs level specified
if shortName == '2t':
grib_param = grbs.select(shortName=shortName, level=2)
# Convention- use short name for less than 5 chars
# Else, use name
elif len(shortName) < 5:
grib_param = grbs.select(shortName=shortName)
else:
grib_param = grbs.select(name=shortName)
data_values = grib_param[0].values
# Need varying returns depending on parameter
grbs.close()
if shortName == '2t':
return data_values, grib_param
else:
return data_values
# This function is used to find the closest lat/lon value to the entered one
def closest(self, coordinate, value):
ab_array = np.abs(coordinate - value)
smallest_difference_index = np.amin(ab_array)
ind = np.unravel_index(np.argmin(ab_array, axis=None), ab_array.shape)
return ind
def get_local_value(data, j, in_lats, in_lons, lats, lons):
lat_ind = closest(lats, in_lats[j])
lon_ind = closest(lons, in_lons[j])
print(lat_ind[0])
print(lat_ind[1])
print(lon_ind[0])
print(lon_ind[1])
if len(lat_ind) > 1 or len(lon_ind) > 1:
lat_ind = lat_ind[0]
lon_ind = lon_ind[0]
dtype = data[lat_ind][lon_ind]
else:
dtype = data[lat_ind][lon_ind]
return dtype
if __name__ == '__main__':
tfile = # Path to grib file
temps, param = grib_data.get_grib_data(tfile, '2t')
lats, lons = param[0].latlons()
j = 0
in_lats = [25.7617, 0 , 0]
in_lons = [-80.198, 0, 0]
temp = grib_data.get_local_value(temps, j, in_lats, in_lons, lats, lons)
当我执行列出的打印时，我得到以下索引：
lat_ind[0]: 182
lat_ind[1]: 1931
lon_ind[0]: 1226
lon_ind[1]: 1756
因此，如果我的纬度/经度是一维的，我只会执行 temp = data[lat[0]][lon[0]]，但在这种情况下，这将给出非代表性数据。我将如何处理纬度/经度为 2 个坐标的事实？我已经验证 lats[lat_ind[0][lat_ind1] 给出了输入纬度，经度也相同。
您无法独立于经度来评估纬度的“接近度” - 您必须评估这对坐标与输入坐标的接近程度。
纬度/经度实际上只是球面坐标。 给定两个点 (lat1,lon1) (lat2,lon2)，接近度（以大圆表示）由这两点之间的球面向量之间的角度给出（将地球近似为球体）。
您可以通过构造两个点的笛卡尔向量并取点积来计算它，即 a * b * cos(theta)，其中 theta 是您想要的。
import numpy as np
def lat_lon_cartesian(lats,lons):
lats = np.ravel(lats) #make both inputs 1-dimensional
lons = np.ravel(lons)
x = np.cos(np.radians(lons))*np.cos(np.radians(lats))
y = np.sin(np.radians(lons))*np.cos(np.radians(lats))
z = np.sin(np.radians(lats))
return np.c_[x,y,z]
def closest(in_lats,in_lons,data_lats,data_lons):
in_vecs = lat_lon_cartesian(in_lats,in_lons)
data_vecs = lat_lon_cartesian(data_lats,data_lons)
indices = []
for in_vec in in_vecs: # if input lats/lons is small list then doing a for loop is ok
# otherwise can be vectorized with some array gymnastics
dot_product = np.sum(in_vec*data_vecs,axis=1)
angles = np.arccos(dot_product) # all are unit vectors so a=b=1
indices.append(np.argmin(angles))
return indices
def get_local_value(data, in_lats, in_lons, data_lats, data_lons):
raveled_data = np.ravel(data)
raveled_lats = np.ravel(data_lats)
raveled_lons = np.ravel(data_lons)
inds = closest(in_lats,in_lons,raveled_lats,raveled_lons)
dtypes = []
closest_lat_lons = []
for ind in inds:
#if data is 2-d with same shape as the lat and lon meshgrids, then
#it should be raveled as well and indexed by the same index
dtype = raveled_data[ind]
dtypes.append(dtype)
closest_lat_lons.append((raveled_lats[ind],raveled_lons[ind]))
#can return the closes matching lat lon data in the grib if you want
return dtypes
编辑：或者使用插值法。
import numpyp as np
from scipy.interpolate import RegularGridInterpolator
#assuming a grb object from pygrib
#see https://jswhit.github.io/pygrib/api.html#example-usage
lats, lons = grb.latlons()
#source code for pygrib looks like it calls lons,lats = np.meshgrid(...)
#so the following should give the unique lat/lon sequences
lat_values = lats[:,0]
lon_values = lons[0,:]
grb_values = grb.values
#create interpolator
grb_interp = RegularGridInterpolator((lat_values,lon_values),grb_values)
#in_lats, in_lons = desired input points (1-d each)
interpolated_values = grb_interp(np.c_[in_lats,in_lons])
#the result should be the linear interpolation between the four closest lat/lon points in the data set around each of your input lat/lon points.
虚拟数据插值示例：
>>> import numpy as np
>>> lats = np.array([1,2,3])
>>> lons = np.array([4,5,6,7])
>>> lon_mesh,lat_mesh = np.meshgrid(lons,lats)
>>> lon_mesh
array([[4, 5, 6, 7],
[4, 5, 6, 7],
[4, 5, 6, 7]])
>>> lat_mesh
array([[1, 1, 1, 1],
[2, 2, 2, 2],
[3, 3, 3, 3]])
>>> z = lon_mesh + lat_mesh #some example function of lat/lon (simple sum)
>>> z
array([[ 5, 6, 7, 8],
[ 6, 7, 8, 9],
[ 7, 8, 9, 10]])
>>> from scipy.interpolate import RegularGridInterpolator
>>> lon_mesh[0,:] #should produce lons
array([4, 5, 6, 7])
>>> lat_mesh[:,0] #should produce lats
array([1, 2, 3])
>>> interpolator = RegularGridInterpolator((lats,lons),z)
>>> input_lats = np.array([1.5,2.5])
>>> input_lons = np.array([5.5,7])
>>> input_points = np.c_[input_lats,input_lons]
>>> input_points
array([[1.5, 5.5],
[2.5, 7. ]])
>>> interpolator(input_points)
array([7. , 9.5])
>>> #7 = 1.5+5.5 : correct
... #9.5 = 2.5+7 : correct
...
>>>
我知道这是一个较旧的问题，但想将其作为一个更简单的答案放在这里，它应该通过最小化 x 距离和 y 距离来满足您的需求：
coordinates = np.unravel_index((np.abs(lats - lat_to_find) + np.abs(lons - lon_to_find)).argmin(), lats.shape)
lats 和 lons 应该是纬度和经度的网格，然后 lat_to_find 和 lon_to_find 是要在数组中查找的纬度/经度对。然后你就有了索引，你可以快速索引你需要的任何变量。