且构网

分享程序员开发的那些事...
且构网 - 分享程序员编程开发的那些事

如何在经度和纬度值的中心位置找到半径范围内的值

更新时间:2022-01-19 05:35:06

就直线距离(或更短的电弧)而言,1度始终是111公里(假设地球是一个完美的球体(*编辑,而不是正方形").

In terms of straight-line (or rather shortest-arc) distance, 1 degree is always 111km (assuming the earth is a perfect sphere (*edited, not "square")).

地球上任意两点之间的最短弧的中心始终是地球的中心.1度=2π/360弧度,因此距离为R(2π/360)= 6371(2π/360)= 111.19.

The center of the shortest arc between any two points on a globe is always the center of the globe. 1 degree = 2π/360 radian, so the distance is R(2π/360) = 6371(2π/360) = 111.19.

更新:

您错过的不是半正弦计算或度-公里转换,而是对 NetCDF 的元数据格式和 NumPy 的网格网格的理解.f.variables['lat'] 为您提供 37 个纬度值,f.variables['lon'] 为您提供 144 个经度值,因此如果您想进行暴力搜索所有这些,你需要使用 np.meshgrid 生成一个 37*144=5328 点的网格.

What you missed is not the haversine calculation or the degree-km conversion, it's the understanding of NetCDF's metadata format and NumPy's meshgrid. f.variables['lat'] gives you 37 latitude values and f.variables['lon'] gives you 144 longitude values, so if you want to brute force search all of them, you need to use np.meshgrid to generate a grid of 37*144=5328 points.

功能代码如下:

import numpy as np

def haversine(lon1, lat1, lon2, lat2):
    # convert decimal degrees to radians
    lon1 = np.deg2rad(lon1)
    lon2 = np.deg2rad(lon2)
    lat1 = np.deg2rad(lat1)
    lat2 = np.deg2rad(lat2)

    # haversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
    c = 2 * np.arcsin(np.sqrt(a))
    r = 6371
    return c * r

# center point
ctr_lon, ctr_lat = 69.7, 19.7

# the lon/lat grids
lon = np.arange(0, 360, 2.5)
lat = np.arange(-45, 46, 2.5)

# get coordinates of all points on the grid
grid_lon, grid_lat = np.meshgrid(lon, lat)
dists_in_km = haversine(grid_lon, grid_lat, ctr_lon, ctr_lat)
dists_in_deg = dists_in_km / 111

# find nearby points
thr = 2.5
for i in range(grid_lon.shape[0]):
    for j in range(grid_lon.shape[1]):
        this_lon = grid_lon[i, j]
        this_lat = grid_lat[i, j]
        dist = dists_in_deg[i, j]
        if dist <= thr:
            print('lon=%.1f  lat=%.1f dist=%.2fdeg' % (this_lon, this_lat, dist))

输出:

lon=70.0  lat=17.5 dist=2.22deg
lon=67.5  lat=20.0 dist=2.09deg
lon=70.0  lat=20.0 dist=0.41deg

这很有意义.