• Print |
• E-mail  |
• Bookmark

/****************************************************************/
/* This program calculates the visual centers, or centers of    */
/* gravity, for counties within a state.  The county names      */
/* are annotated onto the centers of the counties.              */
/* the LIBNAME "MAPS" must point to the SAS data                */
/*          library containing the maps data sets.              */
/*                                                              */
/****************************************************************/


/* The program consists of the following five steps:            */

/* 1.  The first data step creates two data sets: MAP and       */
/*     POINTS. The county points for the selected state are     */
/*     subsetted out of the SAS data set MAPS.COUNTIES (or      */
/*     MAPS.COUNTY for PC SAS/GRAPH), and output to the MAP     */
/*     data set.  One observation for each county, containing   */
/*     the last point and the total number of points used to    */
/*     create the county, is output to the POINTS data set.     */

goptions reset=all border ctext=black colors=(black);

data counties;
   set maps.counties;
   if density < 6 ;
run;

data map(drop=npoints) points(keep=x y npoints
                              rename=(x=xlast y=ylast));
   set counties;
      by state county notsorted;  
   if state=stfips('NV');         /* Change to whatever state. */
   retain flag 'M';               
   output map;
   npoints + 1;
   if last.county then do;
      output points;
      npoints=0;
   end;
run;
 
 /* 2.  Add county names to the MAP data set to create the  */
 /*     NEWMAP data set.                                    */

data newmap;
   merge map(in=inmap) maps.cntyname;
      by state county;
   if inmap;
run;
 
 /* 3.  This data step creates an ANNOTATE= data set called      */
 /*     CENTERS.  The centers of gravity, XCG and YCG, are cal-  */
 /*     culated for each county.  The number of observations in  */
 /*     the POINTS data set corresponds to the number of         */
 /*     counties.  The number of observations in the NEWMAP data */
 /*     set corresponds to the number of points in all the       */
 /*     counties.  For each observation in POINTS, the related   */
 /*     observations in NEWMAP are used to find XCG and YCG for  */
 /*     that county.                                             */
 
data centers(keep=state county flag
                  xsys ysys when function text x y color size);
   length function color $8 text $25;
   retain savptr 1 xold yold 0 when 'a' xsys ysys '2'
          function 'label' size .95 color 'blue' flag 'C';
      set points;                  /* Repeat for each county        */
      xcg=0; ycg=0;                /* xcg,ycg = centers of gravity  */
      aresum=0;                    /* aresum = area of polygon      */
      firstpnt=1;                  /* Flag means first county point */
      endptr=savptr + npoints - 1; /* Set ending pointer for NEWMAP */
      do ptrm=savptr to endptr;    /* Repeat for each county point  */
         set newmap(drop=flag) point=ptrm nobs=nobsm;
         if firstpnt then do;      /* If it is first county point:  */
            xold=x; yold=y;        /* Save first county point       */
            savptr=ptrm + npoints; /* Set next starting pointer     */
            firstpnt=0;            /* Reset flag                    */
            end;    /* xold,yold = previous point in polygon        */
                    /* x, y      = current point in polygon         */
                    /* aretri    = area of triangle formed by last, */
                    /*             previous, and current points     */
         aretri=(xlast - x) * (yold - ylast) +
                (xold - xlast) * (y - ylast);
         xcg + aretri * (x + xold);    /* Update XCG and YCG for    */
         ycg + aretri * (y + yold);    /* each county point         */
         aresum + aretri; /* Accumulate total area of polygon by    */
                          /* summing triangle areas                 */
                          /* (Aresum = double the true area)        */
         xold=x; yold=y;          /* Save previous point            */
         end;
      areinv=1.0 / aresum;        /* Find inverse of total area     */
      x=((xcg * areinv + xlast) * 0.333333); /* Find centers of     */
      y=((ycg * areinv + ylast) * 0.333333); /* of gravity          */
      text=countynm;              /* Text=annotate variable         */

/*  Added this line to subset for only certain counties. */


/*              LOOK      HERE   !!!!!                              */
/*  NC has 100 counties, so it is a very tight fit.                 */
/*  You can subset specific counties, using the following           */
/*  statement....                                                   */
/*    if county in (9 12 13 15 19 21 27);                           */


      output;                             /* Output the observation */
run;
 
 /* 4.  Since NEWMAP and CENTERS contain unprojected coordinates, */
 /*     the points must be projected using PROC GPROJECT.         */

data all;
   set newmap centers;             /* Concatenate NEWMAP to CENTERS */
run;

proc gproject data=all out=allproj;
   id state county;                           /* Project the points */
run;

data mapproj cntrproj;              /* Separate the points into two */
   set allproj;                     /* data sets again              */
   if flag='M' then output mapproj; /* Map data set                 */
   else output cntrproj;            /* Annotate data set            */
run;
 
 /* 5.  PROC GMAP is used to produce the map.  For this applica- */
 /*     tion the same data set (MAPPROJ) is used as both the     */
 /*     response data set and the map data set.                  */

pattern1 v=e c=blue r=100;                /* Number of counties  */
title 'County Map of Nevada';

proc gmap data=mapproj map=mapproj anno=cntrproj;
   id state county;
   choro county / discrete nolegend;
run;
quit;

Operating System and Release Information