Graphing Counts of Physician Types by Zip Code in a Particular state

[HyunJee]

I am trying to develop a map of the state of Oklahoma zip codes that includes color coding to identify how many physician types are located in particular zip codes. I have attached a sample dataset of count of physician types and the corresponding zip codes for the state of Oklahoma. There are some zip codes that are not in Oklahoma, but I planned to exclude them by only keep zip codes where the state = Oklahoma.- I have been having trouble creating the appropriate files to create the exact map I am needing. The below code is what I have written so far, but I know that it is incorrect. I appreciate any help that can be given. Thank you.

proc sort data=zipdata2;

by zip;

run;

data temp (keep = zip count x y city state);

merge zipdata2 sashelp.zipcode;

by zip;

run;

data temp2;

set temp;

if (count ne . and x NE . and y NE .) then do;

      x=atan(1)/45 * x *-1;

      y=atan(1)/45 * y;

      output;

    end;

/* put out a message for bad/missing zipcodes */

else if (count ne . and x=.) then do; 

    put "WARNING: Zipcode " zip " wasn't located.";

    end;

run;

data all;

set maps.states (where=(40)) temp2;

run;

proc gmap data=temp2 map=temp2;

id state zip;

choro count /anno=temp2;

run;

quit;

[GraphGuy]

sashelp.zipcode just has the centroid of the zip codes, not the border/boundary of them.  Therefore you cannot use that data as the map=.  If you want the zip code boundaries, you'll need to download them and import them into a SAS data set that can be used as a map.

http://robslink.com/SAS/democd17/zcta.htm

http://robslink.com/SAS/democd17/zcta_info.htm

If you're going to use the sashelp.zipcodes to annotate points/markers on a map, then you can do that on the SAS county maps, such as ...

http://robslink.com/SAS/democd36/ski_unblinged.htm

http://robslink.com/SAS/democd36/ski_unblinged_info.htm

[HyunJee]

Thank you Robert.

I was unable to open the two last links you provided. The window that pops up says they cannot be found.

The first two links are helpful that you provide. Though I am unsure the code to write to produce the map I need.

I have included the new code below I am using. I should state that I have never created maps in SAS before and this is my first time trying this. I am a bit confused on how to create them and appreciate any help that can be given. Thanks.

*creating a dataset with the correct zip code variables for merging;

data zipdata2;

set zipdata;

length zcta $5.;

zcta = zip_1;

zip = zip_1;

run;

*creating a dataset of the zip code boundaries;

proc mapimport

datafile="./zt40_d00.shp"

out=mymap;

run;

*x/y is really long/lat degrees, so converted them to radians,

so I can gproject them...;

data mymap; set mymap;

  /* Convert long/lat to radians */

  x=atan(1)/45 * x;

  y=atan(1)/45 * y;

run;

*Assigned unique segment values to avoid stray lines between disjoint

3-digit zcta's;

data mymap; set mymap;

segment=(100*ZT40_D00_I)+segment;

run;

*Project the map (note the use of 'eastlong'!);

proc gproject data=mymap out=mymap dupok eastlong project=robinson;

id zcta;

run;

*merging my data with the mymap data;

proc sort data=zipdata2;

by zcta;

run;

proc sort data=mymap;

by zcta;

run;

data combined;

merge zipdata2( in=a) mymap (in=b);

by zcta;

if a and b;

run;

*combining the mymap and my data with te sashelp.zipcode data;

data sas_zipcode;

set sashelp.zipcode;

where state = 40;

length zcta $5.;

zcta = zip;

run;

data combined2;

merge combined sas_zipcode;

by zip;

run;

data mymap2;

merge mymap sas_zipcode;

by zcta;

run;

*creating a graph of the data;

proc gmap data=combined2 map=mymap2;

id state zcta;

choro count;

run;

quit;

[ballardw]

If you are trying to show each type of physician in each ZIP you may want to use the physician type as a BY variable producing one map per physician type.

I have done something similar and unless your final map display is going to be huge a statewide map is unlikely to display ZIP code areas within dense population areas such as Oklahoma City very well. I usually end up at a city level for maps displayed on web sites or in publications. It may also be relatively easy to find a city location SHP file in the same projection as other state features and avoid the projection completely.

[HyunJee]

I am actually trying to show the count of physicians, the type does not matter in this instance. This is a map that is not going to be published, but used a visual tool in presentations given.

