Help using Base SAS procedures

how to get the point of intersection on two plots?

Reply
Contributor
Posts: 36

how to get the point of intersection on two plots?

For example, I used proc gplot to do plots on y vs x and z vs x. How to get get the point of intersection on two plots? Thanks.

data mydata;

input x y z;

datalines;

1          550          100

2          500          200

3          450          300

4          400          400

5          350          500

6          300          600

7          250          700

8          200          800

;

run;

symbol1 v=none INTERPOL=join c=blue width=5;

symbol2 v=none INTERPOL=join c=red width=5;

proc gplot data=mydata;

plot (y z)*x/overlay;

run;

quit;

Respected Advisor
Posts: 3,777

how to get the point of intersection on two plots?

Your question seems somewhat vague.  Clearly the lines intersect at x=4 y,z=400.  Are you asking where the lines intersect if there is no common point in the data? 

You would need to calculate that point, perhaps someone can show us the SAS Procedure for this.  

Contributor
Posts: 36

how to get the point of intersection on two plots?

Thanks a lot for your reply. It's easy to get the intersection point in this example. But for some complicated plots, I need to find one way to  calculate the point where plots intersect.

PROC Star
Posts: 7,366

how to get the point of intersection on two plots?

I think they are addressing a similar problem in the following thread: http://www.listserv.uga.edu/cgi-bin/wa?A2=ind1002c&L=sas-l&D=1&O=A&J=0&P=31586

Hopefully, you can extract what you are looking for from at least one of the methods suggested.

Respected Advisor
Posts: 3,777

Re: how to get the point of intersection on two plots?

I hate to post this program because I "know" there has to be an easy way to do this using a SAS PROC.  Some of you statisticians should be able to help.

This programs plots two lines with the space between filled.  Finding the points where the two lines cross is needed to make the fill work properly.  The program uses JOIN interpolation, more complex interpolation (I don't know).  You will recognize the formula for the intersection of two lines.

I used an associative array so I could do "everything" in one data step.

The Y values are transformed to get one Y*X set for the filled part of the graph.

This does what you want (find the intersection) and then applies it.

goptions reset=all;

data test;

   input x0 yA yB; * had to change your names a bit;

   datalines;

1          550          100

2          500          575

3          450          300

4          400          400

5          350          500

6          300          600

7          850          700

8          200          800

;

run;

data _null_;

   if 0 then do;

      set test;

      by x0;

      end;

   length difx dif1 dif2 xL y1L y2L m1 m2 b1 b2 _m1 _m2 _b1 _b2 x y f1-f3 8;

   length _name_ $32;

   dcl hash d(ordered:'Y');

   dcl hiter di('d');

   do while(1);

      call vnext(_name_);

      put _name_=;

      if upcase(_name_) eq '_NAME_' then leave;

      if _name_ eq: 'FIRST.' then do;

         d.definekey(scan(_name_,-1,'.'));

         continue;

         end;

      else if _name_ eq: 'LAST.' then continue;

      d.definedata(_name_);

      end;

   d.definedone();

   * load data into hash, insures ordered by x0;

   do until(eof);

      set test end=eof;

      rc = d.add();

      end;

   * iterate over the hash, as if it were an array;

   do rc = di.first() by 0 while(rc eq 0);

      *put _all_;

      difx = dif(x0);

      dif1 = dif(yA);

      dif2 = dif(yB);

      xl   = lag(x0);

      y1l  = lag(yA);

      y2l  = lag(yB);

      * compute slope intercept for each pair of points;

      m1   = divide(dif1,difx);

      b1 = yA - (m1*x0);

      m2   = divide(dif2,difx);

      b2 = yB - (m2*x0);

     

      * find the intersection of each pair of lines;

      *m1(x) + b1 = m2(x) + b2;

      _m1=m1; _b1=b1;

      _m2=m2; _b2=b2;

      * factor m2;

      select(sign(m2));

         when(-1) do; _m1=_m1+abs(_m2); _m2=0; end;

         when(1do; _m1=_m1-abs(_m2); _m2=0; end;

         otherwise _m2=0;

         end;

      * factor b1;

      select(sign(b1));

         when(-1) do; _b2=_b2+abs(_b1); _b1=0; end;

         when(1do; _b2=_b2-abs(_b1); _b1=0; end;

         otherwise _b1=0;

         end;

      * Intersection(x,y) of lines;

      x = divide(_b2,_m1);

      y = m1*x + b1;

      * do they cross within the bounds of the points.;

      f1 = min(yA,y1L) lt y lt max(yA,y1L);

      f2 = min(yB,y2L) lt y lt max(yB,y2L);

      f3 = min(x0,xL)  lt x lt max(x0,xL);

      f = d.replace();

      rc = di.next();

      end;

   * Output all these calculations as a check;

   d.output(dataset:'intersection');

  

   length type order 8;

   dcl hash area(ordered:'Y');

   area.definekey('type','order');

   area.definedata('type','order','x0','value','yA','yB');

   area.definedone();

   do type = 1

      do rc = di.first() by 0 while(rc eq 0);

         order = x0;

         value = max(yA,yB);

         rcA = area.add();

         if sum(f1,f2,f3) gt 0 then do;

            if round(x,1e-6) ne round(xL,1e-6) and round(x,1e-6) ne round(x0,1e-6)then do;

               x0   = x;

               order = x;

               value = y;

               call missing(yA,yB);

               rcA = area.add();

               end;

            end;

         rc = di.next();

         end;

      end;     

   do type = 2

      do rc = di.last() by 0 while(rc eq 0);

         order = -x0;

         value = min(yA,yB);

         call missing(yA,yB);

         rcA = area.add();

         if sum(f1,f2,f3) gt 0 then do;

            if round(x,1e-6) ne round(xL,1e-6) and round(x,1e-6) ne round(x0,1e-6)then do;

               order = -x;

               value = y;

               x0   = x;

               rcA = area.add();

               end;

            end;

         rc = di.prev();

         end;

      end;     

   rcA = area.output(dataset:'AREA');

   stop;

   run;

proc print data=intersection;

   run;

proc contents data=area;

proc print data=area;

   run;

symbol interpol=msolid cv=lightgreen co=darkgreen;

symbol2 v=dot  color=black i=join w=2;

symbol3 v=star color=red i=join w=2;

proc gplot data=area;

   plot value*x0=1 yA*x0=2 yB*x0=3 / overlay;

   run;

   quit;

Message was edited by: data _null_

Respected Advisor
Posts: 4,662

Re: how to get the point of intersection on two plots?

The interpolation specification I=JOIN generates line segments, so there could be many intersection between the two curves. You can find all the intersention points with the program (I added an x value to cover a special case) :

data mydata;
input x y z;
datalines;
0          600          600
1          550          100
2          500          200
3          450          300
4          400          400
5          350          500
6          300          600
7          250          700
8          200          800
;
run;

proc sort data=myData; by x; run;

data meet(keep=x1 x2 xi yzi);
set myData;
x1=lag(x); x2=x;
y1=lag(y); y2=y;
z1=lag(z); z2=z;
if _n_ > 1 then do;
  alpha = (z2-y2) / (y1-y2-z1+z2);
  if alpha ge 0 and alpha lt 1 then do;
    xi = alpha*x1 + (1-alpha)*x2;
    yzi = alpha*y1 + (1-alpha)*y2;
    output;
  end;
end;
else if y=z then do;
  xi = x;
  yzi = y;
  output;
end;
run;

/* at this point, intersection points coordinates are in dataset meet */

proc sql;
create table myGraphData as
select A.*, B.xi, B.yzi
from myData as A left join meet as B on A.x = B.x2;

proc sgplot data=myGraphData;
series x=x y=y;
series x=x y=z;
scatter x=xi y=yzi / markerattrs=(symbol=CircleFilled size=8);
run;

PG

Extended by PG :

of if you want to use proc gplot :

proc sql;

create table myGData as

select * from

(select x, y, z, . as yzi from myData)

UNION

(select xi, ., ., yzi from meet)

order by yzi, x;

symbol1 v=none INTERPOL=join c=blue width=5;

symbol2 v=none INTERPOL=join c=red width=5;

symbol3 v=dot INTERPOL=none c=green width=5;

proc gplot data=myGData;

plot (y z yzi)*x/overlay;

run;

quit;

PG
Super User
Posts: 9,691

Re: how to get the point of intersection on two plots?

Just like PGStats . It is only suited for line.

data mydata;
input x y z;
datalines;
1          550          100
2          500          200
3          450          300
4          400          400
5          350          500
6          300          600
7          250          700
8          200          800
;
run;

data mean(keep=x mean_y);
 set mydata(obs=2) end=last;
 x=(z-y)*dif(x)/(dif(y)-dif(z))  + x ;
 mean_y=(z-y)*dif(y)/(dif(y)-dif(z))  + y ;
 if last then output;
run;
data want;
 set mydata mean;run;
symbol1 v=none INTERPOL=join c=blue width=5;
symbol2 v=none INTERPOL=join c=red width=5;
symbol3 v=dot INTERPOL=none c=green width=5;
proc gplot data=want;
plot (y z mean_y)*x/overlay;
run;
quit;

Ksharp

Ask a Question
Discussion stats
  • 6 replies
  • 1281 views
  • 0 likes
  • 5 in conversation