topic Re: Mapping in SAS in Graphics Programming

Mapping in SAS

tmtemples — Fri, 30 Jun 2017 14:47:18 GMT

I am trying to map our data by zip code and overlay the county boundaries. I've been following the code outlined in this document: http://www.misug.org/uploads/8/1/9/1/8191072/rpotter_zip_code_map.pdf.

I got the initial map, but cannot get the county lines to overlay properly. I am new at mapping and greatly appreciate any suggestions!

This initial code gives a map with zip code boundaries and zip code areas shaded accordingly:

/********************READ IN MAP FILES*************************************/
PROC MAPIMPORT DATAFILE = "S:\Infection Control\TinaTemples\SAS_Data\CDIFF\tl_2010_51_zcta510.shp" OUT = VA_ZSHP ;
   ID ZCTA5CE10 ;
RUN ;

DATA VA_ZSHP ;
   SET VA_ZSHP ;
   ZIPCODE = INPUT(ZCTA5CE10, $5.) ;
RUN ;

PROC SORT DATA = VA_ZSHP ; BY ZIPCODE ;
RUN ;

/********************GET SW VIRGINIA ZIP CODES*****************************************/
DATA SW_ZIPS ;
   INFILE "S:\Infection Control\TinaTemples\SAS_Data\Maps\SW_Zips.csv" DSD FIRSTOBS = 2;
   INFORMAT ZIPCODE $ 5. ;

   FORMAT ZIPCODE $ 5. ;

   INPUT ZIPCODE
   ;
RUN ;

PROC SORT DATA = SW_ZIPS ; BY ZIPCODE ;
RUN ;

DATA SW_ZSHP ;
   MERGE SW_ZIPS (IN=SW) VA_ZSHP (IN=VA) ;
   BY ZIPCODE ;
   IF SW AND VA ;
RUN ;

/**********************GET CDIFF CASE ZIPS***************************************************/
DATA ZIP ;
   SET CASE3 (KEEP= ZIP) ;
   FORMAT ZIPCOUNT 1. ;
   ZIPCOUNT = 1 ;
   ZIPCODE = INPUT(ZIP, $5.) ;
RUN ;

PROC SQL ;
   CREATE TABLE ZIP_Y AS
   SELECT DISTINCT SUM(ZIPCOUNT) AS Y, ZIPCODE
   FROM ZIP
   GROUP BY ZIPCODE ;
QUIT ;

TITLE 'C diff Cases in Southwest Virginia by Zipcode' ;
PROC GMAP DATA = ZIP_Y MAP = SW_ZSHP ALL ;
   ID ZIPCODE ;
   CHORO Y ;
RUN ;
QUIT ;

/*************************MAP IS DISTORTED - NEED TO PROJECT*******************************/
PROC GPROJECT DATA = SW_ZSHP DEGREES EASTLONG OUT = SW_ZSHP_PRJ ;
   ID ZIPCODE ;
RUN ;

/*************************MAP THE PROJECTED DATASET****************************************/
TITLE 'C diff Cases in Southwest Virginia by Zipcode' ;
PROC GMAP DATA = ZIP_Y MAP = SW_ZSHP_PRJ ALL ;
   ID ZIPCODE ;
   CHORO Y ;
RUN ;
QUIT ;

TITLE ' ' ;

(copy paste of output did not work)

But when I try the following code, I just get garbage.

/******************CREATE ANNOTATE DATASET*******************************/
DATA SW_VA_ANNO ;
   LENGTH FUNCTION COLOR $ 8. ;
   RETAIN XSYS YSYS '2' COLOR 'BLACK' SIZE 1.75 WHEN 'A' FX FY FUNCTION ;
   SET SW_COUNTIES_PRJ ;
   BY COUNTY SEGMENT ;
   IF FIRST.SEGMENT THEN DO ;
      FUNCTION = 'MOVE' ; *FX = X; *FY = Y ;
   END;
   ELSE IF FUNCTION ^= ' ' THEN DO ;
      IF X = . THEN DO ;
         X = FX ; *Y = FY ; *OUTPUT ; *FUNCTION = ' ' ;
        END ;
       ELSE FUNCTION = 'DRAW' ;
   END ;
   IF FUNCTION ^= ' ' THEN DO ; *OUTPUT ;
      IF LAST.SEGMENT THEN DO ;
         X = FX; *Y = FY ; *OUTPUT ;
       END ;
   END ;
RUN ;

PROC GPROJECT DATA = SW_VA_ANNO OUT = SW_VA_ANNO_PRJ PARALEL1 = 37.053018096 PARALEL2 = 38.078325113 DUPOK ;
   ID COUNTY ;
RUN ;

PROC GMAP DATA = ZIP_Y MAP = SW_ZSHP_PRJ ANNOTATE = SW_VA_ANNO_PRJ ALL ;
   ID ZIPCODE ;
   CHORO Y ;
RUN ;
QUIT ;

/**********************FIX PROJECTION**************************************/
DATA SW_ZSHP_WR ;
   SET SW_ZSHP ;
   X = ATAN(1)/45 * -1 * X ;
   Y = ATAN(1)/45 * Y ;
RUN ;

