<jupyter_start><jupyter_text>7. 進階運算與統計方法 網格內插 (Interpolation)不同的資料通常有不同的時間和網格解析度，根據需求有時需要將資料`A`資料內插到資料`B`的網格上。`xarray.DataArray` 和`xarray.Dataset` 有內建網格內插的方法`xr.interp`，請見以下範例。**Example 1:** GPCP降雨資料網格解析度為2.5度，OLR網格資料解析度為1度，欲將OLR資料內插到GPCP資料網格上。<jupyter_code>import xarray as xr
lats = -20
latn =  30
lon1 =  79 
lon2 = 161

with xr.open_dataset('data/gpcp_precip_1979-2019.pentad.nc') as pcp_ds: 
    pcp = pcp_ds.sel(lat=slice(latn,lats), lon=slice(lon1,lon2)).data
with xr.open_dataset('data/olr.nc') as olr_ds:
    olr = olr_ds.sel(lat=slice(lats,latn), lon=slice(lon1,lon2)).olr

olr_rmp = olr.interp(lon=pcp.lon, lat=pcp.lat)   # 利用interp的method進行網格內插。
olr_rmp<jupyter_output><empty_output><jupyter_text>改變網格解析度時間和空間的網格解析度可以用採用 `xarray.DataArray.coarsen` 更改。**Example 2:** 將每日OLR轉換成侯(pentad)平均資料，意即將時間解析度調整為5天一次。<jupyter_code>olr_noleap = olr.sel(time=~((olr.time.dt.month == 2) & (olr.time.dt.day == 29)))
olr_ptd = (olr_noleap.coarsen(time=5,
                              coord_func={"time": "min"})  
                     .mean())
olr_ptd<jupyter_output><empty_output><jupyter_text>也就是在時間的維度上，5天為一個單位平均。 `coord_func` 預設為 `mean`，也就是重新取樣後的座標 (時間) 軸會選擇5天平均後的時間軸 (如1998-01-03, 1998-01-08, 1998-01-13, ...)，其他選項可選 `min` 設定為最小值 (一侯的開始日) 或 `max` 最大值 (一侯的結束日)。**Note:** 雖然`xarray.DataArray.resample`看似也有同樣功能，但根據API reference，> The resampled dimension must be a datetime-like coordinate.根據以上範例，就是 `olr_ptd = olr_noleap.resample(time='5D')` ，然而`datetime` 物件無法選擇日曆格式，也就是在計算五天平均時，不論原始資料 `olr_noleap` 有沒有2/29，這個方法都會納入2/29 (延伸閱讀：[StackOverflow 的問答](https://stackoverflow.com/questions/64476727/can-we-assign-the-calendar-of-pandas-datetimeindex-as-noleap-or-can-datetime))。 ```{admonition} Exercise:class: seealso請利用coarsen改變空間的解析度 (i.e. 將1˚解析度的OLR資料regrid成2˚解析度)。```  滑動平均 (running mean)可利用 `xarray.DataArray.rolling`，連續選擇3個單位，再進行平均。<jupyter_code>olr_3p_runave = (olr_ptd.rolling(time=3,        # 選擇幾個時間單位
                                 center=False)
                        .mean())
olr_3p_runave<jupyter_output><empty_output><jupyter_text>由於頭尾會有缺失值NaN，可加上`dropna()`移除之。  相關係數地圖 (correlation map)給定兩個DataArray，我們就可以計算兩者之間沿著特定座標軸 (e.g. 沿著時間軸的) 的相關係數。**Example 3:** 計算NDJ [Oceanic Niño Index (ONI)](https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php) 和十二月降雨的相關係數地圖。Step 1: 得到GPCP十二月平均資料 (先所有選12月資料，再用`groupby`)<jupyter_code>pcp_dec = (pcp.sel(time=pcp.time.dt.month.isin([12]))
              .groupby('time.year')
              .mean("time"))
pcp_dec<jupyter_output><empty_output><jupyter_text>Step 2: 製造NDJ ONI index矩陣。<jupyter_code>oni_ndj = xr.DataArray( 
                       data=[0.6, 0.0, -0.1, 2.2, -0.9, -1.1, -0.4, 1.2, 1.1, -1.8, -0.1, 0.4, 
                             1.5, -0.1, 0.1, 1.1, -1.0, -0.5, 2.4, -1.6, -1.7, -0.7, 
                             -0.3, 1.1, 0.4, 0.7, -0.8, 0.9, -1.6, -0.7, 1.6, -1.6, 
                             -1.0, -0.2, -0.3, 0.7, 2.6, -0.6, -1.0, 0.8, 0.5],
                       dims='year',
                       coords=dict(year=pcp_dec.year)
                       )
