Find values within a radius from multiples lat lon centers in a NetCDF

I have a netCDF file with multiple cyclone positions (lat, lon) and air temperature across the southern hemisphere in a particular time.

What I want is to extract the temperature values ​​that are within a radius of 10 geodesic degrees (~1110 km) from the center of each cyclone position. The idea is to identify the temperature values ​​associated with each cyclone (assuming a maximum radial distance of 10º from the cyclone center) and plot one global contourf map with only those temperature values.

I searched a lot here, but I only found codes that applies to distances from just one specific lat lon center (like this one: how to find values within a radius from a central position of latitude and longitude value).

Im stuck in how to apply the Haversine formula for multiple centers at once.

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt

d = xr.open_dataset('cyc_temp.nc')
lat = d['lat']
lon = d['lon']
cyc_pos = d['id'][:,:]
temp = d['temp'][:,:]

# Haversine formula

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

If anyone can help me, I would aprecciate.

This is a fun problem; xarray's automatic broadcasting makes this pretty clean.

I'm not sure how the array of cyclone positions is structured, but I'm going to assume it is structured like the following (or at least could be manipulated to be in this form):

centers = np.array([[12.0, -62.0], [40.0, -80.0], [67.0, -55.0]])
cyc_pos = xr.DataArray(centers, coords={"coord": ["lon", "lat"]}, dims=["cyclone", "coord"])

In other words, each row represents the longitude and latitude values of each cyclone.

With cyc_pos defined in that way, obtaining the distances of each point in the latitude-longitude grid to each cyclone center using the haversine function is fairly straightforward, and from there obtaining the desired mask is only one more line.

distances = haversine(cyc_pos.sel(coord="lon"), cyc_pos.sel(coord="lat"), lon, lat)

If you want a mask for a specific storm you could use:

storm_id = 0
mask = (distances <= 1110.0).isel(cyclone=storm_id)

or if you wanted a mask for all the storms you could use:

mask = (distances <= 1110.0).any("cyclone")