我有一個(gè) netCDF 文件，其中包含多個(gè)氣旋位置（緯度、經(jīng)度）和特定時(shí)間南半球的氣溫。我想要的是從每個(gè)氣旋位置的中心提取 10 測(cè)地線度數(shù)（~1110 公里）半徑范圍內(nèi)的溫度值。這個(gè)想法是確定與每個(gè)氣旋相關(guān)的溫度值（假設(shè)距氣旋中心的最大徑向距離為 10o），并繪制一個(gè)僅包含這些溫度值的全局 contourf 地圖。我在這里搜索了很多，但我只找到了適用于與一個(gè)特定經(jīng)緯度中心的距離的代碼（比如這個(gè)：如何從緯度和經(jīng)度值的中心位置找到半徑內(nèi)的值）。我一直在研究如何同時(shí)為多個(gè)中心應(yīng)用 Haversine 公式。import xarray as xrimport numpy as npimport matplotlib.pyplot as pltd = xr.open_dataset('cyc_temp.nc')lat = d['lat']lon = d['lon']cyc_pos = d['id'][:,:]temp = d['temp'][:,:]# Haversine formuladef 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如果有人可以幫助我，我會(huì)很感激。
1 回答
紅糖糍粑
TA貢獻(xiàn)1815條經(jīng)驗(yàn) 獲得超6個(gè)贊
這是一個(gè)有趣的問題；xarray 的自動(dòng)廣播使它非常干凈。
我不確定氣旋位置數(shù)組的結(jié)構(gòu)，但我假設(shè)它的結(jié)構(gòu)如下（或者至少可以操縱成這種形式）：
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"])
換句話說，每一行代表每個(gè)旋風(fēng)的經(jīng)度和緯度值。
以cyc_pos這種方式定義，使用該函數(shù)獲取經(jīng)緯度網(wǎng)格中的每個(gè)點(diǎn)到每個(gè)氣旋中心的距離haversine非常簡(jiǎn)單，從那里獲取所需的掩碼僅需多一行。
distances = haversine(cyc_pos.sel(coord="lon"), cyc_pos.sel(coord="lat"), lon, lat)
如果您想要針對(duì)特定風(fēng)暴的面罩，您可以使用：
storm_id = 0
mask = (distances <= 1110.0).isel(cyclone=storm_id)
或者，如果您想要一個(gè)可以應(yīng)對(duì)所有風(fēng)暴的面具，您可以使用：
mask = (distances <= 1110.0).any("cyclone")
添加回答
舉報(bào)
0/150
提交
取消