/**********************************************************************/


data rawelev;
   input long lat elev;
   if (lat ge 35 and lat le 45)  /* Western United States        */
      and (long ge 100 and long le 120);
   x=long/57.295;                /* convert longitude to radians */
   y=lat/57.295;                 /* convert latitude to radians  */
   drop long lat;
   cards;
120.000 35 867.02
119.917 35 706.15
119.833 35 692.13
119.750 35 777.04
119.667 35 837.99
119.583 35 902.94
119.500 35 838.97
119.417 35 684.1
119.333 35 740.08
119.250 35 810.05
119.167 35 893.94
119.083 35 488.19
119.000 35 404.19
118.917 35 647.09
118.833 35 1209.86
118.750 35 1592.76
118.667 35 1495.92
118.583 35 1349.98
118.500 35 1161
   /* more data lines */
run;


/**********************************************************************/


proc g3grid data=rawelev out=gridelev;
   grid y*x=elev/join naxis1=250 naxis2=250;
run;


/**********************************************************************/


proc means data=rawelev noprint min max;
   var y x;
   output out=ranges min=latmin longmin
          max=latmax longmax;
run;


/**********************************************************************/


data _null_;
   set ranges;
   call symput('latmin',latmin);
   call symput('latmax',latmax);
   call symput('longmin',longmin);
   call symput('longmax',longmax);
run;


/**********************************************************************/


proc gproject data=maps.states out=map project=none
     latmin=&latmin latmax=&latmax
     longmax=&longmax longmin=&longmin;
   id state;
run;


/**********************************************************************/


data gridelev;
   set gridelev;
   x=-x;
run;


/**********************************************************************/


data anno;
   set map;
   by state notsorted segment notsorted;
   retain xsys ysys '2'
          when 'a'
          color 'black'
          size 5
          xfirst
          yfirst;

   x=-x;               /* multiply X by -1 to match contour data */
   if first.state or first.segment then do;
      function='move';     /* do a 'move' for first point        */
                           /* in each segment                    */
      xfirst=x;            /* save first point to close polygon  */
      yfirst=y;
   end;
   else function='draw';   /* do a 'draw' to intermediate points */
   output;

   if last.state or last.segment then do;
      x=xfirst;            /* connect last point in segment with */
                           /* first one                          */
      y=yfirst;
      function='draw';
      output;
   end;

run;


/**********************************************************************/


proc gcontour data=gridelev anno=anno;
   plot y*x=elev/levels=0 to 3500 by 500
                 pattern
                 join
                 noaxis
                 caxis=black
                 legend=legend1
                 haxis=axis1
                 vaxis=axis2;
   pattern1 color=vig    v=solid;
   pattern2 color=olive  v=solid;
   pattern3 color=vigy   v=solid;
   pattern4 color=daoy   v=solid;
   pattern5 color=vioy   v=solid;
   pattern6 color=vio    v=solid;
   pattern7 color=viro   v=solid;
   pattern8 color=maroon v=solid;

   axis1 length=80 pct width=5 origin=(,25 pct);
   axis2 length=60 pct;

   legend1 label=(position=(top center) f=centb
                 h=2 'Mean Elevation (meters)')
           value=(f=centb h=2 t=1 c=vig  t=2 c=olive
                              t=3 c=vigy t=4 c=daoy
                              t=5 c=vioy t=6 c=vio
                              t=7 c=viro t=8 c=maroon)
           shape=bar(8,2) mode=protect
           frame cshadow=black;

   title1 c=black font=centb h=4
          'Southern Rocky Mountain Region';
   footnote1 c=black f=centb h=2 j=r 'SAS/GRAPH'
           m=(+0,+.5) '02'x m=(+0,-.5) ' Software';
   run;
quit;


/**********************************************************************/