Re: Spectral analysis problems

[a]
You have:
spec = specx_anal(u(lev|kl:kl,lat|nl:nl,lon|ml:ml,time:),d,sm,pct)
the "time:" should be "time|:"
spec = specx_anal(u(lev|kl:kl,lat|nl:nl,lon|ml:ml,time|:),d,sm,pct)
possibly clearer
spec = specx_anal(u(lev|kl, lat|nl, lon|ml ,time|:),d,sm,pct

This will yield the temporal power spectrum the the specified grid point.

===============
As far as finding the major *temporal* periodicities in the data, the
fact that your data are cyclic in longitude is immaterial.

If you are looking at *spatial* wave patterns in longitude, then,
because the data are periodic [cyclic] in longitude, you can use
a simple fft
http://www.ncl.ucar.edu/Document/Functions/Built-in/ezfftf.shtml

perhaps better, see the fourier_info function:
http://www.ncl.ucar.edu/Document/Functions/Built-in/fourier_info.shtml

Good luck
Helen Parish wrote:
> I am trying to do spectral analysis on a netcdf file.
>
> a. Firstly I want to find out what the major periodicities are in the
> data (these may vary from months to several years). I am assuming they
> are periodic in longitude. I would like to plot the range of periods
> which are present, as a function of pressure and / or height at a
> given latitude. I would also like to plot the period present as a
> function of latitude at a given height or pressure. There are a lot of
> possible routines available, and I am not sure if I am using the right
> one(s) (see below).
>
> b. I would also like to be able to find a way to follow a specific
> wave as it travels upwards (or downwards), and to plot this. I am not
> sure what the best routine to use would be, and how to plot this?.
>
>
> in order to do a. above, I have tried to incorporate specx_anal into a
> script, but I am not sure how I should be inputting the arguments, and
> whether this will do what I want if I succeed ?. The error I am
> getting is below, and underneath that is the script I have been trying
> to use (please bear in mind that my pressure levels apply to Venus,
> and the magnitudes are correct).
>
> I am not sure how to do b. above/
>
> Thanks,
> Helen.
>
> Here is the error message:
>
> (0) T: interpolated and written to netCDF
> (0) Q: interpolated and written to netCDF
> (0) U: interpolated and written to netCDF
> (0) V: interpolated and written to netCDF
> fatal:syntax error: line 76 in file four.ncl before or near time
> spec = specx_anal(u(lev
> -------------------------------------------------------^
>
> fatal:Error in subscript, named subscripting is being used, make sure
> each subscript has a name
> fatal:syntax error: function specx_anal expects 4 arguments, got 0
> fatal:error at line 76 in file four.ncl
>
> fatal:Syntax Error in block, block not executed
> fatal:error at line 114 in file four.ncl
>
>
> Here is the script I have been trying to use:
>
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
> begin
> ; input file
> diri = "./"
> fili = "rand09125k.cam2.h0.0026-01-17-00000.nc"
> fi = addfile (diri+fili, "r")
> ; output file
> diro = "./"
> filo = "helen.nc"
> system ("/bin/rm -f "+diro+filo) ; remove any pre-exist file
> fo = addfile (diro+filo, "c")
> ; add any file attributes
> fo_at_title = fi_at_title + ": Selected Variables at pressure levels"
> fo_at_history = systemfunc ("date")
> fo_at_source = fi_at_source
> fo_at_case = fi_at_case
> fo_at_Conventions = fi_at_Conventions
>
> Var = (/ "T" , "Q", "U", "V"/) ; select variable to be interpolated.
> nVar = dimsizes (Var)
> ; desired output levels
> ;lev_p = (/ 91500., 90000., 87500., 80000., 67500., 50000., 30000.,
> 20000., 10000., 5000., 3000., 1500., 750., 390., 200., 100.,
> 50.,25., 20., 15.,10., 8., 5.,3.,1. /)
> ;lev_p = (/ 91500., 90000., 87500., 80000., 67500., 50000., 30000.,
> 20000., 10000., 5000., 3000., 1500., 750., 390., 200., 100.,
> 50.,25., 14., 7., 3.5, 1.9, 0.95, 0.4, 0.1, 0.01 /)
> lev_p = (/ 91500., 90000., 87500., 80000., 67500., 50000., 30000.,
> 20000., 10000., 5000., 3000., 1500., 750., 390., 200., 100.,
> 50.,25., 14., 7., 3.5, 1.9, 0.95, 0.4, 0.1, 0.03 /)
> ;lev_p = (/ 91500., 90000., 87500., 80000., 67500., 50000., 30000.,
> 20000., 10000., 5000., 3000., 1500., 750., 390., 200., 100.,
> 50.,25., 14., 7., 3.5, 2.0 /)
> lev_p!0 = "lev_p" ; variable and dimension name the same
> lev_p&lev_p = lev_p ; create coordinate variable
> lev_p_at_long_name = "pressure" ; attach some attributes
> lev_p_at_units = "hPa"
> lev_p_at_positive = "down"
>
> hyam = fi->hyam ; read hybrid info
> hybm = fi->hybm
> PS = fi->PS
> P0mb = 0.01*fi->P0
>
> do n=0,nVar-1 ; loop over the variables
> X = fi->$Var(n)$
> Xp = vinth2p (X, hyam, hybm, lev_p ,PS, 1, P0mb, 2, True)
> copy_VarAtts(X, Xp)
> fo->$Var(n)$ = Xp ; write to netCDF file
> print (Var(n)+": interpolated and written to netCDF")
> end do
>
> end
> ;***********************
> ; zonal.ncl
> ;***********************
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
> ;***********************
> begin
> f = addfile ("helen.nc","r")
> print(f)
> u = f->U
> printVarSummary( u )
>
> dimu = dimsizes( u )
> ntim = dimu(0)
> klvl = dimu(1)
> nlon = dimu(2)
> mlon = dimu(3)
>
> uavg = dim_avg_Wrap(u)
>
> ;***********************
> ; Create Plot
> ;***********************
> d = 0
> sm = 3
> pct = 0.1
>
> kl = 10
> nl = 95
> ml = 40
> spec = specx_anal(u(lev|kl:kl,lat|nl:nl,lon|ml:ml,time:),d,sm,pct)
>
> wks = gsn_open_wks ("pdf", "rand09125k260117spec" ) ; open workstation
> gsn_define_colormap(wks,"rainbow") ; choose colormap
>
> res = True
> res_at_cnFillOn = True
> res_at_lbLabelAutoStride = True
> res_at_gsnMaximize = True ; if [ps, eps, pdf] make large
> res_at_gsnSpreadColors = True ; span color map
>
> do nt=0,ntim-1
>
> res_at_gsnCenterString = "rand09125k260117spec"
> res_at_tiXaxisString = "Frequency"
> res_at_tiYaxisString = "Variance"
>
> plot = gsn_csm_xy(wks,spec_at_frq,spec_at_spcx,res)
>
> ; plot = gsn_csm_pres_hgt(wks, u(nt,:,:,40), res ) ; (lev,lat)
> ;; plot = gsn_csm_pres_hgt(wks, uavg(nt,:,:), res ) ; (lev,lat)
> ; plot = gsn_csm_contour (wks, uavg(t,:,:), res ) ; (lev,lat)
>
> res_at_trYReverse = True
> ; plot = gsn_csm_contour (wks, u(nt,:,:,0), res ) ; (lev,lat)
> ; plot = gsn_csm_contour (wks, u(nt,:,0,:), res ) ; (lev,lon)
>
> kl = 5
> res_at_mpFillOn = False
> res_at_mpGridAndLimbOn = True
> res_at_mpGridLineDashPattern = 2
> res_at_mpOutlineBoundarySets = "NoBoundaries"
> res_at_mpCenterLonF = 180.
> res_at_gsnCenterString = res_at_gsnCenterString+" p="+u&lev_p(kl)
> ; plot = gsn_csm_contour_map_ce (wks, u(nt,kl,:,:), res )
>
> end do
>
> end
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

-- 
======================================================
Dennis J. Shea                  tel: 303-497-1361    |
P.O. Box 3000                   fax: 303-497-1333    |
Climate Analysis Section                             |
Climate & Global Dynamics Div.                       |
National Center for Atmospheric Research             |
Boulder, CO  80307                                   |
USA                        email: shea 'at' ucar.edu |
======================================================
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk