5

我正在与 Pygrib 合作，尝试使用 NBM grib 数据（如果有帮助，可在此处获得）获取特定纬度/经度坐标的表面温度。

我一直试图获取一个索引值，以用于特定纬度和经度的代表性数据。我能够得出一个索引，但问题是纬度和经度似乎各有 2 个坐标。我将使用佛罗里达州迈阿密（25.7617° N，80.1918° W）作为示例来说明这一点。如果提供了 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_ind 1 ] 给出了输入纬度和经度相同。

4

1 回答 1

6

您不能独立于经度来评估纬度的“接近度”——您必须评估这对坐标与输入坐标的接近程度。

纬度/经度实际上只是球坐标。给定两点 (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
...
>>>
于 2021-01-05T06:34:23.123 回答