从NetCDF中的多个横向中心查找半径内的值

Question

我有一个 netCDF 文件，其中包含特定时间南半球的多个气旋位置（纬度、纬度）和气温。

我想要的是提取温度值​​距离每个气旋位置中心10个测地线度（~1110公里）的半径内。想法是确定温度值​​与每个气旋相关（假设与气旋中心的最大径向距离为10º），并绘制一张仅包含这些温度值的全球等值线图。

我在这里搜索了很多，但我只找到了适用于距离一个特定的纬度中心的距离的代码（比如这个：如何从纬度和经度值的中心位置找到半径内的值）。

我被困在如何将哈弗辛公式同时应用于多个中心。

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

如果有人能帮助我，我愿意。

Answer 1

这是一个有趣的问题；xarray的自动广播让这一切变得非常干净。

我不确定旋风位置阵列的结构，但我会假设它的结构如下（或者至少可以被操纵为这种形式）：

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

换句话说，每一行代表每个气旋的经度和纬度值。

以这种方式定义cyc_pos，使用 haversine 函数获取经纬度网格中每个点到每个气旋中心的距离相当简单，从那里获得所需的掩码只是多一条线。

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

如果你想要一个特定Storm的面具，你可以使用：

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

或者如果你想要一个遮蔽所有Storm的面具：

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