Re: Regression/best fit line

Hi Melissa,
According to the regline documentation the two input timeseries must be
the same length, and you are trying to regress an array of size # of
latitude points between 0:30N and an array of size # of longitude points
between 0:55E. Those two arrays will likely be different in size.

To create a best fit line, you would typically want to regress a 1D
array against a straight line, as is shown in these two examples:
http://www.ncl.ucar.edu/Applications/scatter.shtml#ex4 (example #4)
http://www.ncl.ucar.edu/Applications/regress.shtml#ex1 (example #1)
I am not familiar with the SICZ, so I cannot offer any further advice on
how to calculate/identify this particular index/phenomenon.

If the above does not help please respond back to the ncl-talk email
list with a more detailed explanation about the process needed to
calculate the SICZ.
Adam

On 09/19/2013 04:45 AM, Melissa Lazenby wrote:
> Hi
>
> I am trying to draw a best fit line or regression line for maximum
> precipitation to find the SICZ and I am having trouble with drawing
> the best fit line. The regline function comes up with the error the
> input arrays must be the same length but i want to draw a reg line for
> the latitudes 0 to 30S and longitudes 0 to 55E, is there any other way
> to approach this.
> Below is my code.
>
> ;*************************************************
> ; reglineoverlay.ncl
> ;
> ;
> ;*************************************************
> 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"
>
> begin
> ;************************************************
> ; Read in Precip Data.
> ;************************************************
> fili =
> "/mnt/nfs2/geog/ml382/melphd/regressionline/siczoutputprregline/pr_Amon_ACCESS1-0_historical_safrica_climDJF1.nc"
> ; data
>
> f = addfile (fili , "r") ; add file
> lat = f->lat ; get lat
> lon = f->lon ; get lon
> time = f->time ; get time
> level = f->z ; get level
> pr = f->pr
> pr1 = pr(:,0,:,:) ; get precip
> pr2 = pr(0,0,:,:) ; ignore
> level&time, just want lat, lon
>
> printVarSummary(pr2)
>
> ;************************************************
> ; plotting parameters
> ;************************************************
> wks = gsn_open_wks("X11","ACCESS1-0_contour") ;
> specifies a plot
> gsn_define_colormap(wks,"gsltod") ; choose color map
>
> res = True ; plot mods desired
> res@gsnAddCyclic = False
> res@cnFillOn = True ; color on
> res@lbLabelStride = 2 ; every other label
> res@gsnSpreadColors = True ; use full range of
> color map
> res@gsnFrame = True ; don't advance frame yet
> res@cnLineLabelsOn = False ; no contour line labels
>
>
> res@tiMainString = "ACCESS1-0" ; title
>
> res@vpXF = 0.12 ; default is 0.2 change aspect
> res@vpYF = 0.8 ; default is 0.8 ration
> res@vpHeightF = 0.4 ; default is 0.6
> res@vpWidthF = 0.8 ; default is 0.6
>
> res@cnLevelSelectionMode = "ManualLevels" ; manual levels
> res@cnMinLevelValF = 0 ; min level
> res@cnMaxLevelValF = 400 ; max level
> res@cnLevelSpacingF = 50 ; contour spacing
>
> res@tmXBMode = "Explicit" ; Define own tick mark labels.
> res@tmXBValues = (/ 0.,10.,20.,30.,40.,50. /)
> res@tmXBLabels = (/"0","10E","20E","30E","40E","50E" /)
>
> res@tmYLMode = "Explicit" ; Define own tick mark labels.
> res@tmYLValues = (/ -29.,-20.,-10.,-1. /)
> res@tmYLLabels = (/"30S","20S","10S","0" /)
>
> ; reverse the first two colors
> setvalues wks
> "wkColorMap" : "gsltod"
> "wkForegroundColor" : (/0.,0.,0./)
> "wkBackgroundColor" : (/1.,1.,1./)
> end setvalues
> plot = gsn_csm_contour(wks,pr2({-30:0},{0:55}),res) ; create plot
>
> ;***************************************************
> ;Overlay plot of the regression line
> ;***************************************************
>
> ;pr1 = pr1(time, lat, lon)
>
> xMaxLon = dim_max( pr1(time|:, lat|:, lon|:) ) ; ==>
> xMaxLon(ntim,nlat)
>
> printVarSummary(xMaxLon)
>
> ;xMaxLon = pr1(time,lat)
> xMaxLatValues = xMaxLon(0,:)
>
> printVarSummary(xMaxLatValues)
> print(xMaxLatValues)
>
> pr1@_LatValues = f->lat
>
> lat=where(xMaxLon(0,:).eq. xMaxLatValues, lat, pr1@_LatValues)
>
> print(lat)
>
>
> rc = regline (lon,lat)
>
>
> print(rc)
>
>
> ; Note use of attributes
> df = rc@nptxy-2
> prob = (1 - betainc(df/(df+rc@tval^2), df/2.0, 0.5) )
>
> yReg = rc*x + rc@yintercept ; NCL array notation
> ; yReg is same length as x, y
> print(rc)
>
> print("prob="+prob)
>
> print(yReg)
>
> ;************************************************
> ; Create an array to hold both the original data
> ; and the calculated regression line.
> ;************************************************
> data = new ( (/2,dimsizes(pr3)/), typeof(pr3))
> ;data(0,:) = pr3
> ; y = mx+b
> ; m is the slope: rc returned from regline
> ; b is the y intercept: rc@yave attribute of rc returned from regline
> ;data(1,:) = rc*(x-rc@xave) + rc@yave
> ;************************************************
> ; plotting parameters
> ;************************************************
> wks = gsn_open_wks("X11","ACCESS1-0_regline") ;
> specifies a ps plot
>
> res = True ; plot mods desired
> res@gsnMaximize = True ; maximize plot in frame
> res@xyMarkLineModes = (/"Markers","Lines"/) ; choose which have
> markers
> res@xyMarkers = 16 ; choose type of marker
> res@xyMarkerColor = "red" ; Marker color
> res@xyMarkerSizeF = 0.01 ; Marker size
> (default 0.01)
> res@xyDashPatterns = 1 ; solid line
> res@xyLineThicknesses = (/1,2/) ; set second line to 2
>
> res@tiMainString = "Regline" ; title
>
> plot2 = gsn_csm_xy (wks,pr3,data,res) ; create plot
> end
>
>
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

-- 
______________________________________________________________
Adam Phillips                                asphilli@ucar.edu
NCAR/Climate and Global Dynamics Division       (303) 497-1726
P.O. Box 3000				
Boulder, CO 80307-3000    http://www.cgd.ucar.edu/cas/asphilli

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Sep 20 15:55:36 2013