<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><br></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><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" data-smartmail="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-top:none;border-bottom:none;border-left:none;border-right:1.5pt solid rgb(191,143,0);padding:0cm 5.4pt" width="180" valign="top"><br></td><td style="width:496.1pt;border:none;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><br></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:1px solid 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><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:1px solid 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></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><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><br></div><div><br></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><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><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><b><br></b></div><div><b><br></b></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>;************************************************<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><br></div><div>I'll appreciate any help on this.</div><div><br></div><div><br></div><div>-Lyndz<br></div><div><br></div><div><br></div><div><br></div><div><br></div></div>
</div></div>
_______________________________________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@mailman.ucar.edu" target="_blank">ncl-talk@mailman.ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="https://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">https://mailman.ucar.edu/mailman/listinfo/ncl-talk</a><br>
</blockquote></div>
</blockquote></div>