I am wondering if anyone can look over my code to see what I am missing or adding that is causing the map to not be developed. I do not have any errors in my log, but no map is produced. Thank you for your help.

[GraphGuy]

I don't think you need to sort your map and then merge it with your zip code data.

A map is a map, and data is data ... you do not need to merge them.

You are probably at the stage where Mike Zdeb's book would be very useful to you:

http://www.amazon.com/Maps-Made-Using-Carpenters-Software/dp/1590470931

[GraphGuy]

The link should be fixed now.

[MikeZdeb · Accepted solution]

Hi ... here's some SAS code that uses an OK ZCTA shapefile and your data to produce some maps.  One is just a ZCTA outline map, the other uses your MD data.  If you've never done a map and you'd like to read a bit (for free, though that's my book that Robert cites), try this ...

The Basics of Map Creation with SAS/GRAPH

There's also a paper about importing shapefiles (though it focuses on census tracts and county subdivisions) ...

From Shapefiles to Maps Via PROC MAPIMPORT and PROC GMAP

SAS code (maps attached) ...

* shapefile from ...
http://legacy.jefferson.kctcs.edu/techcenter/GIS%20data/Country/USA/State/OK/Polygons/ ;

proc mapimport datafile="z:\zt40_d00.shp" out=mymap;
select zcta zt40_d00_i;
rename zcta=zip;
run;

* segment fix;

data mymap;
set mymap;
segment=(100*zt40_d00_i)+segment;
drop zt40_d00_i;
run;

* mymap data set contains unprojected X/Y coordinates ... project them;

proc gproject data=mymap out=zproj degrees eastlong project=robinson;
id zip;
run;

* graphics options;

goptions reset=all ftext='calibri' htext=2 gunit=pct;

* create an outline map of OK ZCTAS;

proc gmap data=zproj map=zproj;
id zip ;
choro segment / nolegend levels=1 coutline=white;
run;
quit;

* change ZIP in your data to a character variable (to match imported data);

data okmds;
set z.zipdata2 (rename=(zip=temp));
zip = put(temp,z5.);
keep zip count;
run;

* title plus white space on sides of map;
title1 h=4 'OKLAHOMA ... PHYSICIANS BY ZIP' ls=2;
title2 a=90 ls=2;
title3 a=-90 ls=2;

* some colors from ... http://colorbrewer2.org/ ;

pattern1 v=ms c=cxD0D1E6;
pattern2 v=ms c=cx74A9CF;
pattern3 v=ms c=cx3690C0;
pattern4 v=ms c=cx0570B0;

* relocate LEGEND ... into Texas panhandle;

legend1 mode=share origin=(10,35)pct across=1 label=(position=top 'PHYSICIANS') shape=bar(3,4)pct;

* some ranges for MD counts;

proc format;
value count
1 = '1'
2-3 = '2-3'
4-5 = '4-5'
6-high = '6+'
;
run;

* map your data;
proc gmap data=okmds map=zproj all;
id zip ;
choro count / discrete coutline=blue legend=legend1;
format count count.;
run;
quit;

[MikeZdeb]

Hi ... one more thing, if you'd like to add some city locations to orient folk using the map ...

* add three cities ... use PROC GEOCODE to find lat/long;
data okcities;
retain state 'OK';
input city & $20. @@;
datalines;
Oklahoma City    Tulsa    Norman
;

proc geocode data=okcities out=okcities_xy (keep=city x y);
run;

* combine unprojected map and city locations;

data map_plus;
set mymap okcities_xy;
run;

* project combined data set;

proc gproject data=map_plus out=map_plus degrees eastlong asis project=robinson;
id zip;
run;

* separate project map from projected city locations;

data
map (where=(city is missing))
cities (where=(city is not missing) keep=city x y)
;
set map_plus;
run;

* create an annotate data set to add cities to map (transparent color ... blue, might not be best choice);

data cities;
retain xsys ysys '2' hsys '3' function 'label' color 'a0000ffaa'  style '"calibri/bo"' size 3 when 'a';
set cities ;
text = upcase(city);
run;

goptions dev=png gsfname=gout xpixels=1024 ypixels=768;

filename gout 'z:\ok.png' ;

* map your data with cities;
proc gmap data=okmds map=map all anno=cities;
id zip;
choro count / discrete coutline=blue legend=legend1;
format count count.;
run;
quit;

[poaform]

hey mike, thanks a lot for the post and the info!