BookmarkSubscribeRSS Feed
tmtemples
Fluorite | Level 6

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

5 REPLIES 5
ballardw
Super User

 

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?

 

 

tmtemples
Fluorite | Level 6

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.

ballardw
Super User

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.

tmtemples
Fluorite | Level 6

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!

 

 

 
GraphGuy
Meteorite | Level 14

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

 

 

SAS Innovate 2025: Save the Date

 SAS Innovate 2025 is scheduled for May 6-9 in Orlando, FL. Sign up to be first to learn about the agenda and registration!

Save the date!

How to Concatenate Values

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.

SAS Training: Just a Click Away

 Ready to level-up your skills? Choose your own adventure.

Browse our catalog!

Discussion stats
  • 5 replies
  • 1469 views
  • 0 likes
  • 3 in conversation