• Print |
• E-mail  |
• Bookmark

/* Create a data set name CITIES to contain the X and Y */
/* longitude/latitude values in radians for the cities  */

data cities;

  /* Get the coordinates from the MAPS.USCITY data set */
  set maps.uscity end=last;
  keep x y city long lat flag;	

  /* Subset specific cities */
  where capital='Y' and state in (48 30 23 12);

  /* Count the number of cities in the data set */
  count+1;

  /* Convert degrees to radians */
  x=atan(1)/45 * long;
  y=atan(1)/45 * lat;
  
  flag=1;

  /* Create a macro variable for the number of cities */
  if last then call symput('num',count);
run;

/* Combine the cities data set with an unprojected map data set */

data combo;

  /* Remove Alaska, Hawaii, and Puerto Rico */
  set maps.states(where=(state not in(2 15 72))) cities;
run;

/* Project the combined data set */

proc gproject data=combo out=states dupok;
  id state;
run;
quit;

/* Split the projected data set up into a map data set, ST */
/* and a cites data set, CTY using the flag variable       */

data st cty;

  /* Specify the projected data set. */
  set states;

  /* If observation is from the cities data set */
  if flag=1 then output cty;

  /* Otherwise observation is from the map data set */
  else output st;
run;

/* Create an annotate data set to label the cities */

data ctylbl;
  length function color style $8 text $25;
  retain xsys ysys '2' when 'a';

  /* Use the projected cities data set */
  set cty;

  /* Create an observation to place a symbol at the city */
  function='label';
  style='special';
  text='M';
  color='black';
  position='5';
  size=1.25;
  output;

  /* Create an observation to place a label at the city */
  style='swissb';
  text=city;
  position='E';
  output;
run;

/* Create a data set that contains the distance from        */
/* each city in the data set to every other city            */
/* use the unprojected values for calculating the distances */
/* keep the projected values in the data set for annotating */

data citydist;
  keep fromcity tocity dist frprojx frprojy toprojx toprojy;
  set cty;

  /* Get the values for the city */
  /* Convert degrees to radians  */
  fromx=atan(1)/45 * long;
  fromy=atan(1)/45 * lat;
  fromcity=city;
  
  /* Get the projected values for use in annotating */
  frprojx=x;
  frprojy=y;

  /* Get the observations for each of the cities */
  do i=1 to #
    set cty point=i;

	/* Use the unprojected values for the destination city */
    tox=atan(1)/45 * long;
    toy=atan(1)/45 * lat;
    tocity=city;

	/* Keep the projected values in the data set for annotation */
	toprojx=x;
	toprojy=y;

    /* If fromcity is the tocity, delete the observation */
    if tocity = fromcity then delete;

    /* Calculate the distance between the cities */
    /* using the Great Circle Distance Formula   */
    else
    Dist = round(3949.99 * arcos( sin( fromy ) * sin( toy ) +
            cos( fromy ) * cos( toy ) *
            cos( fromx - tox ) ));
    output;
  end;
run;

/* Create an annotate data set to draw lines between the cities */
/* and label the distance                                       */

data anno;
  length function color style $8 text $25;
  retain xsys ysys '2' when 'a';
  keep function color size style text xsys ysys when x y position angle;

  set citydist;

  /* Move to the first city */ 
  x=frprojx;
  y=frprojy;
  function='move';
  output;

 /* Create an observation to draw from one city to the next */
  x=toprojx;  
  y=toprojy;
  function='draw';
  color='red';
  size=1.75;
  output;

  /* Determine the angle of the line */
  deg=atan(abs(frprojy-toprojy)/abs(frprojx-toprojx));

  /* Convert radians to degrees for the ANGLE variable */
  deg=45/atan(1)*deg;

  /* Change the angle for a downward sloping line */
  if (frprojx>toprojx and frprojy<toprojy) or (toprojx>frprojx and toprojy<frprojy) then deg=-deg;

  /* Use the midpoint formula to determine the center point of the line */
  x=(frprojx + toprojx)/2;
  y=(frprojy + toprojy)/2;

  /* Create an observation to label the distance. */
  /* Angle the text to match the line             */ 
  function='label';
  style='swissb';
  color='black';
  text=dist||' Miles';
  position='E';
  angle=round(deg);
  size=1;
  output;
run;

/* Define the patterns for the map areas */

pattern1 v=me c=green r=50;

title 'Distances Between Selected Cities';

/* Generate a map with cities and the distance between them */

proc gmap data=st map=st anno=anno;
  id state;
  choro state / discrete anno=ctylbl coutline=graycc nolegend;
run;
quit;

Operating System and Release Information