如何查找距离纬度和经度值中心位置的半径范围内的值

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

我使用haversine函数来查找距离。现在我面临的问题是我需要距离中心点 2.5 度半径范围内的值，但我都以公里为单位。因此，如果有人能帮助我说我做错了什么或如何以正确的程序进行，我们将不胜感激

Answer 1

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

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

更新：

你错过的不是haversine计算或度公里转换，而是对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

这是有道理的。