[ncl-talk] Fwd: problem with specx_anal (ntim x nlat x nlon)
Rashed Mahmood
rashidcomsis at gmail.com
Wed Aug 9 03:13:04 MDT 2023
Hi Lyndz,
There are several problems with your code. I do not have time to explain
these issues but you can compare the two codes to see the differences. The
below code should work, however, you need to think about what output you
are expecting and do testing of the results. I am not an expert on spectral
analysis so you will need to make sure if the results are as expected.
begin
ntim = 6954
nlat = 220
nlon = 220
;rain = random_normal(6.5,19.8,(/ntim,nlat,nlon/))
nseg = 49
nlen = ntim/nseg
minlat = -14.875
maxlat = 39.875
minlon = 90.125
maxlon = 144.875
spcavgFin = new((/nlen/2,nlat,nlon/),float)
;onedtond(rain,(/nseg,nlat,nlon/))
r1zsumFin = new((/nlat,nlon/),float)
printVarSummary(spcavgFin)
dm = 1
sm = 1
pct = 0.10
;r1zsum = 0.0
mlon=220
nlat=220
do nl=0,nlat-1
do ml=0,mlon-1
r1zsum =0.
spcavg =new(nlen/2,float)
spcavg =0.
do n=0,nseg-1
;print(n+"")
;dof =
specx_anal(rain(latitude|nl:nl,longitude|ml:ml,time|:),d,sm,pct) ;
current segment spc
dof = specx_anal(rain(n*nlen:(n+1)*nlen-1,nl,ml),dm,sm,pct) ;
assuming rain has dims as: (time, lat, lon)
spcavg = spcavg + dof at spcx ; sum spc of each segment
r1 = dof at xlag1 ; extract segment lag-1
r1zsum = r1zsum + 0.5*(log((1+r1)/(1-r1))) ; sum the Fischer Z
end do
r1z = r1zsum/nseg ; average r1z
r1 = (exp(2*r1z)-1)/(exp(2*r1z)+1) ; transform back
r1zsumFin(nl,ml) = r1 ; this is the mean r1
spcavgFin(:,nl,ml) = spcavg/nseg ; average spectrum
delete([/r1z,r1/])
end do ; ml
end do ; nl
end
Good luck.
Rashed
On Wed, Aug 9, 2023 at 9:16 AM Lyndz via ncl-talk <ncl-talk at mailman.ucar.edu>
wrote:
> Dear NCL-experts,
>
> Just a follow up on my post earlier.
>
> I tried the following:
>
>
> ;************************************************
> ; calculate mean spectrum spectrum and lag1 auto cor
> ;************************************************
>
> ntim = 6954
> nlat = 220
> nlon = 220
> nseg = 49
> nlen = ntim/nseg
> minlat = -14.875
> maxlat = 39.875
> minlon = 90.125
> maxlon = 144.875
> spcavg = onedtond(rain,(/nseg,nlat,nlon/))
> printVarSummary(spcavg)
> spcavg = 0.0
> r1zsum = 0.0
>
> dm = 1
> sm = 1
> pct = 0.10
> r1zsum = 0.0
>
> mlon=220
> nlat=220
>
> do nl=0,nlat-1
> do ml=0,mlon-1
> do n=0,nseg-1
> dof =
> specx_anal(rain(latitude|nl:nl,longitude|ml:ml,time|:),d,sm,pct) ;
> current segment spc
> spcavg = spcavg + dof at spcx ; sum spc of each segment
> r1 = dof at xlag1 ; extract segment lag-1
> r1zsum = r1zsum + 0.5*(log((1+r1)/(1-r1))) ; sum the Fischer Z
> end do
> end do ; ml
> end do ; nl
>
>
> It gives the following error:
>
> ncl 64> do nl=0,nlat-1
> ncl 65> do ml=0,mlon-1
> ncl 66> do n=0,nseg-1
> ncl 67> dof =
> specx_anal(rain(latitude|nl:nl,longitude|ml:ml,time|:),d,sm,pct) ;
> current segment spc
> ncl 68> spcavg = spcavg + dof at spcx ; sum spc of
> each segment
> ncl 69> r1 = dof at xlag1 ; extract
> segment lag-1
> ncl 70> r1zsum = r1zsum + 0.5*(log((1+r1)/(1-r1))) ; sum the
> Fischer Z
> ncl 71> end do
> ncl 72> end do ; ml
> ncl 73> end do ; nl
>
>
> *fatal:Number of dimensions in parameter (0) of (specx_anal) is (3), (1)
> dimensions were expected fatal:["Execute.c":8635]:Execute: Error occurred
> at or near line 67*
>
>
> Any suggestion on how to do this correctly in NCL?
>
> -Lyndz
>
>
>
>
>
>
>
>
>
>
>
>
> ---------------------------------------------------------------------
>
> Dear NCL-experts,
>
> I would like to do a spectral analysis on a *ntim x nlat x lon* data as
> indicated in example 3 from this page:
> https://www.ncl.ucar.edu/Document/Functions/Built-in/specx_anal.shtml
>
> *Specifically, I would like to get the mean spectra so that the final
> output will have meanspectra x lat x lon*
>
>
> *So far, I have the following:*
>
> ufile = addfile("anom_rain_sm_1951-2007_jjas.nc","r") ; open
> netcdf file
> rain = ufile->rAnom_sm(:,:,:)
>
> printVarSummary(rain) ;[time | 6954] x [latitude | 220] x
> [longitude | 220]
>
> d = 0
> sm = 1 ; periodogram
> pct = 0.10
>
>
> ;************************************************
> ; calculate mean spectrum spectrum and lag1 auto cor
> ;************************************************
>
> ntim = 6954
> nlat = 220
> nlon = 220
>
> nseg = 49
> nlen = ntim/nseg
>
> spcavg = new (ntim/141, typeof(rain))
> printVarSummary(spcavg)
> spcavg = new ( ntim/2, typeof(rain))
> spcavg = 0.0
> r1zsum = 0.0
>
> do n=0,nseg-1
> dof =
> specx_anal(rain(latitude|nl:nl,longitude|ml:ml,time|:),d,sm,pct) ;
> current segment spc
> spcavg = spcavg + dof at spcx ; sum spc of each segment
> r1 = dof at xlag1 ; extract segment lag-1
> r1zsum = r1zsum + 0.5*(log((1+r1)/(1-r1)) ; sum the Fischer Z
> end do
> r1z = r1zsum/nseg ; average r1z
> r1 = (exp(2*r1z)-1)/(exp(2*r1z)+1) ; transform back
> ;
> this is the mean r1
> spcavg = spcavg/nseg ; average spectrum
>
> *Questions:*
> (1) How do I loop this across all the gridpoints? I am having a problem
> with the looping and saving the spectra in spcavg(ntim,nlat,nlon). Any
> suggestions on how to do this in NCL? The example in the page is incomplete.
>
>
> I'll appreciate any help on this.
>
>
> -Lyndz
>
>
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at mailman.ucar.edu
> List instructions, subscriber options, unsubscribe:
> https://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20230809/fa2fbdd7/attachment.htm>
More information about the ncl-talk
mailing list