[ncl-talk] Fwd: problem with specx_anal (ntim x nlat x nlon)
Dennis Shea
shea at ucar.edu
Thu Aug 10 19:24:32 MDT 2023
[1] The spectral analysis function uses an FFT. Missing values [_FillValue]
are not allowed.
[2] You should *never* substitute a value [say, 0.0] for missing values.
You will get numbers returned but they are *not* the correct numbers. Value
substitution creates an erroneous spectrum.
[3] My recollection of the way NCL executes the *specx_anal*
<https://www.ncl.ucar.edu/Document/Functions/Built-in/specx_anal.shtml>*
function *is
[a] the *time series for* *each individual grid point *is examined
prior to invoking the FFT
[b] If the grid point time series contains one or more missing values,
an error counter is incremented: ERR=ERR+1
[c] If no missing values are present, the spectral information for
that grid point is computed.
[d] *after* looping through all the grid points, if 'ERR' is non zero,
the error message is issued to the user BUT the grid points with no
_FillValue should have useable numbers.
Good Luck
==========
%> less $NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl
On Thu, Aug 10, 2023 at 1:47 AM Lyndz via ncl-talk <
ncl-talk at mailman.ucar.edu> wrote:
> Dear Rashed, NCL-experts,
>
> Many thanks for this.
> My goal is actually to plot the fractional variance as a map after getting
> the mean spectra,
> Something like this, which I can get from the average spectra:
>
> 5-to-10day variance: var_510 = SUM{*spcx*('f5->f10')*df}
>
> fractional_variance = var_510/total_variance
>
> Attached is the modified script following your suggestion and I got this
> error:
> fatal:specx_anal: 'x' cannot contain any missing values
> fatal:["Execute.c":8637]:Execute: Error occurred at or near line 69 in
> file spectral_analysis.ncl
>
> Im using the Aphrodite rainfall data set with missing values over the
> ocean and some landpoints. spectral analysis cannot be applied to time
> series with missing values.
> To remove this error, I set:
>
> rain at _FillValue = 0.0
> delete(rain at _FillValue)
>
>
> *Is there another way to avoid this? Like ignore the gridpoints with
> missing values and the input/output will still have the 3D dimension?*
>
> If I do something like this at the beginning, the input is no longer 3D.
>
> p1 = ndtooned( rain )
> i = ind(.not.ismissing(p1)) ; all non-missing values
> if (.not.all(ismissing(i))) then
> rainNew = p1(i) ; rainNew is 1D
> copy_VarAtts (rain, rainNew) ; copy attributes
> else
> print("No missing values")
> end if
>
> --Lynz
>
>
>
>
> On Wed, Aug 9, 2023 at 5:13 PM Rashed Mahmood <rashidcomsis at gmail.com>
> wrote:
>
>> 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
>>>
>> _______________________________________________
> 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/20230810/81626efa/attachment.htm>
More information about the ncl-talk
mailing list