{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 利用NCL做FFT-頻譜分析"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Example 1:** 菲律賓中部 (9˚-14˚N, 122˚-127˚E) 在DJF季節區域平均降雨的頻譜分析。"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"vscode": {
"languageId": "ncl"
}
},
"outputs": [
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
";*************************************************\n",
"; spec_3.ncl\n",
";\n",
"; Concepts illustrated:\n",
"; - Calculating confidence intervals\n",
";************************************************\n",
";\n",
"; These files are loaded by default in NCL V6.2.0 and newer\n",
"; load \"$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl\"\n",
"; load \"$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl\" \n",
"; load \"$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl\"\n",
";\n",
"; This file still has to be loaded manually\n",
"load \"$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl\"\n",
";************************************************\n",
"begin\n",
";************************************************\n",
"; variable and file handling\n",
";************************************************\n",
" fn = \"./data/cmorph_sample.nc\" ; define filename\n",
" in = addfile(fn,\"r\") ; open netcdf file\n",
" pcp = in->cmorph(:,{9:14},{122:127}) \n",
" pcpm = dim_avg_n_Wrap(dim_avg_n_Wrap(pcp, 2), 1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"NCL的時間格式必須是yyyymmdd，我們要先做一些轉換。"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"vscode": {
"languageId": "ncl"
}
},
"outputs": [],
"source": [
"rtime = pcpm&time \n",
"delete(pcpm&time) \n",
"pcpm!0 = \"time\"\n",
"pcpm&time = cd_calendar(rtime,2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"接下來一年一年進行頻譜分析，然後再對年做平均。"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"vscode": {
"languageId": "ncl"
}
},
"outputs": [],
"source": [
";************************************************\n",
"; set function arguments\n",
";************************************************\n",
" d = 0 ; detrending opt: 0=>remove mean 1=>remove mean + detrend\n",
" sm = 1 ; smooth: should be at least 3 and odd\n",
" pct = 0.10 ; percent taper: (0.0 <= pct <= 1.0) 0.10 common. \n",
" ys = 1998 ; Start year \n",
" ye = 2020 ; End year\n",
" ntim = 90 ; # of days in DJF\n",
"\n",
";************************************************\n",
"; calculate mean spectrum spectrum and lag1 auto cor\n",
";************************************************\n",
"\n",
"; loop over each segment of length ntim\n",
" \n",
" spcavg = new ( ntim/2, typeof(pcpm))\n",
" spcavg = 0.0\n",
"\n",
" r1zsum = 0.0\n",
" \n",
" do yy = ys, ye\n",
" time1 = yy *10000 + 1201 \n",
" time2 = (yy+1)*10000 + 228\n",
" pcp_5d_ave = runave_n_Wrap(pcpm({time1:time2}), 5, 1, 0)\n",
" dof = specx_anal(pcp_5d_ave,d,sm,pct) ; current segment spc\n",
" spcavg = spcavg + dof@spcx ; sum spc of each segment\n",
" r1 = dof@xlag1 ; extract segment lag-1\n",
" r1zsum = r1zsum + 0.5*(log((1+r1)/(1-r1))) ; sum the Fischer Z\n",
" end do\n",
" \n",
" r1z = r1zsum/(ye-ys+1) ; average r1z\n",
" r1 = (exp(2*r1z)-1)/(exp(2*r1z)+1) ; transform back\n",
" ; this is the mean r1\n",
" spcavg = spcavg/(ye-ys+1) ; average spectrum\n",
"\n",
";************************************************\n",
"; Assign mean spectrum to data object\n",
";************************************************\n",
" \n",
" df = 2.0*(ye-ys+1) ; deg of freedom\n",
" ; all segments\n",
"df@spcx = spcavg ; assign the mean spc\n",
"df@frq = dof@frq\n",
"df@xlag1= r1 ; assign mean lag-1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"接下來繪圖。"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"vscode": {
"languageId": "ncl"
}
},
"outputs": [],
"source": [
";************************************************\n",
"; plotting parameters\n",
";************************************************\n",
" wks = gsn_open_wks(\"png\",\"./images/spec\") ; send graphics to PNG file\n",
"\n",
" r = True ; plot mods desired\n",
" r@gsnDraw = False ; do not draw\n",
" r@gsnFrame = False ; do not advance frame\n",
" r@tiMainString = \"DJF Precip. @ PHL-C\" ; title\n",
" r@tiXAxisString = \"Period (days)\" ; xaxis\n",
" r@tiYAxisString = \"Variance\" ; yaxis\n",
" r@trXMaxF = 80.\n",
" r@trXMinF = 5.\n",
";************************************************\n",
"; plot\n",
";************************************************\n",
"r@xyLineThicknesses = (/2.,1.,1.,1./) ; Define line thicknesses \n",
"r@xyDashPatterns = (/0,0,1,1/) ; Dash patterns\n",
"\n",
";***********************************************\n",
"; Plot: Generate fancier plot showing \n",
"; \"red noise\" confidence bounds\n",
"; (a) solid for spectrum and Markov, \n",
"; (b) dash for bounds \n",
";***********************************************\n",
" p = 1/df@frq ; *highly non-linear*\n",
" p!0 = \"f\"\n",
" p&f = df@frq\n",
" p@long_name = \"period\"\n",
" p@units = \"day\"\n",
"\n",
" ip = ind(p.le.90 .and. p.ge.5) ; all indices for \"short\" periods\n",
" splt = specx_ci(df, 0.05, 0.95) ; calc confidence interval\n",
" plot = gsn_csm_xy(wks,p(ip), splt(:,ip),r)\n",
" draw(plot)\n",
" frame(wks)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"vscode": {
"languageId": "ncl"
}
},
"outputs": [],
"source": [
"end"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "p3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"version": "3.10.10"
},
"orig_nbformat": 4,
"terminalId": 41396
},
"nbformat": 4,
"nbformat_minor": 2
}