oni_ndj<jupyter_output><empty_output><jupyter_text>Step 3: 用`xr.corr`計算相關係數。<jupyter_code>corr = xr.corr(pcp_dec,oni_ndj,dim='year')
corr<jupyter_output><empty_output><jupyter_text>Step 4: 將計算好的結果畫出來看看。<jupyter_code>import numpy as np
from cartopy import crs as ccrs   
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib as mpl
from matplotlib import pyplot as plt
mpl.rcParams['figure.dpi'] = 100

plt.figure()                  
proj = ccrs.PlateCarree()     
ax = plt.axes(projection=proj)
ax.set_extent([lon1,lon2,lats,latn],crs=proj)
ax.set_xticks(np.arange(80,180,20), crs=proj)
ax.set_yticks(np.arange(-20,40,10), crs=proj)
lon_formatter = LONGITUDE_FORMATTER
lat_formatter = LATITUDE_FORMATTER   
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter) 
ax.coastlines()

clevs = np.arange(170,350,20)

# 繪圖
plt.title("Correlation map (Dec precip. and DJF ONI)", loc='left')    # 設定圖片標題，並且置於圖的左側。
olrPlot = (corr.plot.contourf("lon", "lat",  
                              ax=ax,       
                              levels=np.arange(-1,1.1,0.1), 
                              cmap='bwr',  
                              add_colorbar=True, 
                              extend='neither', 
                              cbar_kwargs={'orientation': 'horizontal', 'aspect': 30, 'label': ' '}) #設定color bar
                              )
                              
ax.set_ylabel(' ')             # 設定坐標軸名稱。
ax.set_xlabel(' ')
plt.show()<jupyter_output><empty_output><jupyter_text>`where`條件控制`where`是一個很好用的方法，可以針對資料設定條件進行篩選和過濾。我們直接看以下例子：**Example 4:** 繪製2021年12月平均之Mean Sea Level Pressure，但海拔3000公尺以上的區域不畫。Step 1: 讀MSLP、地形資料。<jupyter_code>lats=0
latn=60

topo_ds = xr.open_dataset('data/etopo5.nc')
mslp_ds = xr.open_dataset('data/mslp.2021.nc')
topo = topo_ds.sel(Y=slice(lats,latn),
                   X=slice(lon1,lon2)).bath
mslp = mslp_ds.sel(time=slice('2021-12-01','2021-12-31'),
                   lat=slice(latn,lats),
                   lon=slice(lon1,lon2)).mslp
mslp = mslp/100.<jupyter_output><empty_output><jupyter_text>`where`的原理是搜尋兩者相對應的網格是否滿足給定的控制條件 （以這個範例來說是搜尋topo >= 3000的網格點），但兩者網格解析度不同時，沒有辦法比較。因此必須先進行網格內插，讓`topo`和`mslp`的解析度相同。<jupyter_code>topo_rmp = topo.interp(X=mslp.lon, Y=mslp.lat)<jupyter_output><empty_output><jupyter_text>接著用`where`設定條件：<jupyter_code>mslp_mask = mslp.where(topo_rmp<=3000)
mslp_mask<jupyter_output><empty_output><jupyter_text>我們設定的條件就是`topo<=3000`，不滿足的會直接被設定為缺失值。也可以在`where`的引數`other`中，設定當不滿足條件時要代入的值。<jupyter_code>plt.figure()                  
proj = ccrs.PlateCarree()     
ax = plt.axes(projection=proj)
ax.set_extent([lon1,lon2,lats,latn],crs=proj)
ax.set_xticks(np.arange(80,180,20), crs=proj)
ax.set_yticks(np.arange(0,70,10), crs=proj)
lon_formatter = LONGITUDE_FORMATTER
lat_formatter = LATITUDE_FORMATTER   
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter) 
ax.coastlines()

# 繪圖
plt.title("MSLP, Dec. 2021", loc='left')    # 設定圖片標題，並且置於圖的左側。
olrPlot = (mslp_mask.mean(axis=0)
                    .plot.contourf("lon", "lat",  
                                   ax=ax,       
                                   levels=np.arange(986,1048,8), 
                                   cmap='coolwarm',  
                                   add_colorbar=True, 
                                   extend='both', 
                                   cbar_kwargs={'orientation': 'horizontal', 'aspect': 30, 'label': ' '}) #設定color bar
                                   )
                              
