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
