[ncl-talk] Fwd: problem with specx_anal (ntim x nlat x nlon)
Lyndz
olagueralyndonmark429 at gmail.com
Thu Aug 10 20:16:18 MDT 2023
Dear Sir Dennis, NCL-experts.
Thank you for this valuable information!
Do you have any suggestions on how to solve this issue?
*fatal:specx_anal: 'x' cannot contain any missing
valuesfatal:["Execute.c":8637]:Execute: Error occurred at or near line 69
in file spectral_analysis.ncl*
*from this line:*
dof = specx_anal(rain(n*nlen:(n+1)*nlen-1,nl,ml),dm,sm,pct) ; assuming
rain has dims as: (time, lat, lon)
I checked the diagnostic functions from the MJO page(
https://www.ncl.ucar.edu/Document/Functions/Diagnostics/index.shtml).
The closest one that I was trying to do is using the mjo_spectra_season
<https://www.ncl.ucar.edu/Document/Functions/Diagnostics/mjo_spectra_season.shtml>
.
But this one is for area averaged spectra. I would like to do the spectral
analysis for all grid points with valid values.
Thank you for all the help.
Sincerely,
Lyndz
On Fri, Aug 11, 2023 at 9:24 AM Dennis Shea <shea at ucar.edu> wrote:
> [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/20230811/ca94b213/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: spectral_analysis.ncl
Type: application/octet-stream
Size: 2248 bytes
Desc: not available
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20230811/ca94b213/attachment.obj>
More information about the ncl-talk
mailing list