Shapefile Import and Census Maps

Click for larger image

I was inspired to write this blog entry after reading the article Here are the fastest growing and fastest shrinking counties in the US by Andy Kiersz. Because I am using SAS 9.2 TS2M3, the maps.uscounty data set does not reflect all county FIPS codes thus requiring creation of my own version from a shapefile. The missing county FIPS values can be identified here. Starting in version 9.3, SAS started using Gfk map data sets that are up to date and do reflect the counties missing in traditional map data sets.

The first step in the SAS program (see below) is to use an URL access method to directly read in the county level data from the Census web site. After that, I used proc mapimport to convert a shapefile to a SAS map data set. Because I only wanted to show the 50 states, FIPS defined territories were removed from the imported data.

The imported data is written to separate data sets for the 48 contiguous states, Alaska and Hawaii. After that, each of those data sets are projected for proper display on a map. Alaska is reduced in size by 60% and then both Alaska and Hawaii are repositioned under the southwestern states. Finally all projected data sets are combined together and then sorted by state.

Proc gremove is used to get the state outlines and then use the annotate facility to draw thicker state outlines. Annotation is also used to write values to the center of each state.

A framed shadowed legend is created and positioned to the right of the map output. I utilized a technique learned from SAS/Graph expert, Robert Allison, to create a blank -90 degree angled title to allow space for the legend.

The last step is to create choropleth maps using proc gmap. Both the state and county output use the same color patterns and levels (10). You can see the full PDF version of the output here.

%macro populationchange ;
  %local states statescount i colors colorscount ;
  %let states = states48 alaska hawaii ;
  %let statescount = %sysfunc( countw( &states. ) ) ;
  %let colors = cb0020 df4949 ef9277 f5c0a9 f6e4dd 
                e0ebf1 b3d5e6 82bbd8 4196c7 0471b0 ;
  %let colorscount = %sysfunc( countw( &colors. ) ) ;

  filename census url 
  "http://www.census.gov/popest/data/counties/totals/2014/files/CO-EST2014-alldata.csv" ;
  data statedata countydata ;
    length state county 5  ;
    infile census dsd dlm = "," truncover firstobs = 2 ;
    input sumlev region division state county stname $ ctyname : $32. 
          pop10 estb10 pest10 pet11 pest12 pest13 pest14 ;
    pctchg = ( ( pest14 - pest13 ) / pest13 ) ;
    format pctchg percent8.2 ;
    if county = 0 then output statedata ;
    else output countydata ;
    keep state county pctchg ;
  run ;
  filename census clear ;

  proc mapimport datafile = 'c:\temp\cb\cb_2013_us_county_500k.shp' 
    out = shapefile ;
    select geoid ;
  run ;

  data alaska hawaii states48 ;
    length state county 5 ;
    set shapefile ;
    where segment < 3 and not missing( x ) ;

    state = input( substr( geoid, 1, 2 ), 2. ) ;
    if state > 56 or state in ( 7, 14, 43, 52 ) then delete ;
    county = input( substr( geoid, 3, 3 ), 3. ) ;

    if sign( x ) = 1 then x = -x ;
         if state =  2 then output alaska;
    else if state = 15 then output hawaii ;
    else                    output states48 ;
    drop geoid ;
  run ;

  %do i = 1 %to &statescount. ;
    /* Project each map data set separately. */
    proc gproject data = %scan( &states., &i. ) 
      out  = projected%scan( &states., &i. )
      eastlong deg ;
      id state county ;
    run;
  %end ;

  data projectedalaska ;
    set projectedalaska ;
    x = ( x * .4 ) - .29629586 ;
    y = ( y * .4 ) - .142574262 ;
  run ;

  data projectedhawaii ;
    set projectedhawaii ;
    x = x -.31 ;
    y = y - .120521998 ;
    if x > -.3726 ;
  run ;

  %do i = 2 %to &statescount. ;
    proc append base=projectedstates48 data=projected%scan( &states., &i. ) ;
    run ;
  %end ;

  proc sort data = projectedstates48 out = counties ;
    by state ;
  run ;

  proc gremove data = counties out = states( where = ( segment <= 8 ) ) ;
    id county ;
    by state ;
  run ;

  data countiesoutline ;
    retain xsys ysys '2' when 'a' color 'black' size 2 style 'empty' ;
    length function $8 ;

    set states ;
      by state segment ;
    if first.state or first.segment then function = 'poly' ;
    else function = 'polycont' ;
  run ;

  data mapanno2 ;
    length function $8 text $16 size 8 ;
    retain xsys ysys '2' hsys '3' when 'a' style "'Albany AMT'" ;

    merge
      statedata  
      maps.uscenter( where = ( fipstate( state ) ne 'PR' ) ) ;
        by state ;

    text = catx( " / ", fipstate(state), strip(putn(pctchg,'percent7.2'))) ;
    lagocean = lag( ocean ) ;
    size = 1.5 ;

    if ocean = 'Y' then do ;
      function = 'label' ;
      position = '6' ;
      output ;
      function = 'move' ;
      output ;
    end ;

    position = '5' ;
    if ocean = 'N' then do ;
      if lagocean = 'Y' then do ;
        function = 'draw' ;
        size = round( 1.5 / 4, .01 ) ;
      end ;
      else do ;
        function = 'label' ;
        position = '2' ;
        text = fipstate( state ) ;
        output ;
        position = '5' ;
        text = putn( pctchg, 'percent7.2' )  ;
      end ;
    end ;
    output ;
  run ;

  options nodate nonumber orientation = landscape ;
  %do i = 1 %to &colorscount. ;
    pattern&i. v = s c = cx%scan( &colors., &i. ) ;
  %end ;
  legend1 position = right across = 1 mode = share
    label = ( position = top font = "albany amt/bold" j = c 'Percent Change' ) 
    shape = bar( .12in, .12in ) cshadow = gray 
    offset = ( 2, 3 ) frame ;
 
  title1 h = 3 "2013 to 2014 Population Change";
  title5 a = -90 h = 4pct " ";
  footnote j=l "Source: US Census Bureau, 2014 Population Estimates" ;

  ods listing close ;
  ods pdf file = 'c:\temp\populationchange.pdf' notoc ;
    proc gmap data = statedata 
      map = maps.us( where = ( fipstate( state ) ne 'PR' ) ) all ;
      id state ;
      choro pctchg / anno = mapanno2 legend = legend1 levels = &colorscount. ;
    run ;
    
    proc gmap data = countydata map = counties all ;
      id state county ;
      choro pctchg / anno = countiesoutline legend = legend1  
        levels = &colorscount.  ;
    run ;
    quit ;
  ods pdf close ;
  ods listing ;
%mend ;
%populationchange