ax.set_ylabel(' ')             # 設定坐標軸名稱。
ax.set_xlabel(' ')
plt.show()<jupyter_output><empty_output><jupyter_text>計算相對渦度 (vorticity)、散度 (divergence)以上的計算都只是一些基本的統計，利用xarray的函數就可以辦到。但如果想要計算一些氣象的函數應該怎麼辦呢？我們可以使用MetPy，目前MetPy對xarray有完整的支援，可以直接將DataArray讀入MetPy的函數中。接下來我們以計算相對渦度為例，渦度和散度是氣象上常用來衡量風的旋轉量和輻合散的方法，[metpy.calc.vorticity](https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.vorticity.html) 可以用來計算渦度。**Example 5:** 計算2017年12月平均之850-hPa 相對渦度。Step 1: 讀風場資料、選擇時空範圍。<jupyter_code>latn = 30 
lats = -20

uds = xr.open_dataset('data/ncep_r2_uv850/u850.2017.nc')
vds = xr.open_dataset('data/ncep_r2_uv850/v850.2017.nc')
u = uds.sel(time=slice('2017-12-01','2017-12-31'), 
            level=850,
            lat=slice(latn,lats),
            lon=slice(lon1,lon2)).uwnd
v = vds.sel(time=slice('2017-12-01','2017-12-31'), 
            level=850,
            lat=slice(latn,lats),
            lon=slice(lon1,lon2)).vwnd<jupyter_output><empty_output><jupyter_text>Step 2: 利用MetPy計算渦度。<jupyter_code>import metpy.calc as mpcalc
vor = mpcalc.vorticity(u, v)
vorm = vor.mean(axis=0)
vorm<jupyter_output><empty_output><jupyter_text>Step 3: 繪圖。<jupyter_code>import cmaps

plt.figure()                  
proj = ccrs.PlateCarree()     
ax = plt.axes(projection=proj)
ax.set_extent([lon1,lon2,lats,latn],crs=proj)
ax.set_xticks(np.arange(80,180,20), crs=proj)
ax.set_yticks(np.arange(-20,40,10), crs=proj)
lon_formatter = LONGITUDE_FORMATTER
lat_formatter = LATITUDE_FORMATTER   
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter) 
ax.coastlines()

plt.title("850-hPa vorticity, Dec. 2017", loc='left')    
vorm = vorm * 10e5
vorPlot = (vorm.plot.contourf("lon", "lat",  
                              ax=ax,       
                              levels=np.arange(-20,24,4), 
                              cmap=cmaps.CBR_coldhot,  
                              add_colorbar=True, 
                              extend='both', 
                              cbar_kwargs={'orientation': 'horizontal', 'aspect': 30, 'label': ' '}) #設定color bar
                              )
                              
ax.set_title(' ')
ax.set_title('850-hPa Vorticity, December 2017', loc='left')
ax.set_ylabel(' ')             # 設定坐標軸名稱。
ax.set_xlabel(' ')
plt.show()<jupyter_output><empty_output><jupyter_text>```{admonition} Exercise:class: tip同Example 5，但計算散度 [metpy.calc.divergence](https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.divergence.html)。```  計算流函數(streamfunction)、速度位(velocity potential)使用`windspharm`套件，範例請見[官方網站](https://ajdawson.github.io/windspharm/latest/examples/sfvp_xarray.html)。`windspharm`也可以用來計算* divergence* vorticity (relative and absolute)* streamfunction* velocity potential* irrotational and non-divergent components of the wind (Helmholtz decomposition)* vector gradient of a scalar function* triangular truncation of a scalar field* magnitude (wind speed)可以試著比較用`MetPy`計算渦度、散度和用`windspharm`計算有沒有什麼不同。  計算經驗正交函數 (EOF)使用`eofs`套件，範例請見[官方網站](https://ajdawson.github.io/eofs/latest/examples/xarray_examples_index.html)。<jupyter_code>from eofs.xarray import Eof<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14543035513253505, 'fraction_numerical': 0.0421655920581312, 'mean_word_length': 3.2450043440486533, 'pattern_counts': {'":': 1, '<': 63, '<?xml version=': 0, '>': 63, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '294370', 'license': 'Unlicense', 'n_tokens_mistral': 4917, 'n_tokens_neox': 4843, 'n_words': 450, 'path': 'pyaos-wks/docs/_sources/07_statistics_calc.ipynb', 'repo_name': 'wyhtsai/wyhtsai.github.io', 'size': '9771'}