从 NetCDF 中的多个经纬度中心查找半径内的值
Find values within a radius from multiples lat lon centers in a NetCDF
我有一个 netCDF 文件,其中包含特定时间南半球的多个气旋位置(纬度、经度)和气温。
我想要的是从每个气旋位置的中心提取10个测地度(~1110公里)半径范围内的温度值。这个想法是确定与每个气旋相关的温度值(假设距气旋中心的最大径向距离为 10º)并绘制一个仅包含这些温度值的全局 contourf 地图。
我在这里搜索了很多,但我只找到了适用于距一个特定经纬度中心的距离的代码(例如这个:)。
我卡在如何一次对多个中心应用 Haversine 公式。
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
如果有人能帮助我,我会很感激。
这是一个有趣的问题; 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.(coord="lat"), lon, lat)
如果您想要针对特定风暴的面具,您可以使用:
storm_id = 0
mask = (distances <= 1110.0).isel(cyclone=storm_id)
或者如果你想要一个可以抵御所有风暴的面具,你可以使用:
mask = (distances <= 1110.0).any("cyclone")
我有一个 netCDF 文件,其中包含特定时间南半球的多个气旋位置(纬度、经度)和气温。
我想要的是从每个气旋位置的中心提取10个测地度(~1110公里)半径范围内的温度值。这个想法是确定与每个气旋相关的温度值(假设距气旋中心的最大径向距离为 10º)并绘制一个仅包含这些温度值的全局 contourf 地图。
我在这里搜索了很多,但我只找到了适用于距一个特定经纬度中心的距离的代码(例如这个:
我卡在如何一次对多个中心应用 Haversine 公式。
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
如果有人能帮助我,我会很感激。
这是一个有趣的问题; 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.(coord="lat"), lon, lat)
如果您想要针对特定风暴的面具,您可以使用:
storm_id = 0
mask = (distances <= 1110.0).isel(cyclone=storm_id)
或者如果你想要一个可以抵御所有风暴的面具,你可以使用:
mask = (distances <= 1110.0).any("cyclone")