missing values in rtest

From: Will Hobbs <whobbs_at_nyahnyahspammersnyahnyah>
Date: Mon, 16 Oct 2006 15:41:55 PDT

('binary' encoding is not supported, stored as-is) I am having problems getting rtest to handle missing/zero values. The
script below calculates the acf of an SST dataset at each gridpoint,
calculates the significance of that acf, and plots an array of the max
lag at which the acf is still significant. Inevitably, because of land
surface thare are a large number of missing values.

When I run the script, i get the following error message,

 ==>SLATEC/NCL: DLBETA: BOTH ARGUMENTS MUST BE GT ZERO: NERR= 1: LEVEL= 2
 ==>SLATEC/NCL: DGAMMA: X IS 0: NERR= 4: LEVEL= 2
 ==>SLATEC/NCL: DBETAI: P AND/OR Q IS LE ZERO: NERR= 2: LEVEL= 2

I have tried forcing zero values of nval to 1, and forcing zero values of
acf to 0.0001, but this doesn't seem to help. What am I doing wrong?

Will

;=====================================================================;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;======================================================================
;checks acf of ocean data, useful for estimating p for an AR(p) process
                                                                                
                                                                                
begin
                                                                                
;************************************************
; Read in test data
;*********************************************
                                                                                
  diri0 = "~/datasets/ocean/"
  FIL0 = "hadisst.187001-200412.nc"
                                                                                
  f0 = addfile(diri0+FIL0, "r")
                                                                          
x = f0->sst(1320:,{-90:0},:)
                                                                                
  lat = f0->lat({-90:0})
  lon = f0->lon
                                                                                
;remove annual cycle
                                                                                
  y = rmMonAnnCycLLT(x(lat|:,lon|:,time|:))
  delete(x)
                                                                                
  dimz = dimsizes(y)
  ntim = dimz(2)
  nlat = dimz(0)
  nlon = dimz(1)
                                                                                
                                                                                
;****************************************************
; Calc acf/sig
;************************************************
                                                                                
  mxlag = 12
  acf = esacr(y, mxlag)
                                                                                
  nval = new((/nlat, nlon/), integer)
                                                                                
  do lt = 0, nlat-1
    do ln = 0, nlon-1
                                                                                
      nval(lt, ln) = num(.not.ismissing(y(lt,ln,:)))
                                                                                
    end do
  end do
                                                                                
  nval1D = ndtooned(nval)
  nval1D(ind(nval1D.eq.0)) = 1
  nval = onedtond(nval1D, dimsizes(nval))
                                                                                
  acf1D = ndtooned(acf)
  acf1D(ind(acf1D.eq.0)) = 0.0001
  acf = onedtond(acf1D, dimsizes(acf))
                                                                                
  siglvl = 0.05
                                                                                
  lags1D = new((nlat*nlon), integer)
                                                                                
  do ml = 1, mxlag
                                                                                
    prob = rtest(acf(:,:,ml), nval, 0) ;prob as dims (nlat, nlon)
    prob1D = ndtooned(prob)
                                                                                
    sigs = ind(prob1D.le.siglvl) ;locn where acf is significant
                                                                                
    lags1D(sigs) = mxlag
                                                                                
    delete(sigs)
                                                                                
  end do
                                                                                
  delete(prob)
  delete(prob1D)
  delete(nval)
                                                                                
  lags = y(:,:,0) ;get metadata
  delete(y)
                                                                                
  lags = (/ onedtond(lags1D, (/nlat, nlon/)) /)
                                                                                
                                                                                
;*******************************************************
;plot lags
;*******************************************************
                                                                                
  wks = gsn_open_wks("ps", "sst_acf_test")
                                                                                
  res = True
  res_at_gsnPolar = "SH"
  res_at_mpMaxLatF = -20
  res_at_mpGridMaxLatF = -90
  res_at_mpGridLatSpacingF = 30.
  res_at_mpGridLonSpacingF = 45.
  res_at_mpFillOn = True
  res_at_mpOutlineOn = True
                                                                                
  res_at_cnLevelSelectionMode = "ExplicitLevels"
  res_at_cnLevels = ispan(1,12,1)
  res_at_cnLinesOn = True
  res_at_cnFillOn = False
  res_at_lbLabelBarOn = False
  res_at_cnLineLabelsOn = True
                                                                                
  plot = gsn_csm_contour_map_polar(wks, lags, res)
                                                                                
end
                                                                                

============================
Room 3221C
UCLA Department of Geography
1255 Bunche Hall
Los Angeles CA 90095-1524

_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Oct 16 2006 - 16:41:55 MDT