<div dir="ltr"><div dir="ltr">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.</div><div dir="ltr"><br></div><div>spcavg = new (ntim/2, typeof(rain))<br>spcavg = 0.0<br>r1zsum = 0.0<br></div><div>nvalid = 0</div><div dir="ltr"><br></div><div dir="ltr">do nl=0,nlat-1</div><div dir="ltr"> do ml=0,mlon-1<br></div><div dir="ltr"> seg = rain(latitude|nl,longitude|ml,time|:) ; extract 1-D grid point</div><div dir="ltr"> if (.not. any (ismissing (seg))) then</div><div dir="ltr"> nvalid = nvalid + 1<br> dof = specx_anal (seg, d, sm, pct) ; current segment spc<br> spcavg = spcavg + dof@spcx ; sum spc of each segment<br> r1 = dof@xlag1 ; extract segment lag-1<br> r1zsum = r1zsum + 0.5*(log((1+r1)/(1-r1))) ; sum the Fischer Z</div><div dir="ltr"> end if<br> end do ; ml</div><div dir="ltr"> print ("nl, nvalid = " + nl + " " + nvalid)</div><div dir="ltr">end do ; nl<br><br>r1z = r1zsum/nvalid ; average r1z<br>spcavg = spcavg/nvalid ; average spectrum</div><div><br></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Thu, Aug 10, 2023 at 8:32 PM Lyndz via ncl-talk <<a href="mailto:ncl-talk@mailman.ucar.edu">ncl-talk@mailman.ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Dear Sir Dennis, NCL-experts.<br></div><div><br></div><div>Thank you for this valuable information!</div><div><br></div><div>Do you have any suggestions on how to solve this issue?</div><div><b>fatal:specx_anal: 'x' cannot contain any missing values<br>fatal:["Execute.c":8637]:Execute: Error occurred at or near line 69 in file spectral_analysis.ncl</b></div><div><b><br></b></div><div><b>from this line:</b></div><div> dof = specx_anal(rain(n*nlen:(n+1)*nlen-1,nl,ml),dm,sm,pct) ; assuming rain has dims as: (time, lat, lon)</div><div><br></div><div>I checked the diagnostic functions from the MJO page(<a href="https://www.ncl.ucar.edu/Document/Functions/Diagnostics/index.shtml" target="_blank">https://www.ncl.ucar.edu/Document/Functions/Diagnostics/index.shtml</a>).</div><div>The closest one that I was trying to do is using the <a href="https://www.ncl.ucar.edu/Document/Functions/Diagnostics/mjo_spectra_season.shtml" target="_blank">
mjo_spectra_season</a>.</div><div>But this one is for area averaged spectra. I would like to do the spectral analysis for all grid points with valid values.</div><div><br></div><div>Thank you for all the help.</div><div><br></div><div>Sincerely,</div><div>Lyndz<br></div><div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Fri, Aug 11, 2023 at 9:24 AM Dennis Shea <<a href="mailto:shea@ucar.edu" target="_blank">shea@ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>[1] The spectral analysis function uses an FFT. Missing values [_FillValue] are not allowed. <br></div><div>[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.<br></div><div>[3] My recollection of the way NCL executes the <a href="https://www.ncl.ucar.edu/Document/Functions/Built-in/specx_anal.shtml" target="_blank"><b>specx_anal</b></a><b> function </b>is <br></div><div> [a] the <b>time series for</b> <b>each individual grid point </b>is examined prior to invoking the FFT<br></div><div> [b] If the grid point time series contains one or more missing values, an error counter is incremented: ERR=ERR+1<br></div><div> [c] If no missing values are present, the spectral information for that grid point is computed.</div><div> [d] <b>after</b> 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.</div><div><br></div><div>Good Luck<br></div><div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Thu, Aug 10, 2023 at 1:47 AM Lyndz via ncl-talk <<a href="mailto:ncl-talk@mailman.ucar.edu" target="_blank">ncl-talk@mailman.ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Dear Rashed, NCL-experts,</div><div><br></div><div>Many thanks for this. <br></div><div>My goal is actually to plot the fractional variance as a map after getting the mean spectra, <br></div><div>Something like this, which I can get from the average spectra: <br><p>5-to-10day variance: var_510 = SUM{<i>spcx</i>('f5->f10')*df} <br></p><p>fractional_variance = var_510/total_variance</p></div><div>Attached is the modified script following your suggestion and I got this error:</div><div>fatal:specx_anal: 'x' cannot contain any missing values</div><div>fatal:["Execute.c":8637]:Execute: Error occurred at or near line 69 in file spectral_analysis.ncl</div><div><br></div><div>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.<br></div><div>To remove this error, I set:</div><div><br></div><div>rain@_FillValue = 0.0<br>delete(rain@_FillValue)</div><div><br></div><div><b><span style="color:rgb(255,0,0)">Is there another way to avoid this? Like ignore the gridpoints with missing values and the input/output will still have the 3D dimension?</span></b><br></div><div><br></div><div>If I do something like this at the beginning, the input is no longer 3D.</div><div><br></div><div>p1 = ndtooned( rain )
<br>
i = ind(.not.ismissing(p1)) ; all non-missing values
<br>
if (.not.all(ismissing(i))) then
<br>
rainNew = p1(i) ; rainNew is 1D
<br>
copy_VarAtts (rain, rainNew) ; copy attributes
<br>
else
<br>
print("No missing values")
<br>
end if
</div><div><br></div><div>--Lynz<br></div><div><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><table style="border-collapse:collapse" cellspacing="0" cellpadding="0" border="0"><tbody><tr><td style="width:134.7pt;border-color:currentcolor rgb(191,143,0) currentcolor currentcolor;border-style:none solid none none;border-width:medium 1.5pt medium medium;padding:0cm 5.4pt" width="180" valign="top"><br></td><td style="width:496.1pt;border:medium;padding:0cm 5.4pt" width="661" valign="top"><br></td></tr></tbody></table></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Aug 9, 2023 at 5:13 PM Rashed Mahmood <<a href="mailto:rashidcomsis@gmail.com" target="_blank">rashidcomsis@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Hi Lyndz,</div><div>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.</div><div><br></div><div>begin<br>ntim = 6954<br>nlat = 220<br>nlon = 220<br>;rain = random_normal(6.5,19.8,(/ntim,nlat,nlon/))<br>nseg = 49<br>nlen = ntim/nseg<br>minlat = -14.875<br>maxlat = 39.875<br>minlon = 90.125<br>maxlon = 144.875<br>spcavgFin = new((/nlen/2,nlat,nlon/),float) ;onedtond(rain,(/nseg,nlat,nlon/))<br>r1zsumFin = new((/nlat,nlon/),float)<br>printVarSummary(spcavgFin)<br>dm = 1<br>sm = 1<br>pct = 0.10<br>;r1zsum = 0.0<br><br> mlon=220<br> nlat=220<br> do nl=0,nlat-1<br> do ml=0,mlon-1<br> r1zsum =0.<br> spcavg =new(nlen/2,float)<br> spcavg =0.<br> do n=0,nseg-1<br> ;print(n+"")<br> ;dof = specx_anal(rain(latitude|nl:nl,longitude|ml:ml,time|:),d,sm,pct) ; current segment spc<br> dof = specx_anal(rain(n*nlen:(n+1)*nlen-1,nl,ml),dm,sm,pct) ; assuming rain has dims as: (time, lat, lon) <br> spcavg = spcavg + dof@spcx ; sum spc of each segment<br> r1 = dof@xlag1 ; extract segment lag-1<br> r1zsum = r1zsum + 0.5*(log((1+r1)/(1-r1))) ; sum the Fischer Z<br> end do<br> r1z = r1zsum/nseg ; average r1z<br> r1 = (exp(2*r1z)-1)/(exp(2*r1z)+1) ; transform back<br> r1zsumFin(nl,ml) = r1 ; this is the mean r1<br> spcavgFin(:,nl,ml) = spcavg/nseg ; average spectrum<br> delete([/r1z,r1/])<br> end do ; ml<br> end do ; nl<br>end<br></div><div><br></div><div>Good luck.</div><div>Rashed<br></div><div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Aug 9, 2023 at 9:16 AM Lyndz via ncl-talk <<a href="mailto:ncl-talk@mailman.ucar.edu" target="_blank">ncl-talk@mailman.ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Dear NCL-experts,</div><div><br></div><div>Just a follow up on my post earlier.</div><div><br></div><div>I tried the following:</div><div><br></div><div>;************************************************<br>; calculate mean spectrum spectrum and lag1 auto cor<br>;************************************************<br><br>ntim = 6954<br>nlat = 220<br>nlon = 220<br>nseg = 49<br>nlen = ntim/nseg<br>minlat = -14.875<br>maxlat = 39.875<br>minlon = 90.125<br>maxlon = 144.875<br>spcavg = onedtond(rain,(/nseg,nlat,nlon/))<br>printVarSummary(spcavg)<br>spcavg = 0.0<br>r1zsum = 0.0<br><br>dm = 1<br>sm = 1<br>pct = 0.10<br>r1zsum = 0.0<br><br>mlon=220<br>nlat=220<br><br> do nl=0,nlat-1<br> do ml=0,mlon-1<br> do n=0,nseg-1<br> dof = specx_anal(rain(latitude|nl:nl,longitude|ml:ml,time|:),d,sm,pct) ; current segment spc<br> spcavg = spcavg + dof@spcx ; sum spc of each segment<br> r1 = dof@xlag1 ; extract segment lag-1<br> r1zsum = r1zsum + 0.5*(log((1+r1)/(1-r1))) ; sum the Fischer Z<br> end do<br> end do ; ml<br>end do ; nl<br></div><div><br></div><div>It gives the following error:</div><div><br></div><div><span style="color:rgb(255,0,0)">ncl 64> do nl=0,nlat-1 <br>ncl 65> do ml=0,mlon-1<br>ncl 66> do n=0,nseg-1<br>ncl 67> dof = specx_anal(rain(latitude|nl:nl,longitude|ml:ml,time|:),d,sm,pct) ; current segment spc<br>ncl 68> spcavg = spcavg + dof@spcx ; sum spc of each segment<br>ncl 69> r1 = dof@xlag1 ; extract segment lag-1<br>ncl 70> r1zsum = r1zsum + 0.5*(log((1+r1)/(1-r1))) ; sum the Fischer Z<br>ncl 71> end do<br>ncl 72> end do ; ml<br>ncl 73> end do ; nl</span></div><div><br></div><div><span style="color:rgb(255,0,0)"><b>fatal:Number of dimensions in parameter (0) of (specx_anal) is (3), (1) dimensions were expected <br>fatal:["Execute.c":8635]:Execute: Error occurred at or near line 67</b></span></div><div><br></div><div>Any suggestion on how to do this correctly in NCL?<br></div><div><br></div><div>-Lyndz<br></div><div><br></div><div class="gmail_quote"><br><div dir="ltr"><div>Dear NCL-experts,</div><div><br></div><div>I would like to do a spectral analysis on a <b>ntim x nlat x lon</b> data as indicated in example 3 from this page: <a href="https://www.ncl.ucar.edu/Document/Functions/Built-in/specx_anal.shtml" target="_blank">https://www.ncl.ucar.edu/Document/Functions/Built-in/specx_anal.shtml</a></div><div><br></div><div><b>Specifically, I would like to get the mean spectra so that the final output will have meanspectra x lat x lon</b></div><div><br></div><div><b>So far, I have the following:</b></div><div><br></div><div> ufile = addfile("<a href="http://anom_rain_sm_1951-2007_jjas.nc" target="_blank">anom_rain_sm_1951-2007_jjas.nc</a>","r") ; open netcdf file<br> rain = ufile->rAnom_sm(:,:,:)<br><br> printVarSummary(rain) ;[time | 6954] x [latitude | 220] x [longitude | 220]<br><br> d = 0<br> sm = 1 ; periodogram<br> pct = 0.10<br><br>;************************************************<br>; calculate mean spectrum spectrum and lag1 auto cor<br>;************************************************<br><br>ntim = 6954<br>nlat = 220<br>nlon = 220<br><br>nseg = 49<br>nlen = ntim/nseg<br><br>spcavg = new (ntim/141, typeof(rain))<br>printVarSummary(spcavg)<br>spcavg = new ( ntim/2, typeof(rain))<br>spcavg = 0.0<br>r1zsum = 0.0<br><br> <span style="color:rgb(255,0,0)"> do n=0,nseg-1<br> dof = specx_anal(rain(latitude|nl:nl,longitude|ml:ml,time|:),d,sm,pct) ; current segment spc<br> spcavg = spcavg + dof@spcx ; sum spc of each segment<br> r1 = dof@xlag1 ; extract segment lag-1<br> r1zsum = r1zsum + 0.5*(log((1+r1)/(1-r1)) ; sum the Fischer Z<br> end do<br> r1z = r1zsum/nseg ; average r1z<br> r1 = (exp(2*r1z)-1)/(exp(2*r1z)+1) ; transform back<br> ; this is the mean r1<br> spcavg = spcavg/nseg ; average spectrum<br></span></div><div><b><br></b></div><div><b>Questions:</b><br></div><div> (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.</div><div><br></div><div>I'll appreciate any help on this.</div><div><br></div><div>-Lyndz</div></div></div></div></blockquote></div></blockquote></div></blockquote></div></blockquote></div>
</blockquote></div></div>