[ncl-talk] Fwd: problem with specx_anal (ntim x nlat x nlon)
Lyndz
olagueralyndonmark429 at gmail.com
Thu Aug 10 01:31:27 MDT 2023
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
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20230810/22d7b690/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: spectral_analysis.ncl
Type: application/octet-stream
Size: 2246 bytes
Desc: not available
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20230810/22d7b690/attachment.obj>
More information about the ncl-talk
mailing list