[ncl-talk] Printing regression formula on plot

[Herb, Jason]

Hello all,

   I am attempting to print the linear regression equation mX+b on my XY scatter plots. Unfortunatly, I have not had much luck in getting it to plot. Code script is below. Any suggestions?

Thankyou


;This will create scatterplots of Modis derived PM2.5 data and surface measured PM2.5 data.

begin
f =  "QNC-XY.out"
data = asciiread(f,(/5856,2/),"float")
x = data(:,1)
; x = new((/2,366/),float)
        x at _FillValue=integertoshort(-999)
        x=where(x.lt.0,x at _FillValue,x)
x at long_name = "PSP PM25"
; x(0,:) = data(:,4)
        y = data(:,0)
y at long_name = "Modis Derived PM25"
        y at _FillValue=integertoshort(-999)
        y=where(y.lt.0,y at _FillValue,y)

;************************************************
; calculate the regression coefficient (slope)
;************************************************
   rcl    = regline(x, y)          ; slope
   rcl at units = "MODIS/PSP"
   print(rcl)

   rcl_anova = regline_stats(x,y) ; linear regression: ANOVA
   print(rcl_anova)

;************************************************
; create an array to hold both the original data
; and the calculated regression line
;     ---------
;     y = mx+b
;     m is the slope:       rc      returned from regline
;     b is the y intercept: rc at yave attribute of rc returned from regline
;************************************************
;   pltarry  := new ( (/2,366/), typeof(x), x at _FillValue)
   pltarry  := new ( (/2,5856/), typeof(x), x at _FillValue)
;************************************************************
   print(sprintf("%4.0f", x)+"  "+sprintf("%4.0f", y))

   pltarry(0,:) = y                                 ; use markers
;   pltarry(1,:) = rcl*x + rcl at yintercept              ; use solid line
   pltarry(1,:) = rcl*(x-rcl at xave) + rcl at yave ; y =m*x + B

;************************************************
; Plotting parameters
;************************************************
        plot = new (1, "graphic")
  wks  = gsn_open_wks("png","scatter_regress")             ; send graphics to PNG file

  resP                     = True                ; modify the panel plot
  resP at gsnMaximize         = True                ; maximize plot in frame
  resP at gsnPanelMainString  = "Modis vs PSP"  ; title
  resP at gsnPanelRowSpec     = True                   ; tell panel what order to plot
;  gsn_panel(wks,plot,(/2,1/),resP)


;  wks  = gsn_open_wks("png","scatter_regress")             ; send graphics to PNG file
  res = True
  res at gsnDraw = False
  res at gsnFrame    = False
  res at xyMarkLineModes     = (/"Markers","Lines"/)  ; choose which have markers
   res at xyMarkers           = 16                     ; choose type of marker
   res at xyMarkerColor       = "red"                  ; Marker color
   res at xyMarkerSizeF       = 0.005                  ; Marker size (default 0.01)
   res at xyDashPatterns      = 1                      ; solid line
   res at xyLineThicknesses   = (/1,2/)                ; set second line to 2
   res at tmYLFormat          = "f"                    ; not necessary but nicer labels

;   res at tiMainString        = "Modis vs PSP"  ; title
;    plt  = gsn_csm_xy (wks,x,pltarry,res)           ; create plot
;   plt  = gsn_csm_xy (wks,x,data,res)           ; create plot
    res  = True
  txres = True                                 ; label BW line
  txres at txFontHeightF = 0.0175                 ; font height
  txres at txJust        = "CenterCenter"         ; Set lable location

  if (rcl_anova at b(0).gt.0) then
      text  = "AOD_PM2.5 = "+ sprintf("%5.3f", rcl_anova at b(1))+"*PM25  + "+sprintf("%5.3f", abs(rcl_anova at b(0)))
  else
      text  = "AOD_PM2.5 = "+ sprintf("%5.3f", rcl_anova at b(1))+"*PM25  - "+sprintf("%5.3f", abs(rcl_anova at b(0)))
  end if
  plt  = gsn_csm_xy (wks,x,pltarry,res)           ; create plot
  gsn_add_text(wks,plt,text, 70, 5.6,txres)
  draw(plt)
  frame(wks)





;  res at gsnMaximize       = True                     ; maximize plot
;  res at tiMainString      = "Scatter Plot"           ; add title
;  res at xyMarkLineMode    = "Markers"                ; choose to use markers
;  res at xyMarkers         =  16                      ; choose type of marker
;  res at xyMarkerColor     = "NavyBlue"               ; Marker color
;  res at xyMarkerSizeF     = 0.01                     ; Marker size (default 0.01)
;
; plot  = gsn_xy (wks,x,y,res)        ; create plot
end


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20201023/a0ad4988/attachment.html>