Re: Segmentation fault with dtrend_msg

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Wed, 11 Oct 2006 17:25:02 -0600 (MDT)

There should be no need to convert the data or coordinates.

Perhaps, there is an error in the interface?

I have asked that the original file be sent to NCAR.

---
I am not sure why the plot code is so complicated.
Why is paneling used for one figure?
---
If "slope" is what ultimately is desired, I think
an  alternative might be regCoef:
http://www.ncl.ucar.edu/Document/Functions/Built-in/regCoef.shtml
The following does not detrend the data. It does return enough info
to test for significance.
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
 begin
  diri  = "/Users/bbaker/data/Vulnerability/"
  fili  = "co_tmax1895-1997.nc"
 
  a     = addfile(diri+fili, "r")
  ta1   = a->tmax                ; get data
  printVarSummary(ta1)           ; Note: NCL imports all meta data
  printMinMax(ta1, True)
  
  ta1   = ta1/100 
  ta1_at_units = ta1_at_units + "/100"        
  printMinMax(ta1, True)
  
  ntim  = filevardimsizes(f, "time") 
  yrs   = ispan(1,ntim,1)     
  rc    = regCoef(yrs, ta1(lat|:,lon|:,time|:) )
  
                      ; add meta data
  rc_at_long_name   = "regression coefficient"
  rc!0   = "lat"      ; name dimensions
  rc!1   = "lon"
  rc&lat = ta1&lat"   ; assign coordinate values to named dimensions
  rc&lon = ta1&lon"
  printVarSummary(rc)           ; variable overview
  printMinMax(rc, True)
                                ; statistics
  tval = onedtond(rc_at_tval , dimsizes(rc))
  df   = onedtond(rc_at_nptxy, dimsizes(rc)) - 2
  
  prob = betainc(df/(df+tval^2),df/2.0,0.5)       ; prob(nlat,nlon)
  prob_at_long_name = "probability"
  copy_VarCoords(rc, prob)
  printVarSummary(prob)
  printMinMax(prob, True)
>The following works to compute the trend and plot it (although, for
>some reason, the map outline is not showing up).
>I had more luck when I converted your double data to float to
>get rid of the segmentation fault.  If you want the slope, then
>use yDtrend_at_slope.
>Regarding the question you emailed just me, the following line:
>ta2 = new((/92,168,103/),float)
>does not work because your data is double and you are trying
>to fit it to float.  You can use dble2flt to convert double to float
>first.
>Michael
>
>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/shea_util.ncl"
>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
>
>begin
>
>a = addfile("co_tmax1895-1997.nc", "r")
>yrs = fspan(1.,103.,103)
>lat   = dble2flt(a->lat)                                  ; get lat
>lon   = dble2flt(a->lon)                                 ; get lon
>ta1   = a->tmax(:,:,:)/100.0                    ; get data
>ta1!0 = "time"
>ta1!1 = "lat"
>ta1!2 = "lon"
>ta1&time=yrs
>ta1&lat=lat
>ta1&lon=lon
>
>ta2 = dble2flt(ta1(lat|:,lon|:,time|:))
>yDtrend = dtrend_msg(yrs,ta2,False,True)     ; Calculate the trend
>slope = onedtond(yDtrend_at_slope, (/92,168/) )
>
>slope!0="lat"
>slope!1="lon"
>slope&lat=lat
>slope&lon=lon
>
>   print(min(slope))
>   print(max(slope))
>
>   wks   = gsn_open_wks ("ps", "script" )            ; open workstation
>   plot = new(1,graphic)
>   gsn_define_colormap(wks,"BlWhRe")
>   res                  = True                ; plot mods desired
>   res_at_gsnDraw         = False
>   res_at_gsnFrame        = False          ; don't advance frame
>   res_at_cnFillOn         = True               ; color Fill
>   res_at_cnLinesOn        =  False             ; Turn off contour lines
>   res_at_cnLevelSelectionMode = "ExplicitLevels" ; set manual contour  
>levels
>   res_at_cnLevels = (/-.02,-.015,-.01,-.005,.005,.01,.015,.02/)
>   res_at_cnFillColors = (/2,20,38,45,52,58,66,83,100/)
>   res_at_tiMainFontHeightF =  0.016
>   res_at_tiMainString         = "Trend"
>   res_at_tmYROn = False
>   res_at_tmXTOn = False
>   res_at_pmLabelBarDisplayMode = True
>   res_at_lbLabelBarOn = True
>   res_at_lbOrientation        = "vertical"
>   res_at_gsnAddCyclic = False
>   res_at_mpOutlineBoundarySets = "AllBoundaries"
>   res_at_mpGeophysicalLineThicknessF = 2.5
>   res_at_mpOutlineOn = True
>   res_at_cnFillDrawOrder      = "PreDraw"
>   res_at_cnLineDrawOrder = "PreDraw"
>   res_at_mpOceanFillColor     = "white"
>   res_at_mpLandFillColor = -1
>   res_at_mpMinLatF = min(lat)
>   res_at_mpMaxLatF = max(lat)
>   res_at_mpMinLonF = min(lon)
>   res_at_mpMaxLonF = max(lon)
>   res_at_mpCenterLonF = avg(lon)
>   plot(0) = gsn_csm_contour_map_ce(wks,slope,res)
>
>   pres             = True                 ; mod panel plot
>   pres_at_gsnDraw=True
>   pres_at_gsnFrame=True
>   pres_at_gsnMaximize = True                        ; blow up plot
>   pres_at_gsnPaperWidth=8.5
>   pres_at_gsnPaperHeight=11.0
>   pres_at_gsnPaperOrientation="portrait"
>   pres_at_gsnPaperMargin=0.04
>   pres_at_gsnPanelLeft = 0.02
>   pres_at_gsnPanelRight = 0.98
>   pres_at_gsnPanelBottom = 0.02
>   pres_at_gsnPanelTop = 0.96
>   pres_at_gsnPanelYWhiteSpacePercent = 3
>   pres_at_gsnPanelXWhiteSpacePercent = 3
>   pres_at_gsnPanelRowSpec=True
>   pres_at_txFontHeightF = .014
>   gsn_panel(wks,plot,(/1/),pres)
>
>end
>_______________________________________________
>ncl-talk mailing list
>ncl-talk_at_ucar.edu
>http://mailman.ucar.edu/mailman/listinfo/ncl-talk
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Oct 11 2006 - 17:25:02 MDT

This archive was generated by hypermail 2.2.0 : Thu Oct 12 2006 - 10:22:33 MDT