<div dir="ltr"><div><div>This question is rather unclear. The following is a 'shot-in-the-dark'.<br></div><div>You have<br></div><div><br> tasP = b->tas(0:364,:,:)<br>
<br>
TIMEP = cd_calendar(tasP&time,0) ; TIMEP(ntim,6)<br>
<br>
tasCp = tasP - 273.15 ; Convert K to C<br>
tasCp = where(tasCp.gt.0.0, tasCp, 0.0)<br></div> copy_VarMeta(tasP, tasCp)<br> tasCp@units = "degC"<br></div> printVarSummary(tasCp)<br><div><br><div>
pddP = dim_sum_n_Wrap(tasCp, 0) ; (lat,lon)<br>
pddP@long_name = "degree days for year="+toint(TIMEP(0,0))<br>
pddP@units = "degC"<br>
<br>
; Calculate ablation (PDD * 5); add meta data<br>
Ab = pddP * 5<br></div><div> copy_VarCoords(pddP, Ab)<br></div><div> Ab@long_name = "Ablation"<br></div><div> Ab@units = "???"<br></div><div>
printVarSummary (Ab) <br><br>;---<br>
<br><div> Orog = c->orog(:,:) ; 2D (lat,lon)<br><br></div><div> Orog_1D = ndtooned(Orog)<br></div><div> Ab_1D = ndtooned(Ab)<br><br></div><div>; convenience, not required, sort into ascending orography<br>
<br> <strong> </strong> iorog =<strong> dim_pqsort_n</strong>(Orog_1D,1,0)<br> x = Orog_1D(iorog) ; clarity<br> y = Ab_1D(iorog)<br></div><div> print(x(0:99)+" "+y(0:99)) ; print 1st 100 values<pre>
rc = regline(x,y)<br></pre></div></div><div>or, if you have 6.2.0 <br><br></div><div> rc = regline_stats(x,y) ; more information<br></div><div> print(rc)<br></div><div><br></div><div><div><pre> df = rc@nptxy-2
prob = (1 - <a href="http://www.ncl.ucar.edu/Document/Functions/Built-in/betainc.shtml"><strong>betainc</strong></a>(df/(df+rc@tval^2), df/2.0, 0.5) )<br><br> yReg = rc*x + rc@yintercept ; NCL array notation
; <em>yReg</em> is same length as <em>x</em>, <em>y</em>
<a href="http://www.ncl.ucar.edu/Document/Functions/Built-in/print.shtml"><strong>print</strong></a>("prob="+prob)
<a href="http://www.ncl.ucar.edu/Document/Functions/Built-in/print.shtml"><strong>print</strong></a>(x+" "+yReg)</pre><a href="http://www.ncl.ucar.edu/Applications/regress.shtml">http://www.ncl.ucar.edu/Applications/regress.shtml</a> <br>
<br></div><div><a href="http://www.ncl.ucar.edu/Document/Functions/Built-in/regline.shtml">http://www.ncl.ucar.edu/Document/Functions/Built-in/regline.shtml</a> <br><a href="http://www.ncl.ucar.edu/Document/Functions/Contributed/regline_stats.shtml">http://www.ncl.ucar.edu/Document/Functions/Contributed/regline_stats.shtml</a> <br>
<a href="http://www.ncl.ucar.edu/Document/Functions/Built-in/dim_pqsort_n.shtml">http://www.ncl.ucar.edu/Document/Functions/Built-in/dim_pqsort_n.shtml</a><br></div></div></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">
On Thu, Aug 7, 2014 at 12:26 PM, Lauren Jean Vargo <span dir="ltr"><<a href="mailto:lvargo@unm.edu" target="_blank">lvargo@unm.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hello,<br>
<br>
I have a code that creates a scatter plot of 2 variables (ablation and surface altitude, both 2d arrays).<br>
<br>
Now I’m trying to 1) include a regression line on the plot, and 2) calculate the slope of that line (the regression coefficient).<br>
I’ve been trying to use regCoef since the arrays are both 2d, but the answers that are being printed are not reasonable.<br>
<br>
The files are both very large, but I can upload them if needed.<br>
<br>
I’m running ncl version 6.1.2, and the system is Darwin Kernel Version 13.3.0<br>
<br>
Here is the script:<br>
<br>
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"<br>
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"<br>
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"<br>
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"<br>
<br>
begin<br>
<br>
; Read in temperature data & calculate PDD sum<br>
b = addfile("tas_day_CCSM4_piControl_r2i1p1_09530101-09871231.nc","r")<br>
tasP = b->tas(0:364,:,:)<br>
<br>
TIMEP = cd_calendar(tasP&time,0) ; TIMEP(ntim,6)<br>
<br>
tasCp = tasP - 273.15 ; Convert K to C<br>
tasCp = where(tasCp.gt.0.0, tasCp, 0.0)<br>
pddP = dim_sum_n_Wrap(tasCp, 0) ; (lat,lon)<br>
pddP@long_name = "degree days for year="+toint(TIMEP(0,0))<br>
pddP@units = "degC"<br>
<br>
; Calculate ablation (PDD * 5)<br>
Ab = pddP * 5<br>
printVarSummary (Ab)<br>
<br>
; Read in topographic data<br>
c = addfile("orog_fx_CCSM4_lgm_r0i0p0.nc","r")<br>
Orog = c->orog(:,:)<br>
<br>
; Calculate regression line<br>
rc = regCoef (Orog,Ab)<br>
print (rc)<br>
<br>
<br>
; Plot<br>
res = True ; plot mods desired<br>
res@gsnMaximize = True ; maximize plot in frame<br>
res@xyMarkLineMode = "Markers" ; choose which have markers<br>
res@xyMarker = 1 ; choose type of marker<br>
res@xyMarkerSizeF = 0.005 ; Marker size (default 0.01)<br>
<br>
res@tiMainString = "Ablation Gradient" ; title<br>
<br>
wks = gsn_open_wks("x11","CCSM_Ablation_gradient")<br>
plot = gsn_csm_xy(wks,Orog(:,:),Ab(:,:),res)<br>
<br>
<br>
end<br>
<br>
Any help would be great!<br>
Thank you,<br>
<br>
Lauren Vargo<br>
<br>
<br>
_______________________________________________<br>
ncl-talk mailing list<br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" target="_blank">http://mailman.ucar.edu/mailman/listinfo/ncl-talk</a><br>
</blockquote></div><br></div>