如何从纬度和经度值的中心位置查找半径内的值

Question

我正在尝试从中心纬度位置计算特定半径内包含的所有值。我使用的代码如下所示：

import numpy as np
import matplotlib.pylab as pl
import netCDF4 as nc
import haversine

f = nc.Dataset('air_temp.nc')


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

# Latitude / longitude grid
#lat = np.linspace(50,54,16)
lat = f.variables['lat'][:]
#lon = np.linspace(6,9,12)
lon = f.variables['lon'][:]
clat = 19.7
clon = 69.7
max_dist = 750      # max distance in km

# Calculate distance between center and all other lat/lon pairs
distance = haversine(lon[:,np.newaxis], lat, clon, clat) 

# Mask distance array where distance > max_dist
distance_m = np.ma.masked_greater(distance, max_dist)

# Dummy data
air = f.variables['air'][0,:,:,:]
data = np.squeeze(air)
data = np.transpose(data)
#data = np.random.random(size=[lon.size, lat.size])
data_m = np.ma.masked_where(distance  >max_dist, data)
# Test: set a value outside the max_dist circle to a large value:
#data[0,0] = 10
#avg = np.nanmean(data_m)-273

我使用过半正弦函数来求距离。现在我面临的问题是我需要距离中心点 2.5 度半径范围内的值，但我得到的都是公里。因此，如果有人可以通过说出我做错了什么或如何以正确的程序来帮助我，我们将不胜感激

Answer 1

就直线（或者更确切地说是最短弧线）距离而言，1 度始终是 111 公里（假设地球是一个完美的球体（*已编辑，而不是“正方形”）。

地球上任意两点之间最短圆弧的中心始终是地球的中心。 1 度 = 2π/360 弧度，所以距离为 R(2π/360) = 6371(2π/360) = 111.19。

更新：

你错过的不是haversine计算或者degree-km转换，而是对NetCDF的元数据格式和NumPy的meshgrid的理解。 f.variables['lat'] 给你 37 个纬度值，f.variables['lon'] 给你 144 个经度值，所以如果你想暴力搜索所有这些值，你需要使用 np.meshgrid 生成一个 37*144=5328 点的网格。

功能代码如下：

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

这是有道理的。