<div dir="ltr"><div>It is not clear to me what:</div><div>"obtain a power spectrum of daily rainfall data for 24 years (for a specific season having 122 days/year)"</div><div><br></div><div>[1] 24 yearly segments of 122 days ?</div><div>[2] soi = in->PAVE <br></div><div> Are the 'PAVE' input as 122*25 days? <br></div><div> As always, ncl-talk asks that users include the output from (here)</div><div> printVarSummary(soi) <br></div><div><br></div><div>--</div><div>Note: You should NOT perform a single spectrum on a series of 24*122 days.</div><div> You must perform segment averaging.</div><div> (a) compute the raw spectrum of each 122 day segment and then average <br></div><div> (b) generally people do not due spectral analysis on series which contain many zeros like you may have with daily precipitation.</div><div><br></div><div><br></div><div> <br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Mon, Oct 19, 2020 at 12:00 AM Nishtha Agrawal via ncl-talk <<a href="mailto:ncl-talk@mailman.ucar.edu">ncl-talk@mailman.ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Hi,</div><div><br></div><div>I am trying to obtain a power spectrum of daily rainfall data for 24 years (for a specific season having 122 days/year). I want to obtain the spectra with respect to period instead of frequency. For this I modified the available code that uses specx_anal and is given below. I am facing a strange issue where the values of period start from 60. I don't get any value before the 60 day period. I don't understand if I am making some mistake in the following code</div><div><br></div><div> fn = "<a href="http://data.nc" target="_blank">data.nc</a>" ; define filename<br> in = addfile(fn,"r") ; open netcdf file<br> soi = in->PAVE ; get data<br></div><div> d = 1 ; detrending opt: 0=>remove mean 1=>remove mean + detrend<br> sm = 21 ; smooth: should be at least 3 and odd<br> pct = 0.10 ; percent taper: (0.0 <= pct <= 1.0) 0.10 common.</div><div> sdof = specx_anal(soi,d,sm,pct)</div><div> splt = specx_ci (sdof, 0.10, 0.90)</div><div> wks = gsn_open_wks("png","psd_rainfall") ; send graphics to PNG file<br><br> r = True ; plot mods desired<br> r@tiMainString = "RegCM" ; title<br> r@tiXAxisString = "Frequency (cycles/day)" ; xaxis<br> r@tiYAxisString = "Variance/frq_interval" ; yaxis<br> r@xyLineColors = (/"foreground","red","blue","blue"/)</div><div> f = sdof@frq/30.5 ; for daily data<br> r@trYMinF = 0.00 ; manually set lower limit<br> r@trYMaxF = 3.0 ; " upper<br> r@gsnFrame = False ; do not advance frame<br> plot = gsn_csm_xy(wks,f, splt,r)<br><br> xf = (/0.40, 0.40+sdof@bw/) ; Create band width line<br> ys = 0.75*max(sdof@spcx) ; 75% up Y axis<br> yv = (/ys,ys/)<br> rpl = True ; resources for polyline<br> rpl@gsLineThicknessF = 2 ; Define line thickness<br> gsn_polyline(wks,plot,xf,yv,rpl) ; Draw BandWidth<br><br> txres= True ; label BW line<br> txres@txFontHeightF = 0.015 ; font height<br> txres@txJust = "CenterLeft" ; Set lable location<br> gsn_text(wks,plot,"BW",0.41+sdof@bw,ys,txres); Label<br> frame (wks)<br><br> r@gsnFrame = True ;reset to default</div><div> p = 1/f ; *highly non-linear*<br> p!0 = "f"<br> p&f = f<br> p@long_name = "period"<br> p@units = "day"<br><br> ;;print("====")<br> ;;print(f+" "+p+" "+splt(0,:) ) ; print<br> ;;print("====")<br><br> r@tiXAxisString = "Period (day)" ; xaxis<br> r@tiYAxisString = "Variance/freq" ; yaxis<br><br> ip = ind(p.le.120) ; all indices for "short" periods<br> plot = gsn_csm_xy(wks,p(ip), splt(:,ip),r)<br></div></div>
_______________________________________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@mailman.ucar.edu" target="_blank">ncl-talk@mailman.ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="https://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">https://mailman.ucar.edu/mailman/listinfo/ncl-talk</a></blockquote></div>