[ncl-talk] problem probability plotting

[sunmin park]

Dear NCL users

I try to plot probability>0.95% on a map
Used regcoef to get trend and betainc to get probability. Attached the
figure from the following script.

I do appreciate any help!

Best,
Sun-
************************************
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"

begin

diri =
"/Users/spark/03_Drought/Season_data/04_TOTALprecipitaton_season_mmperday_trend_21cen/"
fili = "PRE_MAMmean4Dmmperday_21CEN.nc"
filenc       = addfile (diri+fili+".nc", "r")
varname = "pre"
pme = filenc->$varname$
time = filenc->time

pmee = dim_avg_n_Wrap(pme,0)

rc = regCoef(time,pmee(lat|:,lon|:,time|:))
rc at long_name   = "regression coefficient"
rc at units = "mm/day/century"
copy_VarCoords(pmee(0,:,:),rc)

rcc = rc*100
copy_VarCoords(rc,rcc)

tval = onedtond(rc at tval , dimsizes(rc))
df   = onedtond(rc at nptxy-2, dimsizes(rc))
df at _FillValue=-9999
df   = where(df.le.0,df at _FillValue,df)

b=tval
b=0.5
prob = betainc(df/(df+tval^2),df/2.0,b)
prob at long_name = "probability"
copy_VarCoords(rcc,prob)

do i=0, dimsizes(lat)-1
  do j =0, dimsizes(lon)-1
    if(prob(i,j) .ltt. 0.05) then
      print("lat "+lat(i)+" lon "+lon(j))
    end if
  end do
end do

  wks = gsn_open_wks("eps","PRE_US_trend_MAM")
  gsn_define_colormap(wks,"BlRe") ;WhiteYellowOrangeRed") ;WhiteBlue")
 ; choose colormap

  res                     = True               ; plot mods desired
  res at gsnDraw             = False              ; don't draw yet
  res at gsnFrame            = False              ; don't advance frame yet
  res at gsnAddCyclic        = False

  res at cnFillOn            = True               ; turn on color
  res at gsnSpreadColors     = True               ; use full color map
  res at cnLinesOn           = False              ; no contour lines
  res at cnLineLabelsOn      = False              ; no line labels

  res at mpMaxLatF                   = 60     ; choose subregion
  res at mpMinLatF                   = 20
  res at mpMaxLonF                   = -75
  res at mpMinLonF                   = -150

  res at cnLevelSelectionMode = "ManualLevels"
  res at cnMinLevelValF       = -0.004                 ; min level
  res at cnMaxLevelValF       = 0.004     ; max level
  res at cnLevelSpacingF      = 0.0001      ; interval
  res at tmXBLabelFontHeightF     = 0.014         ; adjust some font heights
  res at tmYLLabelFontHeightF     = 0.014
  res at tiMainFontHeightF        = 0.022
  res at txFontHeightF            = 0.017

  res at lbLabelBarOn             = False
  res at mpOutlineBoundarySets = "GeophysicalAndUSStates"
  res at mpGeophysicalLineThicknessF = 1.0
  plot1 = gsn_csm_contour_map_ce(wks,rcc,res)

  res3 = True
  res3 at cnInfoLabelOn       = False
  res3 at cnLinesOn           = False
  res3 at cnLineLabelsOn      = False
  res3 at cnFillDotSizeF       = 0.004
  res3 at cnLevelSelectionMode = "ManualLevels"
  res3 at cnMinLevelValF       = 0.00
  res3 at cnMaxLevelValF       = 1.05
  res3 at cnLevelSpacingF      = 0.05
  llat = probt&lat
  llon = probt&lon
  plot2=gsn_add_polymarker(wks,plot1,llon,llat,res3)

  ppres                  = True
  ppres at gsnPanelLabelBar = True                   ; common label bar
  gsn_panel(wks,plot1,(/1,1/),ppres)

  delete(wks)
end
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190327/a1407b46/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: PRE_US_trend_MAM.eps
Type: application/postscript
Size: 602940 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190327/a1407b46/attachment-0001.eps>