[ncl-talk] Calculating a regression line & slope

[Lauren Jean Vargo]

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