From: Adam Phillips <asphilli_at_nyahnyahspammersnyahnyah>

Date: Tue Aug 02 2011 - 16:52:01 MDT

Date: Tue Aug 02 2011 - 16:52:01 MDT

Hi Shinn,

Your coding is somewhat correct. There is no need to go through a double

do loop nor do you need to use the y-intercept. Your coding could be

condensed to:

trend = data(:,:,0)

temp = dtrend_msg(ispan(0,29,1),data,False,True)

trend = (/ ndtooned(temp@slope,(/dimsizes(data&lat),dimsizes(data&lon))*30/)

To calculate the significance of the trend, you should use the betainc

function:

http://www.ncl.ucar.edu/Document/Functions/Built-in/betainc.shtml

continuing from the coding above:

tval =new((/dimsizes(data&lat),dimsizes(data&lon/),typeof(data))

beta_b = tval

df = new((/dimsizes(data&lat),dimsizes(data&lon/),"integer")

rc = regcoef(ispan(0,29,1),data,tval,df)

df = df - 2 ; regcoef returns N, need N-2

beta_b = 0.5

signif = trend

signif = (/ 1-betainc(df/(df+tval^2),df/2.0,beta_b))*100. /)

Note that in the above coding regcoef sets the degrees of freedom to

simply the number of timesteps. If your data is highly auto-correlated,

then it is suggested that you use equiv_sample_size to compute the degrees

of freedom:

http://www.ncl.ucar.edu/Document/Functions/Built-in/equiv_sample_size.shtml

Hope that helps.. Adam

*> Dear all,
*

*> I would like to get the 30-year linear trend of my data (lat, lon, time)
*

*> in
*

*> each grid box. I have tried to write a script but didn't know whether it
*

*> is
*

*> correct. Could anyone help me to check with that?Also, I would like to
*

*> know
*

*> how can I test the significance of the trend using ncl. Thank in advance
*

*> for
*

*> any help given regarding on that.
*

*> Shinn
*

*>
*

*> trend = new((/dimsizes(data&lat), dimsizes(data&lon)/), "float")
*

*> trend!0 = "lat"
*

*> trend&lat = data&lat
*

*> trend!1 = "lon"
*

*> trend&lon = data&lon
*

*> trend&lat@units = "degrees north"
*

*> trend&lon@units = "degrees east"
*

*>
*

*> do i =0, dimsizes(data&lat)-1
*

*> do j = 0, dimsizes(data&lon)-1
*

*> temp = dtrend_msg(ispan(0,29,1),data(i,j,:),False,True)
*

*> trend(i,j) = (29*(temp@slope) + temp@y_intercept) ; is the trend
*

*> here
*

*> represent the 30-year linear trend?
*

*> end do
*

*> end do
*

*> _______________________________________________
*

*> ncl-talk mailing list
*

*> List instructions, subscriber options, unsubscribe:
*

*> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
*

*>
*

__________________________________________________

Adam Phillips asphilli@ucar.edu

NCAR/Climate and Global Dynamics Division

P.O. Box 3000 (303) 497-1726

Boulder, CO 80307-3000

http://www.cgd.ucar.edu/cas/asphilli

_______________________________________________

ncl-talk mailing list

List instructions, subscriber options, unsubscribe:

http://mailman.ucar.edu/mailman/listinfo/ncl-talk

Received on Tue Aug 2 16:52:06 2011

*
This archive was generated by hypermail 2.1.8
: Fri Aug 05 2011 - 14:53:55 MDT
*