Re: Regline and Locmax

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Mon Sep 23 2013 - 12:08:15 MDT

I think people are confused.

If I understand the question correctly, you are asking
"Why is the regression coefficient positive?"

In your code, you want to regress y (lat of max values)
onto x (lon of max valuess). Yet prior to the
call of regline you do

    x = ispan(0,dimsizes(y)-1,1)*1.

Not sure why you are doing that.

---
Given the figure you sent .... Do you want the following?
Your x (lon) are in 'random order'. I think you want to sort the
x in ascending order. See the following.
        ; your data
     x  = (/ 43.5, 37.5, 37.5, 51, 46.5,30, 22.5, 52.5, 28.5/)
     y  = (/-26.25, -23.25, -20.25, -17.2, -14.25, -12.75 \
           ,-11.25, -11.25, -3.75/)
     ii = dim_pqsort_n(x,1,0)            ; indices of ascending x
     print(ii+"   "+x(ii)+"   "+y(ii))
     print("========================")
     rc = regline(x(ii),y(ii))
     print(rc)
The result is
(0)	6   22.5   -11.25
(1)	8   28.5   -3.75
(2)	5   30   -12.75
(3)	2   37.5   -20.25
(4)	1   37.5   -23.25
(5)	0   43.5   -26.25
(6)	4   46.5   -14.25
(7)	3   51   -17.2
(8)	7   52.5   -11.25
(0)	========================
Variable: rc
Type: float
Total Size: 4 bytes
             1 values
Number of Dimensions: 1
Dimensions and sizes:	[1]
Coordinates:
Number Of Attributes: 7
   _FillValue :	9.96921e+36
   yintercept :	-6.582467
   yave :	-15.57778
   xave :	38.83333
   nptxy :	9
   rstd :	0.2353886
   tval :	-0.9840701
(0)	-0.2316389
On 9/23/13 4:40 AM, Melissa Lazenby wrote:
> Hi All
>
> I am trying to identify a specific phenomena, the South Indian Convergence Zone (SICZ) and I have decided to identify it by picking out the points of maximum precipitation for DJF over southern Africa over a range of longitudes from 0 to 55E.
>
> I first want to find the maximum precipitation values and there corresponding lats and lons so i have used the local_max function on the dataset, which results on 9 values for each lon, lat and max value.
>
> Then what I want to do is draw  a regression line through the max precip values and am having trouble doing that.
> The SICZ slope should be negative but is coming out positive and I am not sure what I am doing wrong.
>
> I have attached an image of what I would potentially like to produce but over southern Africa.
> I have attached my code and the outputs below.
>
> Code:
> ;*************************************************
> ; regline.ncl
> ;
>
> ;*************************************************
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.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)
>
>
>    locmax = local_max(pr2,False,0.)
>     x=lon(locmax@xi)                ; get lat/lon points of maxima
>     y=lat(locmax@yi)
>     z=locmax@maxval
>
>    print(locmax)
>    print(y)
>    print(x)
>    print(z)
>
>    x = ispan(0,dimsizes(y)-1,1)*1.
>    rc = regline(x,y)
>
>    print(rc)                                  ; Note use of attributes
>
>
> ;************************************************
> ; Create an array to hold both the original data
> ; and the calculated regression line.
> ;************************************************
>   data      = new ( (/2,dimsizes(y)/), typeof(y))
>   data(0,:) = y
> ; 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@gsnAddCyclic        = False
>   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,pr2&lon,data,res)        ; create plot
> end
>
>
>
> Output:
> Variable: pr2
> Type: float
> Total Size: 2960 bytes
>              740 values
> Number of Dimensions: 2
> Dimensions and sizes:    [lat | 20] x [lon | 37]
> Coordinates:
>              lat: [-29.25..-0.75]
>              lon: [   0..  54]
> Number Of Attributes: 11
>    z :     0
>    valid_max :    0.001398626
>    valid_min :    -6.214417e-24
>    time :    11
>    date :    01/01/50
>    title :    pr
>    name :    pr
>    source :
>    _FillValue :    2e+20
>    units :
>    long_name :    pr
>
>
> Variable: locmax
> Type: integer
> Total Size: 4 bytes
>              1 values
> Number of Dimensions: 1
> Dimensions and sizes:    [1]
> Coordinates:
> Number Of Attributes: 3
>    maxval :    ( 104.526, 231.4876, 239.5204, 300.4584, 541.7598, 344.9508, 349.5831, 353.6486, 298.2738 )
>    xi :    ( 29, 25, 25, 34, 31, 20, 15, 35, 19 )
>    yi :    ( 2, 4, 6, 8, 10, 11, 12, 12, 17 )
> (0)    9
>
>
> Variable: y
> Type: double
> Total Size: 72 bytes
>              9 values
> Number of Dimensions: 1
> Dimensions and sizes:    [lat | 9]
> Coordinates:
>              lat: [-26.25..-3.75]
> Number Of Attributes: 4
>    axis :    Y
>    standard_name :    latitude
>    units :    degrees_north
>    long_name :    latitude
> (0)    -26.25
> (1)    -23.25
> (2)    -20.25
> (3)    -17.25
> (4)    -14.25
> (5)    -12.75
> (6)    -11.25
> (7)    -11.25
> (8)    -3.75
>
>
> Variable: x
> Type: double
> Total Size: 72 bytes
>              9 values
> Number of Dimensions: 1
> Dimensions and sizes:    [lon | 9]
> Coordinates:
>              lon: [43.5..28.5]
> Number Of Attributes: 4
>    axis :    X
>    standard_name :    longitude
>    units :    degrees_east
>    long_name :    longitude
> (0)    43.5
> (1)    37.5
> (2)    37.5
> (3)      51
> (4)    46.5
> (5)      30
> (6)    22.5
> (7)    52.5
> (8)    28.5
>
>
> Variable: z
> Type: float
> Total Size: 36 bytes
>              9 values
> Number of Dimensions: 1
> Dimensions and sizes:    [9]
> Coordinates:
> (0)    104.526
> (1)    231.4876
> (2)    239.5204
> (3)    300.4584
> (4)    541.7598
> (5)    344.9508
> (6)    349.5831
> (7)    353.6486
> (8)    298.2738
>
>
> Variable: rc
> Type: double
> Total Size: 8 bytes
>              1 values
> Number of Dimensions: 1
> Dimensions and sizes:    [1]
> Coordinates:
> Number Of Attributes: 7
>    _FillValue :    9.969209968386869e+36
>    yintercept :    -25.48333333333333
>    yave :    -15.58333333333333
>    xave :       4
>    nptxy :    9
>    rstd :    0.2009649341599293
>    tval :    12.31558137416341
> (0)    2.475
> (0)    gsn_csm_xy: Fatal: X and Y must have the same dimensions sizes, or one must be one-dimensional and both have the same rightmost dimension.
>
> Many thanks for any help in this regard!
>
> Kind Regards
> Melissa
>
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Sep 23 12:08:25 2013

This archive was generated by hypermail 2.1.8 : Tue Oct 01 2013 - 14:41:43 MDT