[ncl-talk] Calculating a regression line & slope
Dennis Shea
shea at ucar.edu
Fri Aug 8 08:48:09 MDT 2014
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
More information about the ncl-talk
mailing list