[ncl-talk] Calculating a regression line & slope

[Dennis Shea]

This question is rather unclear. The following is a 'shot-in-the-dark'.
You have

        tasP = b->tas(0:364,:,:)

        TIMEP = cd_calendar(tasP&time,0)   ; TIMEP(ntim,6)

        tasCp = tasP - 273.15  ; Convert K to C
        tasCp = where(tasCp.gt.0.0, tasCp, 0.0)
        copy_VarMeta(tasP, tasCp)
        tasCp at units = "degC"
        printVarSummary(tasCp)

        pddP  = dim_sum_n_Wrap(tasCp, 0)       ; (lat,lon)
        pddP at long_name = "degree days for year="+toint(TIMEP(0,0))
        pddP at units = "degC"

; Calculate ablation (PDD * 5); add meta data
        Ab = pddP * 5
        copy_VarCoords(pddP, Ab)
        Ab at long_name = "Ablation"
        Ab at units           = "???"
        printVarSummary (Ab)

;---

        Orog = c->orog(:,:)                               ; 2D (lat,lon)

        Orog_1D = ndtooned(Orog)
        Ab_1D     = ndtooned(Ab)

; convenience, not required, sort into ascending orography

       iorog =* dim_pqsort_n*(Orog_1D,1,0)
       x = Orog_1D(iorog)    ; clarity
       y = Ab_1D(iorog)
       print(x(0:99)+"   "+y(0:99))    ; print 1st 100 values

   rc = regline(x,y)

or, if you have 6.2.0

      rc = regline_stats(x,y)  ; more information
      print(rc)

    df   = rc at nptxy-2
    prob = (1 - *betainc*
<http://www.ncl.ucar.edu/Document/Functions/Built-in/betainc.shtml>(df/(df+rc at tval^2),
df/2.0, 0.5) )

      yReg = rc*x + rc at yintercept   ; NCL array notation
                                    ; *yReg* is same length as *x*, *y*


      *print* <http://www.ncl.ucar.edu/Document/Functions/Built-in/print.shtml>("prob="+prob)

      *print* <http://www.ncl.ucar.edu/Document/Functions/Built-in/print.shtml>(x+"
 "+yReg)

http://www.ncl.ucar.edu/Applications/regress.shtml

http://www.ncl.ucar.edu/Document/Functions/Built-in/regline.shtml
http://www.ncl.ucar.edu/Document/Functions/Contributed/regline_stats.shtml

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


On Thu, Aug 7, 2014 at 12:26 PM, Lauren Jean Vargo <lvargo at unm.edu> wrote:

> Hello,
>
> I have a code that creates a scatter plot of 2 variables (ablation and
> surface altitude, both 2d arrays).
>
> Now I’m trying to 1) include a regression line on the plot, and 2)
> calculate the slope of that line (the regression coefficient).
> I’ve been trying to use regCoef since the arrays are both 2d, but the
> answers that are being printed are not reasonable.
>
> The files are both very large, but I can upload them if needed.
>
> I’m running ncl version 6.1.2, and the system is Darwin Kernel Version
> 13.3.0
>
> Here is the script:
>
> 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
>
> ; Read in temperature data & calculate PDD sum
>         b =
> addfile("tas_day_CCSM4_piControl_r2i1p1_09530101-09871231.nc","r")
>         tasP = b->tas(0:364,:,:)
>
>         TIMEP = cd_calendar(tasP&time,0)   ; TIMEP(ntim,6)
>
>         tasCp = tasP - 273.15  ; Convert K to C
>         tasCp = where(tasCp.gt.0.0, tasCp, 0.0)
>         pddP  = dim_sum_n_Wrap(tasCp, 0)       ; (lat,lon)
>         pddP at long_name = "degree days for year="+toint(TIMEP(0,0))
>         pddP at units = "degC"
>
> ; Calculate ablation (PDD * 5)
>         Ab = pddP * 5
>         printVarSummary (Ab)
>
> ; Read in topographic data
>         c = addfile("orog_fx_CCSM4_lgm_r0i0p0.nc","r")
>         Orog = c->orog(:,:)
>
> ; Calculate regression line
>         rc    = regCoef (Orog,Ab)
>         print (rc)
>
>
> ; Plot
>         res                                   = True                     ;
> plot mods desired
>         res at gsnMaximize           = True                    ; maximize
> plot in frame
>         res at xyMarkLineMode     = "Markers"           ; choose which have
> markers
>         res at xyMarker                  = 1                        ; choose
> type of marker
>         res at xyMarkerSizeF         = 0.005                 ; Marker size
> (default 0.01)
>
>         res at tiMainString        = "Ablation Gradient"  ; title
>
>   wks = gsn_open_wks("x11","CCSM_Ablation_gradient")
>   plot = gsn_csm_xy(wks,Orog(:,:),Ab(:,:),res)
>
>
> end
>
> Any help would be great!
> Thank you,
>
> Lauren Vargo
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20140808/7e3ee46a/attachment.html