[ncl-talk] problem with specx_anal (ntim x nlat x nlon)
Dave Allured - NOAA Affiliate
dave.allured at noaa.gov
Thu Aug 10 21:44:17 MDT 2023
Here is an untested modification of your second posted version. This
should solve the missing value and dimension problems. This avoids using
ndtooned, which I think is unnecessary and complicated here. Please add
the necessary setup and reader code at the top.
spcavg = new (ntim/2, typeof(rain))
spcavg = 0.0
r1zsum = 0.0
nvalid = 0
do nl=0,nlat-1
do ml=0,mlon-1
seg = rain(latitude|nl,longitude|ml,time|:) ; extract 1-D grid
point
if (.not. any (ismissing (seg))) then
nvalid = nvalid + 1
dof = specx_anal (seg, 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 if
end do ; ml
print ("nl, nvalid = " + nl + " " + nvalid)
end do ; nl
r1z = r1zsum/nvalid ; average r1z
spcavg = spcavg/nvalid ; average spectrum
On Thu, Aug 10, 2023 at 8:32 PM Lyndz via ncl-talk <
ncl-talk at mailman.ucar.edu> wrote:
> 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
>>
>>
>> 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
>>>>>
>>>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20230810/017d0d27/attachment.htm>
More information about the ncl-talk
mailing list