DATA ANNO_COUNTY_FLAG ;
   SET SW_VA_ANNO ;
   ANNO_FLAG = 1 ;
RUN ;

DATA ZSHP_COUNTYANNO ;
   SET SW_ZSHP_WR ANNO_COUNTY_FLAG ;
RUN ;

PROC GPROJECT DATA = ZSHP_COUNTYANNO RADIANS WESTLONG OUT = ZSHP_COUNTYANNO_PRJ DUPOK ;
   ID ZIPCODE ;
RUN ;

DATA ZSHP_PRJ COUNTYANNO_PRJ ;
   SET ZSHP_COUNTYANNO_PRJ ;
   IF ANNO_FLAG = 1 THEN OUTPUT COUNTYANNO_PRJ ;
      ELSE OUTPUT ZSHP_PRJ ;
RUN ;

PROC GMAP DATA = ZIP_Y MAP = ZSHP_PRJ ANNOTATE = COUNTYANNO_PRJ ALL ;
   ID ZIPCODE ;
   CHORO Y / NOLEGEND ;
RUN ;
QUIT ;

PATTERN ;
TITLE ' ' ;

Thank you for any suggestions!

Tina Temples

Re: Mapping in SAS

ballardw — Fri, 30 Jun 2017 15:01:19 GMT

Did you get to the part of the paper where it shows that the lines won't align and need fixing? And the extra step to align the data?

Re: Mapping in SAS

tmtemples — Fri, 30 Jun 2017 16:20:48 GMT

Yes. Thank you.

I should add that I got an error when I just used the annotate functions, so I added in the paralel piece.

Re: Mapping in SAS

ballardw — Fri, 30 Jun 2017 17:23:03 GMT

Errors => Show the log with the error messages. Best is to paste into a code box opened using the Forum {i} menu icon to preserve formatting of the error messages. The main message window reformats text that can complicate diagnosing errors.

Re: Mapping in SAS

tmtemples — Mon, 03 Jul 2017 17:32:36 GMT

This portion of the log shows the error I receive after the Create Annotate Dataset portion:

296
297
298  /******************CREATE ANNOTATE
298! DATASET*******************************/
299  DATA SW_VA_ANNO ;
300     LENGTH FUNCTION COLOR $ 8. ;
301     RETAIN XSYS YSYS '2' COLOR 'BLACK' SIZE 1.75 WHEN 'A' FX FY
301! FUNCTION ;
302     SET SW_COUNTIES_PRJ ;
303     BY COUNTY SEGMENT ;
304     IF FIRST.SEGMENT THEN DO ;
305        FUNCTION = 'MOVE' ; FX = X; FY = Y ;
306       END;
307     ELSE IF FUNCTION ^= ' ' THEN DO ;
308        IF X = . THEN DO ;
309           X = FX ; Y = FY ; OUTPUT ;  FUNCTION = ' ' ;
310          END ;
311          ELSE FUNCTION = 'DRAW' ;
312       END ;
313     IF FUNCTION ^= ' ' THEN DO ; OUTPUT ;
314        IF LAST.SEGMENT THEN DO ;
315           X = FX; Y = FY ; OUTPUT ;
316          END ;
317      END ;
318  RUN ;

NOTE: There were 4397 observations read from the data set
      WORK.SW_COUNTIES_PRJ.
NOTE: The data set WORK.SW_VA_ANNO has 4146 observations and 14 variables.
NOTE: DATA statement used (Total process time):
      real time           0.01 seconds
      cpu time            0.01 seconds


319
320  PROC GPROJECT DATA = SW_VA_ANNO OUT = SW_VA_ANNO_PRJ DUPOK ; *PARALEL1
320!  = 37.053018096 PARALEL2 = 38.078325113 DUPOK ;
321     ID COUNTY ;
322  RUN ;

NOTE: PARALLEL1 = -0.502555754.
NOTE: PARALLEL2 = 0.5243522619.
ERROR: Standard parallels lie on opposite sides of the equator.
ERROR: Equator is too close to the median latitude of the map data set for
       default parallel calculation.  Try specifying PARALLEL1 and
       PARALLEL2 values explicitly.
NOTE: The SAS System stopped processing this step because of errors.
NOTE: There were 4146 observations read from the data set WORK.SW_VA_ANNO.
WARNING: The data set WORK.SW_VA_ANNO_PRJ may be incomplete.  When this
         step was stopped there were 0 observations and 0 variables.
NOTE: PROCEDURE GPROJECT used (Total process time):
      real time           0.01 seconds
      cpu time            0.01 seconds

My initial map with the zip codes shaded according to frequency looks correct, and my blank map with county lines looks good. I just cannot get the annotate piece to work in order to overlay the county lines. Once that is accomplished, I can use gremove to remove the zip code boundaries.

Thank you for your reply!

Re: Mapping in SAS

GraphGuy — Thu, 26 Apr 2018 19:12:27 GMT

When you project your map, you need to save the projection parameters using the (fairly new) parmout= option. And then when you project your annotated things, use the parmin= option to have gproject use the exact same projection parameters.

Here's an example where I project a map, and annotated bubbles separately, using the parmout= and parmin= options:

http://robslink.com/SAS/democd96/north_korea_missiles_info.htm