how to find values within a radius from a central position of latitude and longitude value

Question

I am trying to calculate all the values contained within a particular radius from a central lat lon position.The code which I am using is as given:

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

I have used the haversine function for finding the distance. Now what I am facing the problem is I need values within a radius of 2.5 degrees from the center point, but I am getting all in kilometers. So if anyone can help me by saying what I am doing wrong or how to it in the right procedure it will be highly appreciated

Answer 1

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")).

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.

Update:

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.

Functional code below:

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))

Output:

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

which makes sense.