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
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?
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.
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.
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!
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
Are you ready for the spotlight? We're accepting content ideas for SAS Innovate 2025 to be held May 6-9 in Orlando, FL. The call is open until September 16. Read more here about why you should contribute and what is in it for you!
Learn how use the CAT functions in SAS to join values from multiple variables into a single value.
Find more tutorials on the SAS Users YouTube channel.
Select SAS Training centers are offering in-person courses. View upcoming courses for: