[ncl-talk] significant_trend_at_95_percentile

[Kunal Bali]

Dear NCL

I am  trying to get the trend significane plot at 95 percentile.
I used the script given below,
Could any one please let me know that is it correct way to do it.


  a = addfile("/media/Seagate Backup Plus Drive/work//test.nc","r")

   y  = a->TOTEXTTAU

; create x and calculate the regression coefficients (slopes, trends)
;************************************************

;   rc           = regCoef(y&time, y(lat|:,lon|:,time|:))*1000000
 ;  rc at long_name = "regression coefficient (trend)"
 ;  rc at units     = ts at units+"/day"
 ;  copy_VarCoords(y(0,:,:), rc)
 ;  printVarSummary(rc)
 ;  printMinMax(rc,1)

   y_clim=dim_avg_n_Wrap(y,0)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  dim_size   = dimsizes (y)
  aot_dtrend = dtrend_n(y,True,0)
  aot_trend  = onedtond(aot_dtrend at slope,(/dim_size(1),dim_size(2)/))

  copy_VarCoords(y_clim,aot_trend)
  aot_trend at _FillValue = y_clim at _FillValue

 printVarSummary(aot_trend)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;  Significance

  tval = new((/dim_size(1),dim_size(2)/),"float")      ; preallocate tval
as a float array and
  df = new((/dim_size(1),dim_size(2)/),"integer")     ; df as an integer
array for use in regcoef
  line_1 = ispan(0,dim_size(0)-1,1)
  rc = regcoef(line_1,y(lat|:,lon|:,time|:),tval,df)       ; regress z
against a straight line to
                                                                         ;
return the tval and degrees of freedom

  df = equiv_sample_size(y(lat|:,lon|:,time|:),0.05,0)
  df = df-2                                                           ;
regcoef/equiv_sample_size return N, need N-2

  beta_b = new((/dim_size(1),dim_size(2)/),"float")    ; preallocate space
for beta_b
  beta_b = 0.5                              ; set entire beta_b array to
0.5, the suggested value of beta_b
                                           ; according to betainc
documentation
  binc=betainc(df/(df+tval^2), df/2.0, beta_b)
  aot_trend_signif_percent = (1.-binc)*100.            ; significance of
trends ; expressed from 0 to 100%

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  aot_trend_sig = where (aot_trend_signif_percent.gt.95,aot_trend,aot_trend@
_FillValue)
  copy_VarCoords(y_clim,aot_trend_sig)

   plot= gsn_csm_contour_map_ce(wks,aot_trend_sig,res)


Thank You

Regards
Kunal Bali
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170113/f5d77b8d/attachment.html