Regression/best fit line

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
Received on Thu Sep 19 04:47:13 2013